Abstract

Since its development in the last decades, controlled radical polymerization (CRP) has become a very promising option for the synthesis of polymers with controlled structure. The design and production of tailor-made materials can be significantly improved by developing models capable of predicting the polymer properties from the operating conditions. Nitroxide-mediated polymerization (NMP) was the first of the three main variants of CRP to be discovered. Although it has lost preference over the years against other CRP alternatives, NMP is still an attractive synthesis method because of its simple experimental implementation and environmental friendliness. This review focuses on deterministic methods employed in mathematical models of NMP. It presents an overview of the different techniques that have been reported for modelling NMP processes in homogeneous and heterogeneous media, covering from the prediction of average properties to the latest techniques for modelling univariate and multivariate distributions of polymer properties. Finally, an outlook of model-based design studies of NMP processes is given.

1. Introduction

Controlled radical polymerization (CRP), also called reversible-deactivation radical polymerization (RDRP) according to the IUPAC recommendation [1], has emerged in the last decades as a very promising option for the synthesis of polymers with controlled structure. In CRP systems, the high reactivity of the polymer radicals is regulated by the addition of an agent which establishes a balance between active and temporarily inactive (dormant) chains. This dynamic equilibrium is strongly shifted towards the dormant chains. As a consequence, the number of active chains is very small which reduces the interactions between them. For this reason, the overall termination effect is decreased and all chains grow on average at the same rate. CRP systems are distinguished by a linear increase of molecular weight with conversion and by a low amount of dead polymer chains (~1–10%).

There are two general approaches to establishing the necessary balance between active and inactive species: the first one is based on a reversible termination (deactivation), while the second one is based on a reversible (degenerative exchange) transfer. In both cases, the polymer radical propagates a few times during activation periods, before is converted back to the inactive state. The initiation rate, the contribution of the termination reactions, and the dynamics of the activation/deactivation exchange influence the polymer dispersity. As the exchange rate increases, which means that fewer monomer units are added to the growing active chains per activation period, the polymer dispersity decreases. Therefore, CRP is effective in controlling molecular weights. Besides, the dynamic equilibrium between active and dormant radicals spans the lifetime of the living polymer chains from fractions of seconds to several hours. This allows manipulating the chain structure during the polymerization reaction. In this way, it is possible to synthesize complex molecular structures, such as comb or star polymers or copolymers with well-defined structure (block, gradient, or branched copolymers) [2, 3].

The three best-known and most effective variants of CRP are atom transfer radical polymerization (ATRP), nitroxide-mediated polymerization (NMP), and reversible addition-fragmentation chain transfer polymerization (RAFT).

1.1. Atom Transfer Radical Polymerization

The ATRP technique is based on the rapid capture of the propagating radicals in a deactivation reaction. Effective ATRP catalysts () are generally organo-metallic species formed by a transition metal capable of changing their oxidation number, a complexing ligand , and a counterion which can form a covalent or ionic bond with the metal centre. The capture of the active radicals is performed by the counterion upon its release after the reduction of the transition metal. Since the activation in ATRP is a bimolecular process, the inactive species formed by the polymer radical and the counterion is inherently stable and can only be activated sporadically, according to the activity of the transition metal that acts as a catalyst. The characteristic kinetic step of ATRP processes is

The reduction of the transition metal complex allows the release of the atom (), which captures the propagating radical forming the dormant species . The success of control on the polymerization is in the rate of periodic activation of the dormant species [4].

1.2. Nitroxide-Mediated Polymerization

In NMP, the reaction is also controlled through the reversible termination of the active chains. Here, the mediating agent is a nitroxide that traps the propagating radical, forming a stable radical. The dissociation step of this reversible termination releases the active radical from the dormant species, and this radical can then propagate or terminate with another radical before being deactivated again. The reaction is strongly displaced towards the inactive species. Therefore, most of the polymer chains in the reaction medium are temporarily dormant, and the concentration of radicals capable of propagating or terminating is several orders of magnitude lower than that of conventional radical polymerization [5]. The distinctive reaction of these polymerizations technique is

Here, species is the nitroxide mediator which, upon trapping the radical, forms the dormant species with monomer units. The slow thermal dissociation of the dormant species provides a low concentration of radicals, which allows the termination step to be kept to a minimum.

1.3. Reversible Addition-Fragmentation Chain Transfer

RAFT polymerization differs from the two above in that control over the growth of the polymer chains is achieved by a degenerative chain transfer process. Basically, what happens is that a single active radical site is shared among many chains so that these molecules participate in fewer reactions. In these systems, as in conventional radical polymerization, a concentration of radicals in the pseudostationary state is established through the initiation and termination processes. Through bimolecular transfer processes, a minimal number of growing radicals undergoes a chain interchange with dormant species through an intermediate adduct. Control over molecular weight and dispersity is provided by chain transfer agents (CTAs) which exchange the active radical centre between the growing chains. Good control requires that the exchange is rapid compared to propagation [6]. The distinctive kinetic step is

The species acts as a chain transfer agent by transferring the radical centre from the chain with monomer units to the chain with units. In an intermediate step, a radical adduct of 2 branches, , is generated. The success of this polymerization is given by a rapid exchange, that is, a large , to keep the concentration of active radicals low [4].

The three main variants of CRP, ATRP, NMP, and RAFT have their own advantages and disadvantages [7]. Although NMP was initially restricted to a limited number of monomers, the development of new nitroxides and alkoxyamines has allowed using this technique with almost all monomers, with exception of vinyl acetate and vinyl chloride. This can be achieved with few nitroxides such as SG1 (and the BlocBuilder MA alkoxyamine) or TIPNO, which are commercially available [5]. Since the development of the different CRP techniques, NMP lost rapidly the preference against RAFT and ATRP because of several advantages of the latter two, such as the capability of polymerizing more monomers and having more convenient reaction temperatures. However, NMP still remains as an attractive synthesis method because of its thermal activation mechanism, a monocomponent control system, a very simple purification of the product, and no environmental issues.

The development of CRP has made it possible to control key polymer properties such as the molecular weight distribution (MWD), the copolymer composition distribution (CCD), or the length distribution of comonomer sequences (SLD). This feature is outstanding for achieving the production of polymers with prespecified properties [8, 9]. Still, the manipulation of the chain structure is a difficult task since the molecular properties have a strong dependence on the operating conditions, which results in complex interactions between process variables. Thus, the design and production of tailor-made materials can be significantly improved by developing models capable of predicting the polymer properties from the operating conditions [10, 11]. The present review focuses on deterministic methods employed in mathematical models of NMP. It presents an overview of the different techniques that have been reported for modelling NMP processes in homogeneous and heterogeneous media, covering from the prediction of average properties to the latest techniques for modelling univariate and multivariate distributions of polymer properties. Finally, an outlook of model-based design studies of NMP processes is given.

2. Deterministic Methods for the Simulation of NMP

Mathematical modelling is undoubtedly a very valuable tool in polymer science. It helps to gain understanding in the kinetic mechanisms of polymer processes and also complements and may partially substitute expensive and time-consuming laboratory experiments. Besides, the model-based design is becoming a new paradigm in polymer processes [12, 13]. Polymer applications are every day more sophisticated and specific. New synthesis methods like CRP allow producing polymers with controlled, complex microstructures. These include the MWD, the CCD, chain functionality, and topology. The polymer chain microstructure determines several properties of the material. Therefore, the combination of experiments and mathematical models allows designing processes for the production of polymers with precise, tailored characteristics [14].

Most of the mathematical techniques used in the simulation of polymer processes can be grouped into two major approaches: deterministic and stochastic methods. Stochastic approaches are mainly represented by the Monte Carlo (MC) method. The majority of the applications of this method use the kinetic Monte Carlo (KMC) technique proposed by Gillespie [15, 16]. A KMC simulation stochastically follows the evolution in time of a full set of molecules of a sample of the reaction system, based on the probabilities of the individual reaction steps of the kinetic mechanism. A KMC simulation tracks explicitly and individually each molecule in the ensemble of chains. Therefore, it can provide extremely detailed information about the polymer microstructure, which is generally not available with deterministic solvers. Besides, the implementation of these methods is simple. On the other hand, due to the statistical nature of the simulation, the reproducibility and reliability of the results depend on the sample size, which needs to be large enough to guarantee accurate outputs. The problem is that large samples require long computational times and large storage capacities. However, the improvement of the computational efficiency of MC implementations is nowadays an area of active research, which includes parallelization, code optimization, and hybrid deterministic-MC approaches [1720]. Recent reviews summarize relevant works using MC methods in polymer science [21, 22]. In NMP, these methods have been used to predict copolymer SLD [10, 11, 23, 24], chain functionality and the full MWD [25], kinetics of branching [26], and bivariate MWD-CCD in continuous processes with arbitrary residence time distributions [27] or to track the exact position of functional comonomers in the copolymer chains [28]. Besides, the model-based design of copolymer synthesis by NMP using MC models has also been performed [11, 28, 29]. The other variants of CRP have also been studied using MC models [20, 3041].

On the other hand, deterministic methods are based on the solution of the differential-algebraic system of population balance equations (PBEs) that represent the polymer system. PBEs are balance equations that describe the evolution in specific properties, called internal coordinates, of the population of polymer molecules. For instance, it is usual to characterize copolymer molecules in a copolymerization reaction by the content of each comonomer in the chain. Therefore, a PBE for this system may involve a differential balance equation for the time evolution of the concentration of species , a copolymer molecule with units of monomer 1 and units of monomer 2 (the internal coordinates). PBEs are derived from a kinetic mechanism of the polymer process. Deterministic methods provide unique results, as opposed to the random output of stochastic methods, and generally demand shorter computing times and less storage capacity. Besides, they are usually more suitable for identification and optimization purposes where it is necessary to count with smooth, differentiable structures.

Deterministic methods have a long history in polymer reaction engineering [42]. In particular, there are articles that review deterministic approaches used in mathematical models of the different variants of CRP [12, 4345]. The present review deepens on deterministic methods employed in mathematical models of NMP. It starts with models describing polymerization rate and average properties only. Then, it discusses techniques for modelling distributions of polymer properties, like the MWD and the SLD. Afterward, it examines complex techniques for predicting multivariate distributions. This is followed by a review of NMP models in heterogeneous media and hybrid deterministic-stochastic methods. Finally, an outlook of model-based design studies of NMP processes is given.

2.1. Average Properties and Polymerization Rate

PBEs constitute a set of equations that is, in principle, infinite in number, because the internal coordinates are distributed properties that theoretically may take values from zero/one to infinity. Since a balance equation must be set for every possible value of the internal coordinates in order to achieve a complete description of the polymer system, an unbounded system of equations is obtained. For instance, in the case of the population, , a mathematical model would need a balance equation for every possible value of and in order to get the full distribution of chain lengths. Although it is sometimes possible to assume a maximum value for a distributed property, the resulting number of equations is still usually large.

Taking advantage that some of their averages can sometimes represent the property distributions, there are techniques that allow limiting the number of model equations. The method of moments, originally proposed by Bamford and Tompa [46], is the most popular method for modelling averages of distributed polymer properties. Considering that polymer species are described only by their chain length, the following equations define the number average molecular weight (), the weight average molecular weight (), and the dispersity () of a homopolymer in a NMP in terms of the leading moments of the MWD:

In these equations, is the monomer molecular weight, and , , and are the moments of order () of the MWD expressed in molar concentration, of propagating radicals, dead polymer, and dormant radicals, respectively. Besides, (7) expresses the monomer conversion () in terms of the first-order moments: where is the monomer molar concentration.

The moments used above are defined as where is the concentration of a polymer species () of chain length .

The method of moments, as applied in the modelling of average molecular weights and conversion, consists in transforming the univariate PBEs of the polymer species (chain length as the only internal coordinate) into balance equations for the moments of their chain length distributions. This procedure replaces the infinite equations of the chain length distribution with just nine equations for the moments of order zero, one, and two of the MWDs of propagating radicals, dead polymer, and dormant radicals. This method is simple and computationally inexpensive [47]. Its drawback, however, is that it loses information on the full shape of the MWD. The process of transforming the PBEs into the moment equations is not complex, but it can be tedious and lengthy for some cases. Recently, Mastan and Zhu [47] presented a tutorial on the use of the method of moments in different polymer systems.

The method of moments has been used extensively to study CRP systems. Examples on the application of this method to ATRP and RAFT can be found elsewhere [12, 4854]. In particular, several works employed the method of moments for the simulation of NMP processes. One of the pioneering works was presented by Veregin et al. [55]. Assuming a kinetic mechanism consisting only of activation/deactivation equilibrium and propagation reactions, as well as constant nitroxide and living radical concentrations, they obtained an analytical solution for the average molecular weights and dispersity. In spite of the simplicity of the kinetic mechanism, they obtained a good match with experimental data. Later, more comprehensive representations were proposed. Butté et al. [56] presented a moment-based model suitable for handling NMP and ATRP processes. The model predicted polymerization rate and average properties and was validated against experimental data. Zhu [50] developed a kinetic model for NMP processes using the method of moments. This model predicted conversion, average molecular weights, and dispersity. He used the model for investigating the effect of operating variables, such as the concentration of initiator, nitroxide, and monomer, on the process performance. Also, he investigated the influence of the rate constants of initiation, propagation, termination, transfer, and equilibrium between propagating and dormant radicals, highlighting the significant effect of the transfer to monomer reaction in the broadening of the MWD. Bonilla et al. [57] presented a moment-based model for the NMP of styrene. They considered a detailed kinetic mechanism that included chemical initiation, reversible nitroxyl ether decomposition, monomer dimerization, thermal initiation, propagation, the equilibrium between propagating and dormant radicals, alkoxyamine decomposition, rate enhancement, transfer to monomer and dimer, and conventional termination. Nondimensional model equations improved the numerical behavior in the parameter estimation. The model predictions showed good agreement with experimental data. Kruse et al. [58] modeled the NMP of styrene predicting average molecular weights and conversion. Using an Evans-Polanyi description of the activation energy (), the model was fit to experimental data of the reaction at 87°C to obtain parameters of the decoupling reaction of dormant radicals and of the propagation/depropagation reactions. They compared the fit with different side reactions, such as transfer to monomer, monomer thermal initiation, transfer to polymer, and the reaction between active radicals and nitroxide to form a hydroxy amine. They found that the monomer thermal initiation was critical for obtaining a good fit. On the other hand, the effect of transfer to polymer and the reaction between active radicals and nitroxide were found to be negligible. Belincanta-Ximenes et al. [59] used the model developed by Bonilla et al. [57] to simulate polymerization rate, average molecular weights, and concentration of polymeric species in the NMP of styrene over a range of operating conditions. They also performed a parameter sensitivity analysis showing the effects of temperature, activation/deactivation equilibrium constant, and initial concentrations of nitroxide and initiator. The simulated profiles of conversion, number average molecular weight, and dispersity agreed well with experimental data. Roa-Luna et al. [60] updated the model by Bonilla et al. [57] including the diffusion-controlled (DC) effects on bimolecular radical termination, propagation, and activation/deactivation reactions. They employed free volume theory to model the DC effects. The authors found that DC termination enhanced the living behavior of the system, whereas DC propagation, DC activation, and DC deactivation worsened it. By comparing model predictions with experimental data, the overall conclusion was that the slight improvement in model performance by the inclusion of the DC effects did not justify the increased model complexity. A later work by this group [61] postulated that side reactions between the initiator and the nitroxide might reduce the effective level of alkoxyamine produced in the early stages of the reaction. Based on the previous model by Bonilla et al. [57], these side reactions were modelled by introducing a controller efficiency factor that determined the fraction of nitroxide able to produce dormant species. This model refinement improved the model predictions, particularly those of the average molecular weights. Asteasuain et al. [62] used the method of moments to predict average molecular properties in the NMP of styrene in tubular reactors. They obtained a good comparison with experimental data.

Copolymerization processes increase the complexity of the mathematical model. In addition to chain length, copolymer composition and the length of the monomer sequences are important properties to be considered. The extension of the method of moments to copolymerization reactions uses double-index moments to compute the average molecular properties. The balance equations for these moments are derived from PBEs that include two internal coordinates, the number of units of each comonomer in the chain. For these bivariate PBEs, the double-index moments of order may be defined as where is the molar concentration of a copolymer species with units of monomer 1 and units of monomer 2.

Copolymer systems require being careful in the definition of the average weights in terms of the molar concentration of the polymeric chains because there is not a straightforward relationship between chain length and mass as in homopolymers [63]. Based on the appropriate definitions, the following equations express the average molecular weights, in an NMP, in terms of the leading double-index moments:

The dispersity is calculated using (6) as previously. The conversion can also be expressed in terms of the double-index moments as

Equations (10), (11), and (12) involve 18 moments, those of order (2,0), (1,1), (0,2), (1,0), (0,1), and (0,0) for each of the three polymeric species. This is twice the number of moments than in homopolymerization models.

The CCD is also an important property of copolymers. Number and weight average values of this distribution can also be expressed in terms of the double-index moments [64]: where and are the number and weight average copolymer composition of monomer 1 in the copolymer, respectively.

Zhang and Ray [65] used this approach in the development of a comprehensive model for CRP in tank reactors. The model was generic and valid for reversible capping CRP systems (i.e., NMP and ATRP). It was validated using experimental data of NMP of styrene and atom transfer radical copolymerization of styrene and n-butylacrylate. Polymerization in batch, semibatch, and a series of continuous tank reactors was analyzed, getting insight into the strengths and weaknesses of the control of the polymer architecture of each operation mode. These authors later extended this work to analyze polymerizations in tubular reactors [66]. Fortunatti et al. [67] used the double-index moment method to predict average molecular weights and copolymer composition in the nitroxide-mediated copolymerization of styrene and α-methyl styrene.

The SLD, that is, the distribution of the length of sequences of each kind of monomer along the polymer chains, is also a significant copolymer property. MC methods are very well suited for giving a detailed insight into this distribution. These methods can predict the explicit sequence of monomers along each copolymer chain in an ensemble of chains [10]. Deterministic models, on the other hand, are in general less exhaustive. They can predict the SLD of the ensemble of chains as a whole, but not of the individual chains. Deterministic models of the SLD are usually set up by formulating a parallel kinetic mechanism that considers the sequences of monomers as the reacting species [68]. For instance, Table 1 shows an extract of a conventional NMP kinetic mechanism, and Table 2 shows the corresponding parallel kinetic mechanism for modelling the SLD [67].

The kinetic mechanism of Table 1 leads to the usual PBEs for , , and , from which the MWD and its averages can be obtained. The parallel kinetic mechanism of Table 2 leads to PBEs for the monomer sequences , , and , which define the copolymer SLD. It should be noted that this approach regards the different monomer sequences as entities disconnected from the polymer chains to which they belong, so no tracking of the monomer sequences along single chains is available. Only one internal coordinate, the sequence chain length, characterizes the monomer sequences. Therefore, number and weight average sequence length, as well as its dispersity, can be computed with equations analogous to (4), (5), and (6). The only difference is that SLD averages are chain length averages, not mass averages, and therefore, the corresponding equations are not multiplied by the monomer molecular weight. Following this approach, Ye and Schork [69] developed a moment-based model that described the MWD and the SLD. This model was able to compute the average properties of the MWD and of the SLD. They compared model data with experimental information. The model was used to describe the conventional and NMP of styrene and 4-methylstyrene. Fortunatti et al. [67] used the same methodology to compute the averages of the copolymer SLD in the nitroxide-mediated copolymerization of styrene and α-methyl styrene.

The computation of the SLD moments provides an alternative to the calculation of copolymer composition, without resorting to double-index moments. Since monomer sequences are composed of a single monomer, the copolymer composition can be computed as the ratio of the 1st-order moment of the SLD of one monomer and the sum of the 1st-order moments of the SLD of both comonomers [69], that is,

Another approach to modelling copolymer reactions is the so-called “pseudo-kinetic rate constants method” [70]. This method defines “pseudo-kinetic” rate constants in terms of the true kinetic constants and the mole fractions of polymer radicals and monomers in the system. The application of the “pseudo-kinetic” constants to the reaction rate expressions of the copolymerization system allows reducing these rate expressions to those of a homopolymerization. Therefore, methods for homopolymerization models, like the one-index moments, can be applied to compute the copolymer average properties. Hernández-Ortiz et al. [71] used this approach to model the nitroxide-mediated copolymerization with cross-linking of vinyl/divinyl monomers. The model assumed that polymer radicals contained only one active radical centre and one dormant radical centre (monoradical assumption). Predictive capabilities included polymerization rate, average molecular weights, gelation point, the evolution of sol and gel weight fractions, cross-link density, and copolymer composition. The model was validated with experimental data of TEMPO-mediated copolymerization of styrene with divinylbenzene. DC effects were assessed and found unimportant. Later [72], they extended their model to consider a more realistic multiradical approach, which took into account the possibility that polymer radicals contained more than one radical centre. In this model, triple-index moments were applied to PBEs whose internal coordinates were the chain length, the number of active centres, and the number of dormant centres. They also obtained a good comparison with experimental data using this model. Recently [73], the group presented a new work in this line that included a different nitroxide and an extended set of experimental information.

Several authors used the commercial software Predici [74] to develop models of NMP systems. Predici is a simulation package for the modelling and dynamic simulation of polymer processes. It employs the Galerkin h-p method to compute the complete MWD without a prior assumption of its shape. Even though the full MWD is available, many works in which Predici was used reported only the prediction of the MWD averages, such as and . These works presented an analysis of reaction steps, in comparison with experimental data and understanding of the influence of operating conditions [7584]. For instance, Gigmes et al. [75] studied the influence of several kinetic parameters and reaction steps on the livingness and controlled character of a generic NMP. Nicolas et al. [76] presented a comprehensive model of the SG1-mediated copolymerization of methyl methacrylate with small amounts of styrene. They used the model to gain insight into the complex mechanism of the copolymerization and to analyze the influence of different operating conditions. Greszta and Matyjaszewski [77] used a Predici model of the NMP of styrene mediated by TEMPO to find the relevant kinetic steps of this polymerization and to estimate the equilibrium constant of the exchange between active and dormant species using experimental data. Benoit et al. [78] developed a model of the NMP of styrene and of n-butyl acrylate using DEPN as a mediating agent. This model showed good agreement with experimental data. Guillaneuf et al. [79] used a Predici model of the NMP of methyl methacrylate mediated by SG1 to show that the penultimate effect was crucial for the accurate description of this system. Fu et al. [80] developed a model of the high-temperature semibatch NMP of styrene. They studied DC effects, showing that the gel effect was not required to improve the fit of the simulation of polymerizations leading to low molecular weight material. Guillaneuf et al. [81] synthesized a new nitroxide, a dialkylated α-hydrogenated linear compound, and the corresponding 1-phenylethyl alkoxyamine, which yielded enhanced polymerization rate in the NMP of styrene. They used a Predici model of the system to study the effect of temperature and alkoxyamine concentration on the kinetics and on the livingness and controlled character of the polymerization. Payne et al. [82] presented a model of the high-temperature NMP of styrene in batch and CSTRs. Using activation/deactivation coefficients estimated from the batch operation, they analyzed the operation in the CSTR finding explanations for experimental information about the molecular weight of the product. Roa-Luna et al. [84] used a model developed in Predici of the NMP of styrene to study the importance of DC effects. Using a conjunction of model results and experimental data, they concluded that DC effects were weak in this system. Hernández-Meza et al. [83] developed a model in Predici of the NMP using hydroxy-TEMPO as mediating agent and initiated by microwave irradiation. Comparing model outputs of average properties with experimental data, they concluded that microwave-activated production of free radicals from monomer molecules was the most probable initiation mechanism.

2.2. Molecular Weight Distribution and Other Distributions of Polymer Properties
2.2.1. Molecular Weight Distribution

Prediction of the full MWD and other distributions of polymer properties is less frequent than the computation of average values because it demands to solve a more complex mathematical problem, the solution of the theoretically infinite set of PBEs. Different approaches have been reported in the literature for modelling univariate, in most cases the MWD, and multivariate distributions of polymer properties. Some of them have been applied to NMP systems. Usually, authors have resorted to some sort of transformation of the PBEs. One of these is the transformation using probability generating function (pgf). The pgf technique [85, 86] is a powerful method that allows modelling the MWD easily and efficiently in complex polymer systems, without any a priori assumption on the MWD shape. PBEs are transformed into the pgf domain, then the transformed balances are solved for the pgf transform of the MWD, and finally, this transform is numerically inverted to recover the MWD. Asteasuain et al. [62] applied this approach to calculate the full MWD in the NMP of styrene in tubular reactors. The model provided a successful comparison against experimental data. Dias and Costa [87] presented a different transformation-inversion approach to model the full MWD, which uses generating functions of the PBEs. They calculated the MWD in linear and nonlinear CRP, including NMP. They showed that under simplified conditions, it was possible to obtain analytical solutions for the MWD. They also presented solutions based on the numerical inversion of the generating functions for more realistic CRP systems.

Other approaches involve the solution of the original set of PBEs. Naturally, different strategies are applied to reduce its infinite dimensionality. First, they all truncate the chain length dimension at a maximum value beyond which the concentration of polymer chains is assumed zero. Then, different strategies are applied to further reduce the size and complexity of the system of equations. Bentein et al. [88] used a method composed by the integration of moment equations and of modified PBEs. The modification involved the application of the quasi-steady-state approximation (QSSA) to the PBEs of macroradicals and a coarse-graining discretization of the dormant and dead polymer balances. In this discretization, the chain length interval was divided into classes, and the properties of the elements of a class were assigned to all the elements of that class. Convergence was checked by comparing the average properties computed using the method of moments and calculated from the MWD. The authors used this approach to get insight into the complex interplay between reaction rates and polymer properties in the SG1- and the TEMPO-mediated bulk polymerization of styrene at 396 K.

Also, Saldívar-Guerra et al. [89] analyzed the calculus of the full MWD by direct integration of PBEs in addition polymerization. They stated that direct integration was feasible in reasonable computing times using standard computers for a considerable number of realistic systems, including NMP. Apart from the truncation of the maximum chain length, they applied the QSSA to remove the stiffness of the system of DAE. Then, explicit solvers, which are computationally less demanding than implicit ones, were used to solve the remaining differential equations. As an example, the NMP of styrene was presented in this work. In a later work by this group, Zapata-González et al. [90] summarized the methodology for the application of the QSSA in CRP and other polymer systems. They illustrated specific features of the methodology by application to academically and industrially relevant systems of controlled radical polymerization (RAFT and NMP cases) and coordination catalysis polymerization.

2.2.2. Sequence Length Distribution

Some authors also addressed the prediction of the full SLD in CRP of copolymer systems using deterministic models [67, 9193]. They all computed the SLD by formulating a model for the sequences of monomers as reacting species, as was discussed above. In the case of the SLD, it is generally feasible to solve the sequences PBEs by direct integration straightforwardly because the maximum sequence length usually takes small values. Regarding NMP systems in particular, Fortunatti et al. [67] used this approach to predict the SLDs of styrene and α-methyl styrene in the copolymerization of these monomers. Aguiar et al. [92] used this method to develop a model for the nitroxide-mediated radical copolymerization of styrene-divinylbenzene. This model included the cyclic propagation (cyclization) reactions. Monomer sequences, in this model, are those that connect a radical centre and a pendant double bond present in the same polymer chain. They achieved good agreement between the model predictions and experimental data for solution- and suspension-controlled copolymerizations.

2.2.3. Multivariate Distributions

Often, an appropriate characterization of a polymer demands simultaneous information on more than one property distribution, for instance, the joint MWD-CCD in copolymer systems or the MWD-branching distributions in branched polymers. The computation of these joint distributions involves the solution of multivariate PBEs, in which more than one internal coordinate is used to describe the polymer molecules. For instance, the prediction of the MWD-CCD requires the formulation of bivariate PBEs whose internal coordinates are the number of units of each comonomer in the chain. The solution of multivariate PBEs is a highly complex problem, for which seldom solution approaches have been developed. Reported methods include the numerical fractionation technique, the fixed pivot technique, or the moment distribution method [9496]. Regarding applications to NMP processes, one of the few reported works is the one by Fortunatti et al. [67]. These authors developed an extension of the pgf method to bivariate distributions [97, 98]. In this method, the bivariate PBEs are transformed using 2D pgfs leading to balance equations for the 2D pgf of the bivariate distribution. After solving this system of equations, the 2D pgf transforms are inverted to obtain the bivariate distribution. In their work [67], Fortunatti et al. calculated the MWD-CCD in the nitroxide-mediated copolymerization of styrene and α-methyl styrene in tubular and semibatch reactors. They analyzed the influence of different feeding policies and temperature profiles on the polymer properties.

2.3. NMP in Heterogeneous Media

CRP techniques have been extensively studied in homogeneous systems. In addition, over the last few years, there has been great interest in adapting these techniques to aqueous heterogeneous systems because they constitute a good alternative for large-scale production. The reasons are many: aqueous dispersed systems are environmentally friendly, show very good heat transfer, better process control, ease of mixing, flexibility, and ease of handling of the final product [99, 100].

CRP in heterogeneous systems presents some complications that have hampered its straightforward dissemination. Among these difficulties, one could mention partitioning of species between aqueous and organic phases, the exit of radicals from polymer particles, poor colloidal stability, and interactions of species with other components of the recipe. The first attempts to adapt CRP to aqueous dispersed systems were carried out in emulsion. They were mostly unsuccessful because the emulsions did not present enough colloidal stability and had mass-transfer limitations. Further studies focused on miniemulsions, where unlike emulsion systems, there is no need of diffusion of reactants through the aqueous phase from monomer droplets to micelles. Miniemulsions are also able to host the polymerization of highly water-insoluble monomers and to form particles containing additives, such as dyes or pigments [100, 101].

Some of the earliest studies on NMP in aqueous dispersed systems were authored by Ma et al. [102104]. They developed two mathematical models that described the nitroxide-mediated miniemulsion polymerization of styrene initiated by alkoxyamine initiators [103] and by the water-soluble initiator potassium persulfate [102, 104]. The models were based on moment equations and included mechanisms describing reactions in the aqueous and organic phases, particle nucleation, the entry and exit of oligomeric radicals, and the partitioning of nitroxide and styrene between the aqueous and organic phases. Compartmentalization effects were neglected, so a droplet-by-droplet approach was not necessary. Instead, they considered a pseudobulk system and formulated the usual PBEs for the organic and aqueous phases. Model predictions included monomer conversion and average molecular weights. The model outputs matched very well the experimental data. In another work by this group [105], Ma et al. studied specifically the interfacial mass transfer of TEMPO in miniemulsion polymerization. For this purpose, they developed a model consisting of two differential equations that described the concentration of TEMPO in the aqueous and organic phases. The model was used to examine how the diffusivity of TEMPO in the aqueous and organic droplet phases, the average droplet diameter, and the nitroxide partition coefficient influenced the time required for the nitroxide to reach phase equilibrium under nonsteady-state conditions.

In order to apply the pseudobulk approximation, it is necessary that the number of chains per particle is relatively large. This approximation works well for many systems [106]. However, if the number of radicals per particle is low, the compartmentalization effect needs to be taken into account. The well-known Smith-Ewart equations [107] have been traditionally used for this situation. The basic Smith-Ewart equation describes the dynamic evolution of the fraction of particles with radicals, , as where , , and are the rate frequencies of radical entry, desorption, and overall termination, respectively. This is a chemical master equation that, unlike rate equations, is valid for small values of the property (i.e.,  = 0,1,2,…). The basic Smith-Ewart equation can be modified to add more information about the radicals within a particle, such as the chain lengths of one (singly distinguished particle) or two (doubly distinguished particle) of them [108, 109]. This allows the formulation of PBEs describing the length of the polymer chains. When formulating the mathematical model, Smith-Ewart equations are used to describe compartmentalized species, while conventional PBEs are used for the noncompartmentalized species, whose concentration within the particles is high enough. The resulting system of equations can be solved or postprocessed using standard techniques to obtain the full MWD or its averages, such as the method of moments and direct integration.

Different models of NMP in miniemulsion employ modified Smith-Ewart equations. Butté et al. [110] were the first to study compartmentalization effects in NMP in miniemulsion. They used Smith-Ewart equations that were modified to account for the number of living radicals and nitroxide molecules within the particle. Charleux [111] also reported a model of NMP of styrene in miniemulsion using SG1 based on the Smith-Ewart equations, but she considered compartmentalization of living radicals only. Zetterlund and coworkers carried out a number of works on NMP in dispersed media. Zetterlund and Okubo [112, 113] used a model based on Smith-Ewart equations that accounted for the number of living radicals and nitroxides in the particles, to investigate compartmentalization effects in NMP of styrene mediated by TEMPO in miniemulsion at 125°C. They studied the effects of particle size and initial nitroxide concentration on the compartmentalization. The main process variables that were predicted were conversion and reaction rates. Their compartmentalization model, however, was not able to compute the MWD or its averages. In a later work [114], they compared the performance of two mediating agents, TEMPO and TIPNO, on the compartmentalization and livingness of the system. They found that the nature of the nitroxide influenced the minimum particle size beyond which compartmentalization could not be neglected. Subsequently, Zetterlund [115] reviewed and also presented new results on compartmentalization and nitroxide partitioning effects. In this work, they implemented a model in Predici to predict the partitioning of the nitroxide between phases. Using the same model, Zetterlund also studied the NMP of butyl acrylate in miniemulsion mediated by TEMPO [116]. He mentioned that TEMPO-mediated acrylate polymerization is known to be problematic and speculated that this is mainly related to the excessive accumulation of free TEMPO with increasing conversion. The author found out that the maximum particle size for which compartmentalization effects were important was larger for this system than for styrene polymerization. He also postulated that there is a narrow window of particle size of about 110–170 nm for which the miniemulsion polymerization of butyl acrylate is faster than in bulk.

Bentein et al. [117] developed a model of the miniemulsion polymerization of styrene mediated by SG1 at 396 K valid up to high conversion which took into account a targeted chain length and particle diameter. They used three-dimensional modified Smith-Ewart equations that counted compartmentalized macroradicals, initiator radicals, and nitroxide radicals. They used the singly distinguished particle approach to developing PBE that allowed using the method of moments to predict average molecular weights and dispersity. Later, Van Steenberge et al. [118] increased in one the dimensions of the Smith-Ewart equations of the model by Beintein et al. [117] to simultaneously calculate the time evolution of the conversion, number average chain length, dispersity, end-group functionality (EGF), and short-chain branching (SCB) content for the miniemulsion NMP of n-butyl acrylate initiated by poly(nBuA)-(N-tert-butyl-N-(1-diethylphosphono-2,2-dimethylpropyl). They concluded that backbiting could not be neglected for an accurate description of the system, despite the low short-chain branching and that the small loss of end-group functionality at conversions was mainly caused by chain transfer to monomer.

García-Leal et al. [119] studied the NMP of styrene using TEMPO for the production of high molecular weight polystyrene in a mass-suspension process. The main characteristic of the mass-suspension polymerization is that a bulk prepolymerization is first carried out, followed by the suspension of the prepolymer at intermediate conversion in a continuous phase until total conversion. This reduces the inter- and intraparticle mass transfer that occurs at low to moderate conversion in a classic suspension polymerization and permits a better control of the particle size distribution [119]. The authors assumed that the system kinetics behaved as in a conventional polymerization and used the classical PBEs of NMP processed by the method of moments to predict average molecular properties and conversion.

2.4. Hybrid Methods

A KMC simulation provides extremely detailed information of the polymer microstructure because it tracks the evolution in time of a full set of molecules of a sample of the reaction system. However, it usually requires long computational times and large storage capacities. Therefore, a number of authors have been working on the development of hybrid stochastic-deterministic solvers [17, 120, 121]. These solvers try to substitute parts of the stochastic simulation by deterministic techniques in order to accelerate the computation time. Regarding applications to NMP systems, Özçam and Teymour [17] developed a hybrid deterministic-probabilistic method that builds chains one-by-one or chain-by-chain and was named “Chain-by-Chain Monte Carlo” method. A deterministic component of the algorithm, based on the method of moments, is used to obtain information about the global concentration of polymer species and of small molecules, such as monomer or initiator. This information is passed to the MC part of the algorithm. They tested this hybrid model on the synthesis of styrene/methyl methacrylate linear gradient copolymers via NMP and methyl methacrylate/methyl acrylate linear hyperbolic gradient copolymers via ATRP. Results were compared with a conventional KMC model and with a deterministic method-of-moments model and confirmed that full information regarding the microstructure of chains can be obtained using the hybrid method with reduced simulation times and smaller sample sizes.

2.5. Model-Based Design of NMP Processes

Model-based design plays a fundamental role in polymer reaction engineering. Polymers are “process products”: the microstructure of the material, and therefore, its application properties are strongly dictated by the conditions of the manufacturing process. Model-based design allows for a faster exploration of the search space, reduced uncertainty, and better, faster, and safer design and operation decisions through a deeper understanding of processes. The controlled growth of polymer chains in CRP processes makes them particularly suitable for synthesizing polymers with tailor-made molecular architecture. Model-based design is a very helpful tool in obtaining the specific operating policies required to achieve the desired molecular structure [14]. Deterministic modelling methods are particularly suitable for model-based design applications due to their fast running times, in comparison with stochastic approaches, and the availability of derivatives for gradient-based optimization solvers. Nevertheless, the model-based design of NMP processes using MC models has also been performed [11, 28, 29].

Different authors performed model-based design studies of NMP processes using deterministic models. Faliks et al. [122] applied an optimal control methodology to the optimization of an NMP of styrene in a plug flow reactor. The main purpose of the work was to show the benefits of optimal control in NMP reaction engineering, for which a typical case study was selected. The optimization problem aimed at finding the best nitroxide radical flux along the length of the reactor, with the objective of matching a target chain length and keeping conversion high. The mathematical model of the process was based on moment equations and predicted conversion and average molecular properties. The authors found that it was possible to reduce the reaction time significantly for a given conversion. However, this came at the expense of a modest increase in dispersity.

Lemoine-Nava et al. [123] used a moment-based model to find optimal operating policies for a simulated industrial reactor of NMP of styrene. They posed a dynamic optimization problem to find time-dependent optimal operating policies that maximized conversion and minimized dispersity.

Asteasuain et al. [124] used a model of the NMP of styrene in tubular reactors based on the pgf technique [62], for the optimization of the process aiming at producing tailor-made MWDs. Optimization variables included the feed rates of styrene, initiator (BPO), and TEMPO at the reactor inlet and lateral feeds, the location of the lateral feeds, and the operating temperature. They showed that it was possible to obtain bimodal and trimodal MWDs, with predetermined peak locations, through proper manipulation of the operating conditions. The authors highlighted that the pgf technique used to model the full MWD presented the advantage that it allowed tailoring the size of the mathematical model as a function of the particular details of the MWD that would be specified in the optimization problem. This permitted to reduce significantly the running time of the optimization tasks.

Zitlalpopoca-Soriano et al. [125] optimized grade transitions in the NMP of styrene in tubular reactors. They investigated the effect of the reactor temperature, monomer feed rate, and cooling water flow rate as manipulated variables for the grade transition. The reactor model involved moment equations that could predict average molecular properties and conversion. The DAE optimization problem was solved using a simultaneous approach wherein the differential and the algebraic variables were fully discretized leading to a large-scale nonlinear programming (NLP) problem. The resulting optimization problem was solved using an interior point algorithm capable of handling large-scale NLPs.

When it comes to process design, optimization objectives are often conflictive. For instance, polymers synthesized by CRP (e.g., NMP) are usually required to have the narrowest possible MWD and highest end-group functionality (EGF). This goal is usually contradictory to the economic objective of minimizing the reaction time for a given conversion. These conflicting objectives are best dealt with by means of multiobjective optimization (MOO) techniques. Recently, Fierens et al. [126] showed the benefits of applying MOO techniques to NMP of styrene initiated by MAMA-SG1 in a semibatch reactor. For this purpose, they coupled a moment-based model of the process with the nondominated sorting genetic algorithm-II (NSGA-II) in FORTRAN code to obtain the Pareto set of optimal solutions. The Pareto set is a collection of optimal solutions in which none of the objective functions can be improved in value without degrading some of the other ones. The design objectives of this study were the minimization of reaction time and of the dispersity and the maximization of the EGF. Process optimization variables included the temperature path, initial reactor charge, and semibatch feeding flows.

Bifurcation analysis is a widely known technique for studying the operability of chemical processes. It helps to assess steady-state multiplicity in chemical reactors and identify stable and unstable operating regions. A few works reported the use of this technique in NMP processes. Lemoine-Nava et al. [127] studied the nonlinear behavior of the NMP of styrene in a CSTR. They selected the cooling water flow rate, feed-stream temperature, cooling water feed temperature, monomer feed-stream concentration, and residence time as the bifurcation parameters. They found input multiplicities, disjoint bifurcations, and isola and hysteresis behaviors, which made the reactor prone to operability and control problems. Later, Zitlalpopoca-Soriano et al. [128] extended this analysis to the NMP of styrene in a tubular reactor. They studied the effect of the cooling flow rate, the feed-stream temperature, the cooling water temperature, the feed-stream monomer concentration, the flow rate of monomer feed, diameter, and global heat transfer. They used moment-based models [57, 125] in both works.

Another very useful application of mathematical models to the process development is the model-based design of experiments. This technique uses process models to design experiments in an optimal way in order to maximize the information that can be extracted from them and minimize the experimental effort required to estimate model parameters. Nabifar et al. [129] demonstrated the benefits of the Bayesian design methodology in complex polymerization processes using an NMP system as a case study. Scott et al. [130] carried out a model-based D-optimal design of experiments in the NMP of styrene. They used a full, moment-based model developed previously by the group [57], as well as a refined model which included new model parameters and revised kinetic steps, and a reduced model for faster simulation runs.

3. Conclusions

Mathematical modelling is a very helpful tool in polymer reaction engineering. It helps to gain understanding in the kinetic mechanisms of polymer processes, complements expensive and time-consuming laboratory experiments, and aids in obtaining the specific operating policies required to achieve the desired molecular structure. This applies especially to CRP, such as NMP, because these synthesis routes are particularly meant to manufacture polymers with controlled, precise molecular structures. This manuscript presents a review that deepens on the main deterministic modelling techniques employed for the simulation of NMP processes. It is described that most of the reported models focus on the prediction of average properties, such as average molecular weights and conversion. The method of moments is an overwhelming most selected option for this task, although commercial packages like Predici have also been employed. Other works applied more sophisticated techniques for modelling distributions of polymer properties, including recent developments for the calculation of multivariate distributions. A description of mathematical models for the simulation of NMP in heterogeneous systems is then presented, which shows that Smith-Ewart equations have been most used for analyzing compartmentalization effects in living radicals and dormant species. Finally, the manuscript gives an outlook of model-based design studies of NMP processes.

Conflicts of Interest

The author declares that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

The author wishes to thank the National Research Council of Argentina (CONICET) and the Universidad Nacional del Sur (UNS) for the financial support.