Reactor Modeling and Recipe Optimization of Ring-Opening

Oct 23, 2013 - ABSTRACT: This paper addresses reactor modeling and recipe optimization of semibatch ring-opening polymerization processes for...
0 downloads 0 Views 1MB Size
Article pubs.acs.org/IECR

Reactor Modeling and Recipe Optimization of Ring-Opening Polymerization: Block Copolymers Yisu Nie,† Lorenz T. Biegler,*,† Carlos M. Villa,‡ and John M. Wassick¶ †

Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, United States The Dow Chemical Company, Freeport, Texas 77541, United States ¶ The Dow Chemical Company, Midland, Michigan 48674, United States ‡

ABSTRACT: This paper addresses reactor modeling and recipe optimization of semibatch ring-opening polymerization processes for making block copolymers. Two rigorous reactor models are developed on the basis of the population balance and method of moments, respectively. The complete polymerization process model also includes vapor−liquid equilibrium equations from applying Flory−Huggins theory. The accuracies of both reactor models are validated against historical plant data by adjusting model parameters such as kinetic rate constants. The recipe optimization problem is formulated to design the optimal reactor operating policy to minimize polymerization time and incorporate additional process constraints in accordance with final product properties and process safety requirements. The resulting dynamic optimization problem is translated to a nonlinear program by using the simultaneous collocation method, and further solved by the interior point method. In the case study example, both reactor models show satisfactory matches between their predictions and the historical plant data. The recipe optimization with both models demonstrates significant process improvement and reductions in batch operating time. Moreover, the moment model shows superiority over the population balance model in terms of computational efficiency.



INTRODUCTION A great variety of cyclic monomers have been successfully polymerized by the ring-opening polymerization process.1 Many commercial polymers produced this way are statistical or block copolymers with molecular weights (MWs) from a few hundred to several million. Conventionally, a block copolymer can be synthesized by a semibatch polymerization process by successive feeds of respective monomers.2 For instance, a block copolymer made from monomer A and B is shown in Figure 1 with a cascade structure. The initiator has two branches, and each branch consists of a linear internal A block and external B block in tandem. The external block may contain a small amount of A in practice. Modeling such a polymerization process requires sufficient knowledge of the reaction kinetics, and process safety regulations are also important considerations for process optimization.3 In this example, monomer B is potentially explosive such that the reactor for making this polymer needs to be blanketed with nitrogen at the start. Population balance equations are often used in macro-scale polymerization reactor modeling, which correspond to a set of differential mole balance equations that can reveal the evolution of the complete chain length distribution over time. The size of the resulting balance equation system is decided by the breadth of the polymer chain length distribution, where modeling high MW polymers can lead to a considerably large-scale model that is computationally demanding. On the contrary, the method of moments,4 as a classical modeling approach, represents the average polymer properties based on statistical quantities, namely, moments. The moments can be applied to derive commonly used quality indices such as number/weight average molecular weight, polydispersity index (PDI), and monomer consumption/polymer production rates. The size of moment models is usually small and is not influenced by the chain length of polymers. As a result, solving a moment model requires much less computational effort than the full population balance model. In this study, we derive both population © 2013 American Chemical Society

balance and moment models; the population balance model is appropriate for process recipe simulation, while the moment model facilitates recipe design optimization.5−8 Our previous study9 investigates the optimal design of a homopolymerization process by solving dynamic optimization problems with population balance models. In the dynamic optimization formulation, the participating constraints include the reactor model as a differential-algebraic equation (DAE) system, and process constraints described via algebraic equalities or inequalities. The dynamic optimization formulation is translated into a continuous nonlinear program (NLP) by applying the simultaneous collocation method, where both state and control variables of a DAE system are discretized using orthogonal collocation on finite elements (OCFE). This discretization scheme corresponds to an implicit Runge−Kutta method with high order accuracy and excellent stability properties; and also allows direct enforcement of state and control variable constraints.10 The translated NLP problem is solved with the interior point optimizer (IPOPT) which implements the barrier approach coupled with filter line-search methods.11 In this study, we follow the same methodology in dealing with the recipe optimization task. In the next section, a detailed reactor model based on population balances is first derived, followed by a reformulation procedure aiming to reduce the stiffness of the resulting DAE system. Next, the corresponding moment model is presented. Then we discuss the recipe optimization problem and its solution algorithm with a case study example of the block polymer in Special Issue: John Congalidis Memorial Received: Revised: Accepted: Published: 7434

August 22, 2013 October 14, 2013 October 23, 2013 October 23, 2013 dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446

Industrial & Engineering Chemistry Research

Article

Figure 1. General structure of a block copolymer.

Figure 1. Our models are first calibrated against historical plant data for validation, and both models, based on the population balances and moments, are optimized. Lastly, we give concluding remarks and future directions.

To summarize the reaction scheme, the following notation is introduced. Let M denote the monomer, PnA,nB denote a linear branch of the product polymer chains and UnA,nB represents the byproduct chains. Depending on the different functional end groups, we define the following.



REACTOR MODEL DEVELOPMENT An anionic ring-opening polymerization process can be carried out in conventional stirred-tank reactors coupled with heat exchangers. These polymerization processes are often catalyzed by alkali metal hydroxides.12 Molecules with active hydrogen atoms are used as the initiation site (starter) for polymerization. Initially, the polymerization starter and catalyst are mixed in the reactor tank in an appropriate ratio. To produce monomer blocks, monomer A is first fed into the reactor continuously to grow the internal block. Next, monomer B is allowed to enter the reactor to form the external block. A degassing step can be performed to eliminate the unreacted A by vacuum distillation before feeding B to obtain a purer external block. During polymerization, external heat is usually required in the start-up stage, and soon after the polymerization reactions have been kicked off a significant amount of heat is released from the reactions and needs to be removed from the tank. In addition, a venting control system is installed on the reactor to prevent extreme reactor pressures. Reaction Mechanism. Two stages are defined for the block copolymerization process. In the first stage of polymerizing A, four primary reactions are considered: the initiation, propagation, transfer, and exchange reactions. The initiation reaction is carried out by ring-opening of cyclic monomer A, and the formed chains are extended with monomer units. The reaction rate of the initiation reaction is in general higher than the propagation rate. The transfer reaction creates undesired byproduct molecules (impurities) which undertake the same set of reactions as the product polymer chains. The side reaction is favored by higher temperatures. Impurity levels may increase with the number average molecular weight of the polymer in some circumstances, which suggests that the rate of polymerization decreases relative to the rate of chain transfer. During anionic polymerization, the growing polymer chains ending with the alkali metal cation are in equilibrium with the dormant polymer chains ending with the hydrogen atoms. These acid−base proton exchange reactions are the fastest reactions in the polymerization scheme, and the equilibrium constants are around unity since the acidity of the participating species is similar. In the secondary step, B is added to form terminal blocks by the corresponding initiation and propagation steps with higher reactivity than A. However, B does not participate in transfer reactions. If A is not completely digested or degassed, then A and B coexist in the reactor in the second stage. Under this circumstance, there are four propagation reactions, shown as follows: propagation scheme rate constant ...A* + A → ...A* kAA p ...A* + B → ...B* kAB p ...B* + A → ...A* kBA p ...B* + B → ...B* kBB p Two associated reactivity ratios are defined as r BA r = kBB p /kp . B

A

GnA,nB to denote the growing product chains DnA,nB to denote the dormant product chains QnA,nB to denote the growing byproduct chains RnA,nB to denote the dormant byproduct chains Moreover, superscripts are introduced to the polymer species to indicate the terminal repeating unit or monomer type, except for the initiators that do not have any repeating units (nA = nB = 0). The reaction scheme is summarized in Table 1. We assume the byproduct chains share the same set of kinetic parameters as the product chains in chain initiation, growth, exchange, and transfer. The four propagation rates are nonidentical. All the exchange reactions are reversible equilibrium reactions, and some of them are denoted with a single reaction rate because of the symmetry between reactants and products. First-Principles Model Equations. According to the reaction scheme described earlier, the population balance equations for individual polymeric species can be established as a set of first-order ordinary differential equations, with generation and consumption rates on the right-hand side: dG0,0 dt

= V −1[−(k iAM A + k iBMB)G0,0 − ktG0,0M A ∞

− keG0,0



∑ ∑

(DmA , mB + R mA , mB)

mA = 0 mB = 0 ∞ ∞

+ keD0,0

∑ ∑

(GmA , mB + Q m

A dG1,0

dt

)]

A , mB

mA = 0 mB = 0

(1a)

A = V −1[k iAG0,0M A − (k pAAM A + k pABMB)G1,0 ∞



A k tG1,0 MA



A keG1,0



∑ ∑

(DmA , mB + R mA , mB)

mA = 0 mB = 0 ∞ A + keD1,0



∑ ∑

(GmA , mB + Q m

B dG0,1

dt

)]

A , mB

mA = 0 mB = 0

(1b)

B = V −1[k iBG0,0MB − (k pBAM A + k pBBMB)G1,0 ∞



B k tG0,1 MA



B keG1,0



∑ ∑

(DmA , mB + R mA , mB)

mA = 0 mB = 0 ∞ B + keD1,0

AB = kAA p /kp and



∑ ∑ mA = 0 mB = 0

7435

(GmA , mB + Q m

)]

A , mB

(1c)

dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446

Industrial & Engineering Chemistry Research dGnAA , nB dt

Article B dQ 0,1

= V −1[(k pAAGnAA − 1, nB + k pBAGnBA − 1, nB)M A

B = V −1[k iBQ 0,0MB − (k pBAM A + k pBBMB)Q 0,1

dt



− (k pAAM A + k pABMB)GnAA , nB − k tGnAA , nBM A ∞



∑ ∑

− keGnAA , nB

B B M A − keQ 0,1 − ktQ 0,1

∑ ∑



(DmA , mB + R mA , mB)

B + keR 0,1



∑ ∑

(GmA , mB + Q m

(GmA , mB + Q m

A

)], ,m

dQ nA , n

B

A

nA ≥ 2, nB ≥ 0

A

B

A

dt

=V

[(k pABGnAA , nB − 1

+



A

+

∑ ∑

A

)],

+





∑ ∑

B−1

+ k pBBQ nB , n B

A

dR nSA , nB

)],

dt

A , mB

mA = 0 mB = 0

B



∑ ∑

(DmA , mB + R mA , mB)

mA = 0 mB = 0 ∞ ∞

∑ ∑

(GmA , mB + Q m

)],

A , mB

nA ≥ 0, nB ≥ 2

(DmA , mB + R mA , mB) (GmA , mB + Q m

)M B

B−1

A

mA = 0 mB = 0

mA = 0 mB = 0 ∞ ∞

− keDnSA , nB

B

+ keR nBA , nB

= V −1[ktGnSA , nBM A

∑ ∑

(1j)

A

(1e)



)],

A , mB

− (k pBAM A + k pBBMB)Q nB , n − ktQ nB , n M A A

keGnSA , nB

(GmA , mB + Q m

A

− keQ nB , n

dt

(DmA , mB + R mA , mB)

= V −1[(k pABQ nA , n

B

dt

A , mB

nA ≥ 0, nB ≥ 2

dDnSA , nB

B

nA ≥ 2, nB ≥ 0 dQ nB , n

mA = 0 mB = 0

A

mA = 0 mB = 0

(DmA , mB + R mA , mB) (GmA , mB + Q m

B

mA = 0 mB = 0 ∞ ∞

∑ ∑

+ keR nAA , nB

mA = 0 mB = 0 ∞ ∞

keDnBA , nB

B



∑ ∑

B



∑ ∑

− keQ nA , n

k pBBGnBA , nB − 1)MB

− (k pBAM A + k pBBMB)GnBA , nB − ktGnBA , nBM A − keGnBA , nB

A

− (k pAAM A + k pABMB)Q nA , n − ktQ nA , n M A ∞

−1

∞ A B

− keR nSA , nB

(1f)

(1k)

= V −1[ktQ nS , n M A + keQ nS , n ∞

S ∈ {A, B}, nA ≥ 0, nB ≥ 0

A B



∑ ∑

dt



∑ ∑

(GmA , mB + Q m

)],

A , mB

mA = 0 mB = 0

+ kt



∑ ∑

− keQ 0,0

A

(GnA , nB + Q n

)M − ktQ 0,0M

xnA , nB =

A

A , nB

nA = 0 nB = 0 ∞

∑ ∑

(GmA , mB + Q m

)]

A , mB

mA = 0 mB = 0

(2)

The monomer balance equations are defined for both monomers. Monomer A is fed into the reactor at rate FA, and consumed in the initiation, propagation, and chain transfer reactions resulting in byproduct chains. The balance equation is similar for monomer B, except for the different initiation and propagation rates and absence of the proton transfer reaction term.

(DmA , mB + R mA , mB)

mA = 0 mB = 0 ∞ ∞

+ keR 0,0

xnSA , nB , x ∈ {G , D , Q , R }

∑ S ∈ {A,B}



∑ ∑

(1l)

Here the polymer species without a superscript index refer to the total population, regardless of the different terminal units:

= V −1[−(k iAM A + k iBMB)Q 0,0 ∞

(DmA , mB + R mA , mB)

mA = 0 mB = 0

S ∈ {A, B}, nA ≥ 0, nB ≥ 0

dQ 0,0

(1i)

= V −1[(k pAAQ nA − 1, n + k pBAQ nB − 1, n )M A

B

dt

(1d)

)]

A , mB

mA = 0 mB = 0

mA = 0 mB = 0

dGnBA , nB

(DmA , mB + R mA , mB)

mA = 0 mB = 0

mA = 0 mB = 0 ∞ ∞

+ keDnAA , nB



∑ ∑

(1g)

dM A = F A − V −1[k iA(G0,0 + Q 0,0) dt A dQ 1.0

dt



=V

−1

[k iAQ 0,0M A



(k pAAM A ∞

A A A M − keQ 1,0 − ktQ 1,0

+

+ k pAA

A k pABMB)Q 1,0

(DmA , mB + R mA , mB)

+ k pBA

mA = 0 mB = 0 ∞ A + keR1,0

mA = 0 mB = 0

A

∑ ∑

(GmA , mB + Q m

)]

A , mB

+ kt (1h)

∑ ∑ nA = 0 nB = 0

7436

B

(GnBA , nB + Q nB , n ) A

nA = 0 nB = 1 ∞ ∞



∑ ∑

(GnAA , nB + Q nA , n )

nA = 1 nB = 0 ∞ ∞



∑ ∑



∑ ∑

(GnA , nB + Q n

B

)]M A

A , nB

(3a)

dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446

Industrial & Engineering Chemistry Research

Article

dM B = F B − V −1[k iB(G0,0 + Q 0,0) dt ∞

+ k pAB



∑ ∑

Initiation:

(GnAA , nB + Q nA , n ) A

nA = 1 nB = 0 ∞ ∞

+ k pBB

Table 1. Reactions in Anionic Ring-Opening Copolymerization

∑ ∑ (GnB ,n A

nA = 0 nB = 1

B

k iA

B

A G0,0 + M A ⎯→ ⎯ G1,0 k iA

+ Q nB , n )]MB A

B

A Q 0,0 + M A ⎯→ ⎯ Q 1,0

(3b)

k iB

B G0,0 + MB → G0,1

For the total mass balance, as the monomers consecutively enter the polymerization system, we have: dm = F AM WA + F BM WB dt

k iB

B Q 0,0 + MB → Q 0,1

Propagation:

(4)

k pAA

GnAA , nB + M A ⎯⎯⎯→ GnAA + 1, nB

where MWA and MWB denote the molecular weights of the two monomers, respectively; and m is the total mass of the polymerization system. The liquid density is solely dependent on the reactor temperature since the effects of molecular weight on density are found to be minor. The liquid volume is calculated by:9 V = m[10−6 + 7.576 × 10−10(T − 298.15)]

Q nA , n + M ⎯⎯⎯→ Q nA + 1, n A B

∑ Pi i

Pi = aiPisat· Bi = Ai − T + Ci

Q nA , n + M ⎯⎯⎯→ Q nB , n

(5)

⎛ 1⎞ ln a sB = ln ϕs + ⎜1 − ⎟ϕB + χ B ϕB2 B ⎝ l⎠

A B+1

A B

k pBA

GnBA , nB + M A ⎯⎯⎯→ GnAA + 1, nB k pBA

A

B

k pBB

GnBA , nB + MB ⎯⎯⎯→ GnBA , nB + 1

(6a)

k BB B p

Q nB , n + M ⎯⎯⎯→ Q nB , n

A B+1

A B

(6b)

(nA + nB ≥ 1)

(nA + nB ≥ 1)

Q nB , n + M A ⎯⎯⎯→ Q nA + 1, n A B

(nA + nB ≥ 1) (nA + nB ≥ 1)

k AB B p

(nA + nB ≥ 1)

(nA + nB ≥ 1) (nA + nB ≥ 1)

Transfer: kt

GnAA , nB + M A → DnAA , nB + Q 0,0

(nA , nB ≥ 0)

kt

GnBA , nB + M A → DnBA , nB + Q 0,0

(7)

(nA , nB ≥ 0)

A kt

The activities ai are given by the Flory−Huggins theory.13 In the polymerization process, possible volatile components include the monomer A, B, initiators, etc. As a result, multiple components exist in the liquid and gaseous phase simultaneously. However, as a rational method to deal with polymer−solvent equilibrium, the Flory−Huggins theory is in principle applicable to binary mixtures of a polymer and solvent. To manage the VLE calculation, we make the following assumptions. First, the polymer chains are totally in the liquid phase; second the polymer solution can be viewed as a ternary mixture where A and B (denoted in index sA and sB, respectively) are two types of solvents and all the other species are treated as polymer (denoted in index p); last, the interaction between the two solvent molecules is ignored. Hence, the activities of A and B are determined via: ⎛ 1⎞ ln a sA = ln ϕs + ⎜1 − ⎟ϕA + χ A ϕA 2 A ⎝ l⎠

B

GnAA , nB + M ⎯⎯⎯→ GnBA , nB + 1

The vapor pressures by the Antoine equation: log10 Pisat

A

k AB B p

Finally, a set of VLE equations is incorporated as a means to predict the reactor pressure. The total reactor pressure P is calculated as the sum of the partial pressures of volatile components, determined from the liquid phase activities and saturated vapor pressures:

P=

(nA + nB ≥ 1)

k AA A p

Q nA , n + M → R nAA , nB + Q 0,0

(nA , nB ≥ 0)

A B

A kt

Q nB , n + M → R nBA , nB + Q 0,0

(nA , nB ≥ 0)

A B

Exchange: ke

GnAA , nB + DmAA , mB → DnAA , nB + GmAA , mB

(nA , nB , mA , mB ≥ 0)

ke

GnAA , nB + DmBA , mB ⇌ DnAA , nB + GmBA , mB

(nA , nB , mA , mB ≥ 0)

ke

ke

GnBA , nB + DmBA , mB → DnBA , nB + GmBA , mB

(nA , nB , mA , mB ≥ 0)

ke

Q nA , n + R mAA , mB → R nAA , nB + Q mA , m A B

B

(nA , nB , mA , mB ≥ 0)

A , mB

(nA , nB , mA , mB ≥ 0)

A

ke

Q nA , n + R mBA , mB ⇌ R nAA , nB + Q mB A B

ke ke

Q nB , n + R mBA , mB → R nBA , nB + Q mB

A , mB

A B

(8a)

ke

GnAA , nB + R mAA , mB ⇌ DnAA , nB + Q mA , m ke

B

(nA , nB , mA , mB ≥ 0)

A , mB

(nA , nB , mA , mB ≥ 0)

A

ke

GnAA , nB + R mBA , mB ⇌ DnAA , nB + Q mB

(8b)

ke

In eq (8), the lattice fractions ϕSA, ϕSB and ϕp are used instead of mole fractions. The interaction parameters χA and χB are nondimensional and account for the energy of interdispersing polymer and solvent molecules. Here the interaction effect of monomer A and B molecules is not considered for simplicity as well as the fact that a successive feeding pattern is used for block polymers. A more rigorous treatment is proposed by Favre et al.14 at the cost of additional complexity in activity calculation with more model

ke

GnBA , nB + R mAA , mB ⇌ DnBA , nB + Q mA , m ke

B

(nA , nB , mA , mB ≥ 0)

A , mB

(nA , nB , mA , mB ≥ 0)

A

ke

GnBA , nB + R mBA , mB ⇌ DnBA , nB + Q mB ke

(nA , nB , mA , mB ≥ 0)

parameters. The number average chain length l is taken into account when calculating the lattice fraction: 7437

dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446

Industrial & Engineering Chemistry Research

Article

n sA

ϕs =

dX 0,0

n sA + n sB + npl

A

n sB

ϕs =

A dX1,0

n sA + n sB + npl

B

= V −1( −k iAM A − k iBMB)G0,0

dt

(9a)

A ] = V −1[k iAG0,0M A − (k pAAM A + k pABMB)G1,0

dt

(9b)

(14a)

(14b)

n pl

ϕp =

n sA + n sB + npl

B dX 0,1

· (9c)

Here, nsA, nsB and np are the numbers of moles of the solvents and polymer. On the other hand, the activities of other volatile components can be treated as constants for simplicity, and their influence on the total reactor pressure is often found to be minor. This simplified VLE model is valid because the dominating liquid components are monomers A and B and polymer for the majority of time. The reactor model developed above is a DAE system, consisting of population, monomer, and total mass balances, as well as volume and VLE calculations. The size of the model can be considerably large when high MW polymers are modeled, since the detailed chain length distribution is calculated. Reformulated Reactor Model. The foregoing DAE model of the copolymerization process is very stiff due to the presence of the very fast exchange reactions, compared with the rest of the reaction rates in the differential population balance equations. This poses great challenges in the numerical solution procedure, where simultaneous stable propagation of both fast and slow dynamic modes is required. To overcome this, we proposed a nullspace projection method9 to systematically reformulate the reaction equation system by isolating the fast equilibrium reactions from differential balance equations and modeling their quasi-steady states with algebraic equations. The reformulation procedure can be carried out with matrix manipulations, since the reaction rates (dynamic modes) are linearly connected by stoichiometry relations. As a result of the reformulation procedure, two pseudospecies X and Y are introduced: X nSA , nB = GnSA , nB + DnSA , nB ,

S ∈ {A, B}, nA , nB ≥ 0

Y nSA , nB

S ∈ {A, B}, nA , nB ≥ 0

=

Q nS , n A B

+

R nSA , nB ,

A dX1,1

dt

nc =

B dX1,1

dt

(GnA , nB + Q n

A

B

dt

ni =

dX nBA , nB dt

dY0,0 dt A dY1,0

dt B dY0,1

dt A dY1,1

(10)

dt

∑ ∑

= V −1[−(k iAM A + k iBMB)Q 0,0 + ktncM A ]

A , nB

(14h)

A ] = V −1[k iAQ 0,0M A − (k pAAM A + k pABMB)Q 1,0

B ] = V −1[k iBQ 0,0MB − (k pBAM A + k pBBMB)Q 0,1

B A M A − (k pAAM A + k pABMB)Q 1,1 ] = V −1[k pBAQ 0,1

B dY1,1

dt

A B B M − (k pBAM A + k pBBMB)Q 1,1 = V −1[k pABQ 1,0 ]

(14l)

(11)

dYnAA , nB dt −

= V −1[(k pAAQ nA − 1, n + k pBAQ nB − 1, n )M A A

(k pAAM A

+

B

A

k pABMB)Q nA , n ], A B

B

nA ≥ 2, nB ≥ 0 (14m)

(12)

(Q n

(14g)

(14k)

dYnBA , nB



nA = 0 nB = 0

nA ≥ 0, nB ≥ 2

(14j)

and the total moles of the byproduct chains nu =

(14f)

(14i)

(GnA , nB + DnA , nB)

nA = 0 nB = 0



nA ≥ 2, nB ≥ 0

= V −1[(k pABGnAA , nB − 1 + k pBBGnBA , nB − 1)MB

− (k pBAM A + k pBBMB)GnBA , nB],



∑ ∑

= V −1[(k pAAGnAA − 1, nB + k pBAGnBA − 1, nB)M A

− (k pAAM A + k pABMB)GnAA , nB],

and additionally, we introduce ni as the total number of moles of the initiator ∞

A B MB − (k pBAM A + k pBBMB)G1,1 = V −1[k pABG1,0 ]

dX nAA , nB

) ,n

nA = 0 nB = 0

B A M A − (k pAAM A + k pABMB)G1,1 ] = V −1[k pBAG0,1

(14e)



∑ ∑

(14c)

(14d)

We introduce the following notation: first, the total amount of catalyst is equal to the amount of chains attached with base metal ions, denoted as ∞

B = V −1[k iBG0,0MB − (k pBAM A + k pBBMB)G0,1 ]

dt

dt

+ R nA , nB)·

= V −1[(k pABQ nA , n A

B−1

+ k pBBQ nB , n A

− (k pBAM A + k pBBMB)Q nB , n ],

(13)

A

After reformulation, the differential population balance equations for X and Y are expressed in eqs 14, as shown below:

B

)M B

B−1

nA ≥ 0, nB ≥ 2

(14n)

Meanwhile, the algebraic equations giving the quasi-steady states of the exchange reactions are shown as: 7438

dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446

Industrial & Engineering Chemistry Research X nSA , nBnc = GnSA , nB(n i + nu),

S ∈ {A, B}, nA , nB ≥ 0

Y nSA , nBnc = Q nS , n (n i + nu),

S ∈ {A, B}, nA , nB ≥ 0

A

B

Article k−1 ⎛ ⎛k ⎞ dλkA = V −1⎜⎜k iAG0M A + k pAA ∑ ⎜ ⎟ζiAM A ⎝ ⎠ dt ⎝ i=0 i k ⎞ ⎛k ⎞ + k pBA ∑ ⎜ ⎟ζiBM A − k pABζkAMB⎟⎟ , ⎝ ⎠ ⎠ i=0 i

(15)

Note that the dormant species D and R are no longer explicitly included in the reformulated model. In sum, the reformulated model comprises three major building blocks: population balances eq 14 quasi-steady states eq 15 and definitions in eqs 11−13 additional equations monomer balance (eq 3) volume determination (eqs 4 and 5) VLE relations (eqs 6−9) The reformulated model has fewer differential equations and variables, and it is less stiff since the fast exchange reaction rates are eliminated from the differential population balances. Nevertheless, the entire chain length distribution is recorded. In practice, the distribution needs to be truncated by giving upper bounds to the number of repeating units in a polymer chain. In the context of copolymerization, appropriate upper bounds are required for both monomers. The reformulated copolymerization model is also of the same order of magnitude of the size as the original one.

(17d)

⎛ ⎛k ⎞ dλkB = V −1⎜⎜k iBG0MB + k pBB ∑ ⎜ ⎟ζiBM A ⎝ ⎠ dt ⎝ i=0 i k−1

k ⎞ ⎛k ⎞ + k pAB ∑ ⎜ ⎟ζiAMB − k pBAζkBM A ⎟⎟ , ⎝ ⎠ ⎠ i=0 i

dY0 = V −1( −k iAQ 0M A − k iBQ 0MB + k tncM A ) dt dμ0A dt dμ0B dt dμkA

MOMENT MODEL The method of moments is a very well-known method for solving polymerization systems with very large number of individual species. The moment equations are derived from aggregating population balances with different weights. The obtained moment model retains information for tracking average polymer properties. The method of moments can be applied well to linear polymers. The following notation is introduced to develop the moment model for the copolymerization process: ζ moment of growing product chains (G) v moment of growing byproduct chains (Q) λ moment of product chains (X) μ moment of byproduct chains (Y) The kth moment of a polymer species (e.g., G) is defined as

dt

dλ 0A = V −1(k iAG0M A + k pBAζ0BM A − k pABζ0AMB) dt

= V −1(k iBQ 0MB + k pABν0AMB − k pBAν0BM A )

(17h)

k−1 ⎛ ⎛k ⎞ = V −1⎜⎜k iAQ 0M A + k pAA ∑ ⎜ ⎟νiAM A ⎝ ⎠ ⎝ i=0 i

⎛ ⎛k ⎞ = V −1⎜⎜k iBQ 0MB + k pBB ∑ ⎜ ⎟νiBM A ⎝ ⎠ dt ⎝ i=0 i k−1

k ⎞ ⎛k ⎞ + k pAB ∑ ⎜ ⎟νiAMB − k pBAνkBM A ⎟⎟ , ⎝ ⎠ ⎠ i=0 i

k = 1, 2, ... (17j)

In addition, the algebraic equations for the quasi-steady states of the exchange reactions are obtained:

(16)

X 0nc = G0(n i + nu)

(18a)

Y0nc = Q 0(n i + n u)

(18b)

λkSnc = ζkS(n i + n u),

S = {A, B}, k = 0, 1, ...

(18c)

μkS nc = νkS(n i + n u),

S = {A, B}, k = 0, 1, ...

(18d)

and the definitions of the total amounts of the catalyst, initiator, and byproduct are rewritten as follows: nc = G0 + Q 0 +

∑ S ∈ {A,B}

(17a)

ni = X0 +



n u = Y0 +

(17b)



ζ0S + ν0S (19a)

λ 0S (19b)

S ∈ {A,B}

S ∈ {A,B}

dλ 0B = V −1(k iBG0MB + k pABζ0AMB − k pBAζ0BM A ) dt

k = 1, 2, ... (17i)

where, n represents the total number of repeating units irrespective of monomer types, viz. n = nA + nB, and Gn is the population of living product chains of length n. The rest of the moments can be defined analogously. In this study, the moment model for the copolymerization of A and B needs to account for the effect of the chain-ends on propagation. Therefore, superscripts A and B are also used for the moment notation to designate ending units. In eqs 17, the moment balances are derived: dX 0 = V −1( −k iAG0M A − k iBG0MB) dt

(17g)

dμkB

k = 0, 1, ...

n=1

(17f)

= V −1(k iAQ 0M A + k pBAν0BM A − k pABν0AMB)

k ⎞ ⎛k ⎞ + k pBA ∑ ⎜ ⎟νiBM A − k pABνkAMB⎟⎟ , ⎝ ⎠ ⎠ i=0 i



∑ nkGn ,

k = 1, 2, ... (17e)



ζk =

k = 1, 2, ...

μ0S

(19c)

Lastly, the monomer balance equations for A and B are shown as

(17c) 7439

dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446

Industrial & Engineering Chemistry Research

Article

dM A = F A − V −1[k iA(G0 + Q 0) + k pAA(ζ0A + ν0A ) dt +

k pBA(ζ0B

+

ν0B)

+ k tnc]M

A

λ1 λ0

(21a)

M w = MW

λ2 λ1

(21b)

(20a)

dM B = F B − V −1[k iB(G0 + Q 0) + k pAB(ζ0A + ν0A ) dt + k pBB(ζ0B + ν0B)]MB

M n = MW

The molecular weight of the repeating units is calculated as a weighted average of the molecular weights of A and B: (20b)

MW = w MWA + (1 − w)MWB

The constructed moment model for copolymerization is capable of predicting typically used polymer property indices. This assembles a reactor model that is significantly smaller than the full reactor model based on species balances. Although the detailed chain length distribution information is lost, the moment model is well suited for applications that only require average polymer properties. This becomes especially useful when computation resources are constrained and computation times are critical, such as (online) process optimization. Process Recipe Optimization. The polymer industry has extensive interests in developing optimal process operation policies with desired product properties. Optimizing a (semi)batch polymerization process, such as anionic ring-opening polymerization, generally requires the determination of an optimal control trajectory to drive the process to a desired final state. The obtained solution of the off-line optimal control problem renders the optimized process recipe. The essential building blocks to fulfill this task include (1) an accurate process model, (2) a selected number of control variables, (3) an objective function, and (4) a suitable numerical method for solving the specific optimization problem.15 The population balance and moment models are able to provide predictions of the reactor condition, if the associated model parameters, the kinetic constants, and thermodynamic parameters are known with sufficient accuracy. The parameter values used in this study are collected from different open literature sources. However, these published parameters were obtained under many different experimental conditions, and they do not necessarily reflect the true operating conditions of industrial polymerization reactors, because industrial reactors are operated over a wider temperature range, compared with laboratory environments. To address the inconsistency, we take advantage of historical plant data and improve the model fidelity by adjusting the corresponding model parameters to fit plant measurements under the same input conditions. This validation process can be done systematically by solving nonlinear parameter estimation problems.16 However, confidentiality considerations prevent us from presenting more detailed information here. Nevertheless, we show that the developed reactor models are able to adequately represent the polymerization process through adjustments of key kinetic and VLE parameters for a known set of plant data in the case study example. In this study, the control variables are the polymerization temperature and monomer feed rates. The objective function can be formulated to optimize the process economics using a wide array of indices, for example, the minimum batch time. Several process constraints are introduced to specify the demanded final product properties and the safety operation regulations. Commercial polymers have a set of widely used property indices. First, the number average molecular weight Mn and weight average molecular weight Mw can be obtained with the moments:

(22a)

t

w=

∫t F A(τ )dτ − M A (t ) 0

t

∫t [F A(τ ) + F B(τ )]dτ − M A (t ) − MB(t ) 0

(22b)

In eq 22b, w represents the mole fraction of A in the polymer chains. PDI is used as a measurement of the heterogeneity of the polymer to characterize the spread of the polymer chain length distribution, calculated by Mw/Mn. Narrow distribution, corresponding to low PDI, is preferred in many applications. On the other hand, the byproduct chains created by the transfer reaction are measured in milliequivalents per total mass: n α = 1000 u (23) m The percentage of polymer chains ending with monomer B in the external block is also an important property for the block copolymer, which can be calculated from the zeroth moments: ε (%) =

λ 0B + μ0B ∑S ∈ {A,B} (λ 0S + μ0S )

× 100 (24)

Also, the unreacted monomer in the product mixture at final time should stay below proper upper limits. Conventionally, it is measured in parts per million (ppm), defined as β S = MWS

MS × 106 , m

S ∈ {A, B}

(25)

A simplified heat balance is employed to derive the process safety constraints by considering the propagation reaction as the only source of energy: d(mHb) = dt



[F SΔHfS + rpS( −ΔHpS)]MWS − q

S ∈ {A,B}

(26)

where rp is the lumped rate of propagation reactions (superscripts represent different monomers added), and the heat of reaction −ΔHp for both monomer additions is assumed to be constant. Also, q denotes the heat removal rate from the heat exchanger that can be determined by using the overall heat transfer coefficient U and area A: q = UA(T − Tw )

(27)

where Tw is the temperature of the water used by the heat exchanger, and it is assumed to be constant. In eq 26, the enthalpies of the feed flow and bulk liquid are Hf and Hb, respectively, which are obtained as averages of corresponding heat capacities over the selected temperature range. For the first safety constraint, the heat removal duty cannot exceed the allowed maximum cooling capacity of the heat exchanger attached to the reactor: 7440

dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446

Industrial & Engineering Chemistry Research



Article

that creates a fully discretized version of the original dynamic system via OCFE, where the discretized problem can be optimized with of fthe-shelf nonlinear optimizers, such as IPOPT. In sum, the solution methodology of this study is similar to our previous work.9 Case Study. We test our modeling and optimization framework with an example study of two-thousand dalton A−B block copolymers. The internal A block is produced in tandem with the external B block without the degassing step, but there is a specified minimum percentage of B ended chains that needs to be met. A real-world process recipe obtained from plant data is used as the base case in the following study. We validate the population balance model (PBM) and moment model (MM) with the base case plant data, and recipe optimization is performed by using both models. For model implementation, both PBM and MM are implemented in the General Algebraic Modeling System (GAMS)19 platform after discretization into NLP form. The continuous time domain is separated into two stages, where A and B are consecutively used as feed monomers. The lengths of the two stages are decisions to make for the optimizer. Furthermore, the discretization settings for the two models are different: a finer mesh of 40 finite elements along with a twopoint Radau collocation is applied to MM for each stage, while a coarser one of 10 finite elements with three-point Radau collocation roots is introduced for PBM due to its large model size. To record the entire chain length distribution of the PBM, the upper bounds for repeating units are specified as NA = 30 and NB = 14. It is worth noting that a lower order collocation is used in MM, in conjunction with a larger number of finite elements; this in general excels in handling path constraints, particularly for high index ones, in the optimization problem.

rpS( −ΔHpS)MWS

S ∈ {A,B}





F S( −ΔHfS)MWS + UA(T − Tw ) (28)

S ∈ {A,B}

Here, we assume the monomer feed enters at a constant ambient temperature Tm = 25 °C, which is lower than the reactor temperature, offering extra cooling capability. Second, for the exothermic polymerization process, the amount of unreacted monomers should be carefully controlled to minimize the risk of product decomposition, when cooling capability is lost during polymerization. To carry out such a task, it is conducive to add a constraint on the adiabatic end temperature, which equals the summation of the current reactor temperature and the potential adiabatic temperature rise due to the latent heat existing in the unreacted monomers. When heat exchange stops at time tc with the reactor temperature noted as Tc, it is clear that q = 0 and also reasonable to assume F = 0 (stop feeding monomers) after tc. Therefore, integrating eq 26 starting from tc to the steady state (infinity) gives m[Hb(Tad) − Hb(Tc)] =



M S( −ΔHpS)MWS

S ∈ {A,B}

(29)

where Tad is the adiabatic end temperature and we assume all the monomers are consumed in the propagation reactions. For the safety limit of the adiabatic end temperature, a sufficient safety margin is also needed to tolerate uncertainties. Meanwhile, monomer B is potentially explosive in nitrogen, and a linear limit extrapolation of the explosive region boundary is given by yB = 1.3865 − 0.001764T − 0.0003568P

(30)



and yB represents the vapor phase concentration of B. The recipe optimization problem can be cast into a generic dynamic optimization formulation given by min Φ(z(tf ), tf ) u ∈

MODEL VALIDATION In the polymerization process, the reactor pressure is a measurement that is relatively easy to access. We use the total pressure in the reactor as the criterion to illustrate the model validation procedure in our study. Besides the VLE equations described in eqs 6−9, two additional issues need to be addressed in order to obtain pressure predictions correctly. With a known initial amount of nitrogen nN2, the partial pressure over nitrogen can be calculated by using the ideal gas law. Second, when reactor venting occurs at critical reactor pressures, we assume a negligible loss of the monomer molecules in the vapor phase. Moreover, the escaped nitrogen amount can be deduced as the gas phase volume in the reactor decreases, and while the reactor pressure stays constant at the maximum allowed by the vent system. The batch time of the base case recipe is normalized to 1, including an A feeding period followed by a B feeding period. The reactor temperature and monomer feed rates are recorded periodically. To compare with, we consider the model predictions from the PBM and MM under the same operating conditions, and the result is shown in Figure 2. The match between the prediction curves and data trajectory is considered satisfactory. Parameter adjustments were made with respect to reaction kinetics. All the kinetic constants are first obtained from published articles. Next the pre-exponential factors are tuned to best fit the model predicted pressure to the plant data. The same kinetics and VLE parameters are used in the PBM and MM. In fact, additional model validation work has been carried out to compare polyol property indexes. The result is also satisfactory but is omitted here for confidentiality reasons.

(31a)

z ̇ = f (z(t ), y(t ), u(t )), z(0) = z 0 s.t.

g (z(t ), y(t ), u(t )) = 0 h(z(t ), y(t ), u(t )) ≥ 0 t ∈ [0, tf ]

(31b)

Here, z and y are differential and algebraic state variables respectively, and u denotes the control variables that in this study include the polymerization temperature and monomer feed rates. The reactor model is represented as a DAE system with differential equations f(·) and algebraic equations g(·). Note that the reactor model can be either the large-scale population balance model or succinct moment model. The process constraints are noted by h (·). The constraint h (·) can be enforced either at final time tf (usually corresponding to the end product properties) or along the evolution path of the dynamic system t ∈ [0,tf] (referring to the process safety constraints). In this study, the initial condition of the differential variables z0 is given by the charge condition of the batch. The total length of the operation equals tf. The objective function is noted in the Mayer form, and minimizing the final batch time tf is a typical choice. A comprehensive investigation of numerical solution methods for dynamic optimization problems can be found in Biegler.18 In particular, we are interested in the simultaneous collocation method 7441

dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446

Industrial & Engineering Chemistry Research

Article

Figure 2. Reactor pressure profiles for model validation.

Figure 3. Optimal control profiles of the process from the moment model.



RECIPE OPTIMIZATION RESULTS

Table 2. Recipe Optimization Statistics and Results from the Moment Model

Recipe optimization is carried out with respect to the same polymerization system, with the following constraints added: (1) product property constraints including the minimum numberaverage molecular weight, maximum byproduct ratio, maximum unreacted monomers, and minimum level of B ended chains. The specified threshold values are determined from the base case recipe; (2) the maximum heat removal duty of the heat exchange is expressed as UA(T − Tw) kW; (3) the upper limit of the adiabatic end temperature is (Tb + 100 °C); (4) the allowed maximum reactor temperature is (Tb + 80 °C). Here, the threshold values on product quality are obtained from simulating the base case recipe, in which the polymer model is solved with the batch time and controls fixed to their recipe values. Among them, UA and Tb are specified constant parameters. The safety constraints are feasible in the simulation. We first show the result for optimizing the moment model, followed by the result from the population balance model. The optimization problems are solved with GAMS/IPOPT to proved local optimality with appropriate initialization of the participating variables; IPOPT runs with linear solver MA86. All computations

optimized soln

no. of variables

no. of constraints

CPU(s)

0.580

11,248

11,121

22

are performed on a desktop with an 8 core 2.80 GHz Inteli7 processor, 9 GB memory, and installed Linux kernel 3.2.0−34. Optimization with Moment Models. The moment model is well suited for recipe optimization owing to its small scale and ability to predict product properties. The optimized recipe is able to reduce the copolymerization processing time by 42.0%. Details on the statistics and solution of the model can be found in Table 2. The computation load is quite manageable for optimization with MM. The operating strategies for the copolymerization process are depicted in Figure 3. For the base case recipe, the polymerization temperature is designed to remain constant over time for the A and B addition stages. The recipe design pattern for the monomer feed policy is as follows: A and B are fed consecutively with a noticeable gap in between. Each feeding window begins with a ramping period, where the feed rate increases roughly linearly, and then the rate is kept at a desired constant level. The polymerization is continued 7442

dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446

Industrial & Engineering Chemistry Research

Article

Figure 4. Process constraint profiles from the moment model.

Figure 5. Polymer properties from the moment model: molecular weights and polydispersity indices.

concentration in the vapor phase. The value of the adiabatic end temperature indicates the potential effect of the latent heat existing in unreacted monomers on the reactor temperature. This constraint is not active since the monomer inventories are well controlled. On the contrary, the reactor cooling capacity is the primary limiting factor for further improvement of the process performance and the corresponding constraint is active during most of the operation time horizon, especially for the A feeding stage. The total cooling capacity also includes a portion provided

when there is no monomer feed entering for a certain time interval; these are called digestion periods. The optimized recipe redesigns the temperature and feed rate profiles. The reactor temperature rises sharply during the A digestion and B feeding period, reaching its upper limit. In the feed profile, B is gradually added after A, and the digestion period of B is almost negligible by virtue of its high reactivity. Figure 4 shows the three important constraining factors of the process: the adiabatic end temperature, heat removal rate, and B 7443

dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446

Industrial & Engineering Chemistry Research

Article

Figure 6. Polymer properties from the moment model: important quality indices.

Figure 7. Chain length distributions of the optimized recipe of the PBM.

by the monomer feed streams that are of a lower temperature than the reactor. Finally, we inspect the safety constraint regarding the vapor phase composition distribution, and the system is found to be within the safety zone predicted by eq 30. To demonstrate the product properties, we first show the number average molecular weight and PDI as functions of polymerization time in Figure 5. The product polymers produced by the optimized recipe give the same number average molecular weight and a PDI very close to the base case products, which means the MWDs of the two cases are rather similar. For the byproduct, the differences between the two recipes are also minor. In addition, Figure 6 shows the evolution of the two quality indices: impurity level and B-ended chain ratio. The product from the optimized recipe shares the same quality level as the base case product. It is worth noting that the constraint on the impurity level is active for the optimized recipe, and relaxing the limit value may further reduce the batch time. Optimization with Population Balance Models. Recipe optimization over the PBM is a challenging task that requires significant computation effort. To facilitate the optimization, we adopt the optimized recipe from the MM and start the optimizer at that solution. The initial guesses for population states are

Table 3. Recipe Optimization Statistics and Results from the Population Balance Model optimized soln (h)

no. of variables

no. of constraints

CPU(s)

0.588

191,408

185,147

5221

obtained through dynamic simulation beforehand, and it can be also verified that the optimal recipe generated by the MM remains a good feasible solution for the PBM. By this means, a large-scale optimization problem with the PBM is able to be solved by IPOPT, and the detailed information is listed in Table 3. The optimal batch time is slightly worse than the one obtained with the MM, and the end-point values of the final time process constraints are similar to the previous case. However, there is a huge expansion of the model scale in terms of variable and constraint numbers, despite a less dense discretization mesh being used. Consequently, the optimization problem becomes slow to solve. It should be noted that the base case recipe and the set of process constraints are chosen to illustrate the use of dynamic optimization and do not necessarily reflect the true capability or limitations of the plant. However, the actual potential reaction time saving can still be significant. 7444

dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446

Industrial & Engineering Chemistry Research

Article

n N P q Q r rA rB R t tc T Tad Tb Tc Tw U V w X yB Y ζ λ μ ν ϕ χ

The optimal control recipe generated by using the PBM is very similar to the optimal MM solution. In addition, the PBM model is able to reveal the chain length distribution with three-dimensional plots shown in Figure 7 for the product and byproduct chains. The complete distribution information is helpful in predicting physical properties such as viscosity of the polymer product.



CONCLUSIONS In this study, we have developed a modeling and optimization methodology for semibatch ring-opening polymerization processes for block copolymers. We followed two different approaches to model the polymer species: the population balance approach and method of moments. In the case study example, we first demonstrated that both reactor models were able to match real plant data adequately. The following recipe optimization results showed the potential of our modeling and optimization framework to enhance the process performance by redesigning the reactor operating policy. In particular, the moment model excelled with respect to computational performance. The population balance model can still be appreciated for its ability to uncover the chain length distribution. We adopted the same numerical solution strategy applied in our previous work9 in solving the recipe optimal design problem, featured by the simultaneous collocation method that converted DAEs to nonlinear algebraic equations and interior point method that solved the resulting nonlinear programs.



AUTHOR INFORMATION

number of moles, mol maximum number of repeating units pressure, kPa heat removal rate, J/s growing byproduct chains, mol reaction rate, mol/S AB reactivity ratio kAA p /kp , dimensionless BB BA reactivity ratio kp /kp , dimensionless dormant byproduct chains, mol time, s time (loss of cooling), s reactor temperature, K adiabatic end temperature, K base temperature, K reactor temperature (loss of cooling), K water temperature, K overall heat transfer coefficient, W/m2·K liquid volume, m3 mole fraction of A in polymers, dimensionless total product chains, mol vapor phase B concentration, dimensionless total byproduct chains, mol growing product polymer moment, mol product polymer moment, mol byproduct polymer moment, mol growing byproduct polymer moment, mol lattice fraction, dimensionless interaction parameter, dimensionless

Corresponding Author

Subscripts Substances

Notes

b c f i m n p s S u

*E-mail: [email protected]. The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS We would like to thank Dr. Paul Witt (The Dow Chemical Company) for his valuable comments on this work. DEDICATION We would like to dedicate this work to the memory of John Congalidis, whose contributions to Polymer Reaction Engineering at the intersection of academic and industrial research have inspired and shaped our work and that of many others in the field.

bulk catalyst feed initiator repeating units repeating units polymer solvent monomer species byproduct chains

Reactions

i p e t



NOMENCLATURE a liquid phase activity, dimensionless A heat transfer area, m2 Ar pre-exponential factor, m3/mol s Cp specific heat capacity, J/g K D dormant product chains, mol Er activation energy, J/mol f functionality, dimensionless F monomer feed rate, mol/s G growing product chains, mol H enthalpy, J/g ΔH reaction heat, J/g k reaction rate constant, m3/mol s l number average chain length, dimensionless m total mass, g M monomer, mol Mn number-average molecular weight of polymer, g/mol Mw weight-average molecular weight of polymer, g/mol MW molecular weight, g/mol

initiation propagation exchange transfer

Superscripts

A monomer A B monomer B S monomer species



REFERENCES

(1) Odian, G. Principles of Polymerization; John Wiley and Sons: New York, NY, 1991. (2) Hungenberg, K. In Handbook of Polymer Reaction Engineering; Wiley Online Library: New York, NY, 2008. (3) Abel, O.; Helbig, A.; Marquardt, W.; Zwick, H.; Daszkowski, T. Productivity optimization of an industrial semi-batch polymerization reactor under safety constraints. J. Process Control 2000, 10, 351−362. (4) Ray, W. H. On the Mathematical Modeling of Polymerization Reactors. J. Macromol. Sci., Part C: Polym. Rev. 1972, 8, 1−56. (5) Zavala, V.; Biegler, L. Optimization-based strategies for the operation of low-density polyethylene tubular reactors: nonlinear model predictive control. Comput. Chem. Eng. 2009, 33, 1735−1746.

7445

dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446

Industrial & Engineering Chemistry Research

Article

(6) Chatzidoukas, C.; Kiparissides, C.; Perkins, J.; Pistikopoulos, E. Optimal grade transition campaign scheduling in a gas-phase polyolefin FBR using mixed integer dynamic optimization. Comput.-Aided Chem. Eng. 2003, 14, 71−76. (7) Cervantes, A.; Tonelli, S.; Brandolin, A.; Bandoni, J.; Biegler, L. Large-scale dynamic optimization for grade transitions in a low density polyethylene plant. Comput. Chem. Eng. 2002, 26, 227−237. (8) Flores-Tlacuahuac, A.; Biegler, L.; Saldívar-Guerra, E. Dynamic optimization of HIPS open-loop unstable polymerization reactors. Ind. Eng. Chem. Res. 2005, 44, 2659−2674. (9) Nie, Y.; Biegler, L.; Villa, C.; Wassick, J. Reactor modeling and recipe optimization of polyether polyol processes: polypropylene glycol. AIChE J. 2013, 59, 2515−2529. (10) Biegler, L. An overview of simultaneous strategies for dynamic optimization. Chem. Eng. Process.: Process Intensif. 2007, 46, 1043−1053. (11) Wächter, A.; Biegler, L. On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming. Mathem. Program. 2006, 106, 25−57. (12) Villa, C. Reactor modeling for polymerization processes. Ind. Eng. Chem. Res. 2007, 46, 5815−5823. (13) Flory, P. Principles of Polymer Chemistry; Cornell Univ. Press: New York, NY, 1953. (14) Favre, E.; Nguyen, Q.; Clement, R.; Neel, J. Application of Flory− Huggins theory to ternary polymer-solvents equilibria: A case study. Eur. Polym. J. 1996, 32, 303−309. (15) Kiparissides, C. Polymerization reactor modeling: A review of recent developments and future directions. Chem. Eng. Sci. 1996, 51, 1637−1659. (16) Bindlish, R.; Rawlings, J.; Young, R. Parameter estimation for industrial polymerization processes. AIChE J. 2004, 49, 2071−2078. (17) Mayo, F.; Lewis, F. Copolymerization. I. A basis for comparing the behavior of monomers in copolymerization; the copolymerization of styrene and methyl methacrylate. J. Am. Chem. Soc. 1944, 66, 1594− 1601. (18) Biegler, L. Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes; Society for Industrial and Applied Mathematics: Philadelphia, PA, US, 2010. (19) Brooke, A.; Kendrick, D.; Meeraus, A.; Raman, R.; Rosenthal, R. GAMS a User’s Guide. GAMS Development Corporation: Washington DC, US, 2006.

7446

dx.doi.org/10.1021/ie402770k | Ind. Eng. Chem. Res. 2014, 53, 7434−7446