Modeling the Distribution of Functional Groups in Semibatch Radical

Subscriber access provided by TUFTS UNIV

Kinetics, Catalysis, and Reaction Engineering

Modeling the Distribution of Functional Groups in Semi-Batch Radical Copolymerization: An Accelerated Stochastic Approach Amin Nasresfahani, and Robin A. Hutchinson Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b01943 • Publication Date (Web): 26 Jun 2018 Downloaded from http://pubs.acs.org on June 27, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Modeling the Distribution of Functional Groups in Semi-Batch Radical Copolymerization: An Accelerated Stochastic Approach Amin Nasresfahani, Robin A. Hutchinson* ––––––––– Department of Chemical Engineering, Queen’s University, Kingston, ON K7L 3N6, Canada E-mail: [email protected] ––––––––– Abstract While Kinetic Monte Carlo (KMC) techniques provide a powerful means to model polymer microstructure, the associated computational cost has been a barrier to its wide-spread adoption. The case of radical semi-batch polymerization under starved-feed policy is a particularly challenging application: at the initial stage, a large simulation volume is required to accurately represent the low concentration of radicals generated at the start of the reaction, while the reactant feed dictates the further increase of the simulation volume with time. A combination of approaches is implemented in a stepwise fashion to greatly accelerate the KMC representation of this system. First, a correction factor is developed to maintain a constant simulation volume in order to improve the efficiency of the solution, followed by scaling of the reaction rates to preserve accuracy at low control volumes and further reduce computational effort. A novel strategy for storing the explicit chain sequences and parallel analysis of the stochastic data is also implemented, with the computational time required to accurately represent a semi-batch radical copolymerization test case reduced from 50 to less than 2 minutes. The accelerated stochastic approach provides a foundation for future optimization of feeding strategies to minimize the fraction of non-functionalized chains formed during the production of low molar mass copolymers.

-1ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1. Introduction Free radical polymerization (FRP) remains the predominant commercial route for the production of functionalized dispersants and surfactants due to its capability to produce a broad spectrum of materials from relatively inexpensive monomers in aprotic or protic solvents in the presence of impurities such as monomer stabilizers and oxygen.1 Solvent-borne automotive coatings and surfactant-free compositions used in pharmaceutical industries are examples of higher value materials industrially manufactured through FRP using operating strategies developed to ensure robust control of average copolymer composition and polymer molecular weights.2–6 Uniform operating conditions must be maintained, as the final product consists of a mixture of chains with short lifetimes (< 1s) produced at different times throughout the overall reaction, which often proceeds for several hours. To improve product uniformity, the synthesis strategy often adopted is a semi-batch starvedfeed reactor policy. The system is referred to as “starved” when the concentrations of reactants in the reaction media are kept low by setting the feed rates of the reactive species to roughly match their consumption rates.7–11 Production of low molecular-weight coatings is a common example of using this approach at industrial scale. In contrast to a batch process, variations in copolymer composition and molecular weight are substantially avoided by maintaining a nearly constant concentration ratio of the monomers under semi-batch operation.12 Although average copolymer composition is controllable, the distribution of comonomer units among the chains is based on their random incorporation, restricting the opportunity of tailoring microstructure. Chains of higher molecular weight (MW) contain a sufficient number of repeat units (e.g., >100) such that their individual compositions closely match the average composition of all chains produced at that instant in time. However, the random placement of functional groups creates more of a problem when producing low MW material by FRP, as the probability of generating chains without the desired functionality increases with decreasing chain-length in the distribution.13–17 While average composition can be modeled using deterministic methods, stochastic modeling tools such as the Kinetic Monte Carlo (KMC) methodology are required to track the


Page 2 of 29

Page 3 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


placement of comonomer (often functional groups) across the polymer chain-length distribution.16,18– 21

Utilizing a probability function based on rate equations to predict reaction events, KMC is capable

of predicting the explicit sequences of the individual chains. Such detailed information allows the estimation of the full polymer weight distribution under conditions of stationary and non-stationary radical generation, as well as the specific locations of cross-links, pendant double bonds, and comonomer units (among other features) in the population of chains that can be compared to measured properties such as swelling behavior and gel fractions, and bivariate copolymer composition–chain length distributions inferred by crystallization behavior, or electron spray ionization mass spectra.18,22–24 However, the time-consuming computational calculations required to obtain this detailed information often limit the utilization of the KMC method. To maintain accuracy, the simulation volume must be large enough such that the number of molecules of each species (e.g., monomers, radicals, polymer chains, etc.) remains representative of its actual population, and the number of species required correlates directly with the computational effort. As their concentration is orders of magnitude lower than the other components, it is the treatment of radicals that becomes critical for an accurate representation of FRP. In fact, the average number of radicals in the system must be kept at two or greater in order to properly address bimolecular termination reactions and predict precisely the reaction rates.16,17,25,26 The objective of this work is to reduce the computational workload required to represent semi-batch FRP by introducing a correction factor to maintain a constant simulation volume despite the addition of fresh feed. This approach is outlined in the following sections, and combined with other strategies suggested in recent literature16,21,26–29 to reduce the computing time while preserving the precision of the KMC calculations. Further discussion of these techniques is provided below as we describe the stepwise improvements implemented to accelerate the solution of a case study previously developed to represent the distribution of functional groups in a low MW copolymer produced by semi-batch operation.16 The final implementation is capable of reducing the computational time from 50 to less than 2 minutes while maintaining a high level of accuracy.



Page 4 of 29

2. Description of the Copolymerization Procedure As described previously,16 the semi-batch test-case recipe was developed to produce a copolymer consisting predominantly of n-butyl methacrylate (BMA) with just enough glycidyl methacrylate (GMA) comonomer to yield chains with an average of one GMA unit per chain. With a number-average chain-length of 20, thus, the monomer mole fraction is 5%. Adapted from a patent filed by Barsotti et al.,30 a solution of tert-butyl peracetate (TBPA), BMA, and GMA is continuously fed to a 1.0 L semi-batch reactor for a period of 4 hours. The reactor, well-controlled at 138 °C, is filled with 230 mL (200 g) xylene prior to feeding a total mass of 500 g consisting of 0.9082 g/g BMA, 0.0478 g/g GMA, and 0.044 g/g TBPA. The addition of feed, as well as the formation of polymers (of higher density than monomer), alters the reaction volume continuously during polymerization. The density of the copolymer formed is considered to be 1.078 g/mL.16

Table 1. The reaction mechanism and kinetic rate coefficients implemented to describe the radical copolymerization of BMA and GMA in a semi-batch reactor.16,31 Reaction Type

Reaction Scheme

Initiation

2

%

' = 4.69 × 10

$

67 = 0.266

# + $

#&

Chain Transfer

# + $

+ 2#

%

45

# + 3

+ 2# Termination

$

$

Propagation

01

= 1.32 × 10 ! = 0.515

+

Kinetic Parameters

$ ;2

0.00 0

20

40

60

80

100

120

Chain-Length

Figure 3. Time evolution of average molecular weights (A right axis), cumulative mole fractions of chains containing n GMA group denoted as F=n (A left axis), and monomer concentration profiles (inset B) computed for the ensembles attributed to S= 40 and 50. The number distributions (B) of the final polymer product are presented for all chains (Overall) as well as for chains with specific GMA functionality.

- 10 ACS Paragon Plus Environment

Page 11 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2. The total computational time requirements and the total cumulative average errors computed for KMC simulation of semi-batch BMA/GMA copolymerization as a function of sampling number. Errors are calculated using results for S=50 as reference “true” values.

3

AKLMN × 10?

Cumulative Average Errors* _à &_`R F

50 40 30 20 16 12 10 9 8 7 6 5 4 3 2 1

3.8569 3.0856 2.3142 1.5428 1.2342 0.9256 0.7714 0.6942 0.6171 0.5399 0.4628 0.3857 0.3085 0.2314 0.1542 0.0771

m kLnopN 6m k 6

∗ jk = l

m kLnopN 6

0.34 1.01 2.24 3.43 5.34 6.70 7.70 9.13 10.7 13.1 16.7 21.8 30.9 45.3 68.3

a

_`b &_`c

b

F

_de &_da &_dR &_dfg

0.05 0.21 0.66 0.85 1.53 1.89 2.23 2.56 2.93 3.65 4.66 6.21 7.55 8.42 26.0

:

0.18 0.32 0.57 0.97 1.38 1.79 2.06 2.35 2.85 3.89 5.26 7.44 10.1 14.9 27.2

c

h,Z ∗∗ ;,hH= 50 40 30 20 16 12 9.0 8.2 7.2 6.3 5.2 4.1 3.3 2.5 1.6 0.8

l × 100 for monomer concentrations,a average polymer MWs,b and chain

composition distributionsc ** Measured internally through MATLAB® stopwatch timer

Table 2 summarizes the reduction in computational time required as the control volume was reduced by lowering S from 50 to 1, as well as the associated increase in errors, calculated as a cumulative average of the least absolute errors with respect to the monomer concentrations, molecular-weight averages, and mole fractions of chains containing 0, 1, 2, or above 2 comonomer units (i.e., K] , K , KF , Kr ). Variations from the “true” values (as represented by the simulation results obtained from S=50) are significant for S≤10. Thus, the advantage of the reduced simulation time is accompanied by an increased error, as illustrated by Figure 4B. However, a significant amount of the error stems from the initial stages of the semi-batch reaction, as shown in Figure 4A. To preserve the desired accuracy, therefore, the sampling number should be large enough to capture the initial dynamic behavior of the reactor when radical concentrations are at their lowest.



0.20

0.014 GMA GMA GMA GMA

Concentration (mol/L)

0.18

0.16

BMA [Sampling Number = 5] BMA [Sampling Number = 10] BMA [Sampling Number = 20] BMA [Sampling Number = 50]

0.012

0.010 0.14

0.008

0.12

0.10

0.006

0.08 0.004 0.06 0

25

50

75

100

125

150

175

200

225

Time (min)

70

50 Time (min) Total Cumulative Ave. Error

60

50 30

40

30 20 20

Total Cumulative Ave. Error

40

Time (min)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 29

10 10

0

0 1

5

9

13

17

21

25

29

33

37

41

45

49

Sampling Number (S)

Figure 4. The monomer concentration profiles calculated for S = 5, 10, 20, and 50 (A). The mean time-averaged errors and computational times for KMC simulation of semi-batch copolymerization (B). See Table 2 for further details.

In Figure 4, the slope attributed to the total cumulative average errors begins to exponentially increase for S values below 20 while the curve associated with the computational time


Page 13 of 29

experiences almost no alteration in its slope. The increased error is due to the erroneous calculation of reaction rates executed by the stochastic technique at the lower simulation volumes. These inaccuracies can be rationalized by observing the patterns of radical birth during the simulation. Once the size of simulation volume is lowered below S=20, a zero-one radical condition becomes more frequent in the simulated ensemble, as indicated in Figure 5. Thus, the concentration of the radicals in the small ensemble becomes unrealistically high, embedding a significant error in the calculated rates. The simulation allows the radical to grow until the next radical arrives, at which point the termination reaction becomes dominant. Using a larger S value prevents the simulation from becoming trapped in situations where only one radical exists in the system (Figure 5). A scaling method to address this issue will be introduced in Section 3.4, after first presenting a means to address the increased calculational time associated with increasing simulation volumes. 20 [S=50]

[S=20]

[S=5]

[S=10]

[S=1]

18

Instantaneous Number of Radicals

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


16 14 12 10 8 6 4 2 0 0

50

100

150

200

250

Time (min)

Figure 5. The instantaneous number of radicals contained within different control volumes as a function of reaction time for the simulations selected from Table 2.

3.3 Time-dependent Control Volume with Adjusted Simulation Volume The volume of polymerization mixture increases continuously during semi-batch operation. Correspondingly, the simulation volume is being enlarged during the KMC calculation, having a



Page 14 of 29

direct influence on the computational workload. However, as shown by the set of simulations summarized in Table 2 and Figure 4A, it is the initial stage of the process, corresponding to the semi-batch start-up with no initiator contained in the initial reactor charge, that is the most difficult to accurately represent. With initiator molecules and thus radicals at their lowest values, the initial simulation volume must be chosen to be large enough to converge the stochastic iterations and preserve the accuracy of the results. However, the remainder of the simulation does not require such a large ensemble size. Once the initiator molecules have had a chance to accumulate in the system, the simulation volume becomes inevitably larger than the optimal size required, lengthening the calculation time. To balance the need for a larger initial control volume for accuracy with the reduction in computational effort, the ensemble has been discretized. A correction factor ;s ≤ 1= is introduced to adjust the simulation volume to its initial value throughout the calculation at set intervals. Thus, s is defined as the ratio between a prior simulation volume and the current one, with Δv chosen to be 60 seconds in order to maintain an almost-constant simulation volume: s=X

X`w ;6=

`w ;6&x6=

(2)

As the simulation volume is adjusted back to its initial value at every minute, the number of species is reduced by multiplying with s. In addition, the resizing factor AKLMN is no longer constant but decreases over the course of the simulation. Thus, the treatment of the physical system (constant volume reactor with increasing reaction volume) is reversed in this KMC implementation, which uses an (almost) constant simulation volume and a control volume that decreases from its high initial value (needed at the start of reaction due to the low number of radicals) gradually and continuously during the simulation. The two approaches are contrasted in Figure 6. With a constant control volume (as considered in Section 3.2), the simulation volume increases as new feed molecules are added to the system (Figure 6A), a straightforward implementation that mimics the physical process. The modified methodology keeps the simulation volume constant (reset to its initial value every minute), such that the control volume (representing the 1 L reactor) decreases continuously, as


Page 15 of 29

shown in Figure 6B. As also shown in the figure, both simulation and control volumes are 10 times smaller for S=5 compared to with S=50.

×10-17

×10-18 Control Volume

Simulation/Control Volume (L)

4.0

4.0

3.5

3.5

Sampling Number = 5

3.0

3.0

Sampling Number = 50 2.5

2.5

2.0

2.0

Simulation Volume

1.5

1.5

1.0

1.0 0

25

50

75

100

125

150

175

200

225

Time (min) ×10-17

×10-18

4.0

Sampling Number = 50 Sampling Number = 5

3.5

Simulation/Control Volume (L)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.0

4.0 3.5 3.0

Control Volume

2.5

2.5

2.0

2.0

1.5

1.5

1.0

1.0 Simulation Volume 0.5

0.5

0.0

0.0 0

25

50

75

100

125

150

175

200

225

Time (min)

Figure 6. Simulation volumes used to represent semi-batch reactor operation using a constant control volume with time-varying simulation volume (A), and a constant simulation volume



with time-dependent control volume (B). Results are shown for S=50 (left axis) and S=5 (right axis). Figure 7 illustrates the flowchart implemented, with the effect of these changes on the AKLMN values summarized in Table 3. Note that AKLMN changes by a factor of 2.5 for all values of S, with a corresponding change in 8TLMN . As seen by comparing the results to those in Table 2, the correction factor (s) introduced to keep simulation volume constant accelerates the solution considerably, with the time required to complete the simulation for the sampling numbers of 50, 40, 30, and 20 reduced by 13, 12, 9, and 6 minutes, respectively.

Figure 7. The flowchart implemented to maintain a constant KMC simulation volume to represent semi-batch operation.


Page 16 of 29

Page 17 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


As expected, the reduced computational time is accompanied by an increase in error. However, the error for S ≥ 20 increases only slightly compared to the Table 2 values, with solutions remaining accurate within 5%. The relative errors increase with decreasing S, as also observed for the “constant control volume, increasing simulation volume” results reported in Table 2. Thus, the correction factor provides a simple and effective approach in lowering the calculational cost, but only provided that the initial control volume is sufficiently large to suppress the errors associated with low radical numbers at start-up. Extending the approach to the lower values of S can be accomplished using recently-proposed acceleration methods discussed below.

Table 3. The total computational time requirements and the total error percentages computed for a series of sampling numbers via maintaining the reaction volume at a constant level. Errors are calculated using results for S=50 from Table 2 as reference “true” values.

3 50 40 30 20 16 12 10 9 8 7 6 5 4 3 2 1

Cumulative Average Errors*

AKLMN × 10?

_à &_`R

3.8569→1.5086 3.0856→1.2071 2.3142→0.9055 1.5428→0.6040 1.2342→0.4834 0.9256→0.3628 0.7714→0.3025 0.6942→0.2724 0.6171→0.2422 0.5399→0.2121 0.4628→0.1820 0.3857→0.1519 0.3085→0.1218 0.2314→0.0917 0.1542→0.0615 0.0771→0.0314 m kLnopN 6m k 6

∗ jk = l

m kLnopN 6

F

0.94 1.59 2.68 4.94 6.59 9.63 12.3 13.9 16.3 19.3 23.9 29.4 37.8 48.1 62.2 79.5

a

_`b &_`c

b

F

0.04 0.17 0.50 0.90 1.44 2.14 2.91 3.33 3.87 4.81 5.78 6.72 8.32 9.10 13.3 33.2

_de &_da &_dR &_dfg :

0.21 0.29 0.31 0.82 1.28 2.07 2.84 3.36 3.69 4.59 6.11 7.46 10.5 13.4 19.1 32.0

c

h,Z ∗∗ ;,hH= 37 28 21 14 8.7 6.4 6.3 4.6 4.0 3.5 3.0 3.0 2.0 1.5 1.2 0.7





Page 18 of 29

3.4 Improved Accuracy through Scaling of Reaction Rates With low control/simulation volumes, the KMC solution fails to accurately represent the macroscopic system due to consecutive inaccurate and divergent KMC decisions originating from an erroneous probability axis (related to the “zero-one” radical problem). The methodologies introduced to address this issue can be categorized into two approaches. The first, adopted in this work, attempts to reconcile the reaction rates in parallel with the KMC simulation,26,35 while the second category, known as hybrid methods, utilizes predetermined reaction rates calculated through a deterministic modeling strategy such as the method of moments.36–38 Both approaches allow a successful KMC solution at low simulation volumes by constructing a reliable probability axis. Herein, we implement the scaling method introduced to reconcile the reaction rates (in particular, the termination rates) and greatly improve accuracy for S ≤ 20. Combined with the time-dependent control volume (s correction factor from Equation 2), significant savings in computational time are achieved. The implementation is adapted from the scaling method developed by Gao et al.26 to greatly reduce the number of molecules in a KMC simulation volume while maintaining an accurate representation of reaction rate, accelerating simulation of free radical batch copolymerization by two orders of magnitude. The procedure starts with the idea that the simulation volume requires a minimum of two radicals throughout the course of polymerization to be representative of the termination phenomena; while decreasing the simulation volume further does not change the number of radicals produced from an initiation step, the sudden corresponding increase in the radical concentration within the simulation volume affects the ratio of termination to propagation rates such that the radicals have a higher propensity for termination compared to the “real” macroscopic system. As a result, the average molecular weights become lower than the true values. As developed by Gao et al.,26 z is a dimensionless, time-dependent scaling factor ( ≥ 1) defined by: F

z = V

0U0

k% 0 ∑$ | LQR N %

P

(3)

where L F N, C$ , }$ , and B6J6 are the initiator concentration, the mole fractions and molar densities of the monomers and solvent, and the total number of unreacted molecules (i.e., xylene, BMA, GMA, and initiator) within the simulation ensemble, respectively. 6 and are the continuum rate - 18 ACS Paragon Plus Environment

Page 19 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


coefficients. This scaling neutralizes the artificial influence of the erroneous radical concentration at simulation volumes below the minimum (i.e., the volume corresponding to two radicals) via dividing the reaction rates by z or z F , depending on whether the reaction is first order (e.g., propagation rates), or second order (e.g., termination rates) with respect to the radical concentration. Thus, z indicates the species population ratio of an optimal ensemble containing two radicals to the simulated one.26 Since the scale factor requires an update at each KMC iteration, the related calculations must be embedded prior to the reaction selection shown in the Figure 7 algorithm. Table 4 summarizes the results of the series of KMC simulations with the scaling factor implemented for both the “constant control volume” cases (no s correction factor) reported in Table 2, and the “constant simulation volume” (with s correction factor) cases of Table 3. Table 4 contains results only for S ≤ 20, as the larger control volumes associated with S > 20 do not require a “confined space” correction for radical concentrations (and thus reaction rates). The scaling method successfully suppresses the errors by resolving the compartmentalization of radicals; similar patterns in radical number are observable for S=1 with scaling compared to the ensemble where the scaling factor was not applied (see supplementary information – Figure S1-S6). Indeed, the errors associated with applying both the correction factor s and the scaling factor z are comparable to those for the corresponding cases for which only the scaling factor was implemented, with the small differences attributable to random fluctuations around the expected profiles rather than a constant bias. However, the solution is accelerated effectively through introducing s to maintain the simulation volume constant, reducing the computational times by approximately 30-40%. Figure 8 provides a graphical comparison of the different techniques; while the scaling methodology substantially reduces the simulation error, its combination with the volume correction volume accelerates the simulation, as seen by the lowered slope in the plot of simulation time against sampling number. In fact, utilizing both methods simultaneously further accelerate the algorithm to converge the stochastic data toward the true values.



Page 20 of 29

Table 4. The total computational time requirements and the total error percentages calculated based on whether only the scaling method z is applied or accompanied by implementation of the correction factor s .

z

s

3

AKLMN × 10?

Cumulative Average Errors* _à &_`R F

20 1.5428 1.5428→0.6037 16 1.2342 1.2342→0.4830 12 0.9256 0.9256→0.3624 10 0.7714 0.7714→0.3021 9 0.6942 0.6942→0.2719 8 0.6171 0.6171→0.2417 7 0.5399 0.5399→0.2116 6 0.4628 0.4628→0.1814 5 0.3857 0.3857→0.1512 4 0.3085 0.3085→0.1211 3 0.2314 0.2314→0.0909 2 0.1542 0.1542→0.0607 1 0.0771 0.0771→0.0306

m kLnopN 6m k 6

∗ jk = l

m kLnopN 6

a

_`b &_`c

2.00 1.99 2.13 2.07 2.10 2.01 1.95 2.01 1.84 2.11 2.03 2.14 1.79 2.03 1.85 2.08 1.86 2.06 2.15 1.76 1.85 1.87 1.84 1.19 1.83 0.80

F

b

_de &_da &_dR &_dfg

0.52 0.66 0.52 0.49 0.40 0.71 0.49 0.37 0.46 0.42 0.49 0.52 0.61 0.36 0.31 0.49 0.38 0.30 0.93 0.63 0.92 0.82 0.71 0.68 1.13 1.93

:

0.45 0.65 0.38 0.56 0.70 0.51 0.54 0.50 0.46 0.30 0.58 0.43 0.60 0.30 0.14 0.30 0.75 0.30 0.48 0.57 0.80 0.81 0.54 0.60 0.53 1.76

c

h,Z ∗∗ ;,hH= 21 12 16 10 12 7.4 10 5.9 9.3 5.4 8.2 4.8 7.1 4.2 6.2 3.6 5.2 3.0 4.1 2.4 3.1 1.9 2.1 1.2 1.1 0.6




Page 21 of 29

55

55

50

Without Scaling/Without Correction Without Scaling/Without Correction Factor Factor Without Scaling/With Correction Without Scaling/With Correction Factor Factor With Scaling/With Correction With Scaling/With Correction Factor Factor With Scaling/Without Correction With Scaling/Without Correction Factor Factor

45 40

50 45 40

35

35

30

30

25

25

20

20

15

15

10

10

5

5

0

Time (min)

Total Ave. Cumulative Error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0 1

6

11

16

21

26

31

36

41

46

51

Sampling Number (S)

Figure 8. Mean time-averaged errors (left axis) and computational times (right axis) for KMC simulation of semi-batch copolymerization. The total errors are calculated based on the mean of the cumulative average errors reported in Table 2 (without scaling and correction factors), Table 3 (without scaling, with correction factors), and Table 4 (with scaling, with and without correction factors).

3.5 An Optimal Solution with Time-dependent Control and Simulation Volumes Although the error is small, the results attributed to S=1 (Table 4 entry with scaling and correction factor applied) is distinguishable in the initial stage of the reaction compared to the “true” solution calculated for S=50, in terms of the mole fraction of chains containing discrete numbers of GMA functional units (Figure 9). The error can be attributed to the randomness of the stochastic technique, meaning that the population of polymer chains formed in the extremely small simulation volume is insufficient to represent the macroscopic system. A classic example of this type of error is the stochastic estimation of the heads/tails probability of a coin, as the coin must be repeatedly flipped for a large number of trials to reach an accurate value.39 Similarly, an accurate estimation of the area underneath of a curve by MC sampling requires a large number of random points.27 Thus, even with a precise probability axis, a sufficient number of chains must be produced in the simulation to converge the results toward a satisfactory precision. As mentioned earlier, the concentrations of species are very low in the start-up of the semi-batch process operated under the - 21 ACS Paragon Plus Environment


starved-feed policy with only solvent present initially. As this stage of the process demands a much larger control volume than at later times in order to capture the details correctly, we introduce a procedure that initializes the simulation with a larger ensemble that decreases in size with time to lessen the computational cost.

×10-17

0.8

4.0

Sampling Number=[50-1] Sampling Number=1

0.7

Cumulative Mole Fraction

3.5 3.0

Volume (L)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 29

2.5 2.0

Sampling Number=50

0.6 0.5

F=0 0.4 0.3 0.2

1.5

0.1

1.0

0.0

F=1

F>2 F=2

0

50

100

150

0.5

200

Simulation Volume Control Volume

0.0 0

25

50

75

100

125

150

175

200

225

Time (min)

Figure 9. Variation in control and simulation volumes with time for “Sampling Number=[501]” case (main plot). Cumulative mole fractions with respect to the number of GMA group denoted as F for this optimized case compared to simulations with S=50 (Table 2) and with S=1 (Table 4 with correction factor) (inset plot).

To resolve the initial discrepancies associated with the case where S=1, the initial sampling number is set to 50. The control volume is then gradually reduced to approximately one-fiftieth of its initial value from the first minute after the reaction starts to the time corresponding to the maximum concentration of the initiator in the reactor; i.e., at 35 min (see supplementary information – Figure S7-S9). Analogous to the correction factor, a constant factor is applied to impose a g~1/50 reduction in the ensemble size at one-minute intervals, reaching 8TLM\N at 35 min from the beginning of the copolymerization. As this reduction occurs, the simulation volume is also adjusted,


Page 23 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


always taking into account the fresh feed added to the process (Figure 9). Afterward, the simulation proceeds normally according to the flowchart shown in Figure 7, with a constant simulation volume and a decreasing control volume. To maintain the reaction probabilities precisely, the scaling factor z is applied to the entire simulation. The accuracy is improved by adopting this strategy, as seen in the prediction of the distribution of GMA units in the early stages of the semi-batch reaction (Figure 9 inset plot); the cumulative average error over the initial 35 minutes is reduced from 5.2% to 0.3%. While the accuracy of the “true” (S=50) solution is approached, the simulation time is greatly reduced, from 50 to 2 min, as the majority of the calculations are performed with S=1. This technique of manipulating the control volume, while developed specifically for this example, can be adapted to accurately represent other situations for which radical concentrations are particularly low. As mentioned earlier, the concept of binary trees is used to simultaneously store and analyze the stochastic data. Since the systematic volume reduction occurs at every minute till the 35th minute of the reaction, the simulation exhibits thirty-six distinguishable ensembles having different initial volumes. For instance, the chains produced within the first minute of the reaction cannot be exposed to the population of chains generated within the second minute. To derive a cumulative property, the volume ratio of the ensembles should be considered as a coefficient to weight the corresponding property derived from a different ensemble. The postprocessing of the explicit chain sequences can be also achieved similar to this concept. In this work, however, the results are obtained in parallel with the KMC simulation and directly from the storage scheme shown graphically in Figure 2.

4. Conclusions Investigation of semi-batch starved-feed systems through stochastic simulations provides the opportunity of improving not only the quality of materials but the cost-efficiency of this common industrial approach for producing solvent-borne automotive coatings. As these polymers are often functionalized for end-use applications, knowledge of the average composition and molecular weights of the product is often not sufficient, as the distribution of the reactive moieties among the



chains affects performance. Minimizing the fraction of functional monomer used in the production of low molecular-weight resins for solvent-borne coatings is a challenge that can be tackled via accelerated stochastic modeling approaches. Since this opportunity is strongly connected with the computational time of these simulations, development of acceleration techniques is pivotal. The improved KMC solution developed in this work is demonstrated through consideration of a previously published case-study describing the radical copolymerization of glycidyl methacrylate (GMA) and butyl methacrylate in which the average number of GMA units per chain is unity. The solution provides equivalent accuracy for computation of the complete copolymer composition distribution along with the standard model outputs (i.e., average MW, monomer concentration profiles, etc.) in a fraction of the simulation time previously required. Coupled with the utilization of binary trees for storage of output, an improved implementation of the Gillespie’s algorithm for accelerating a semi-batch operation was established stepwise to reduce the computational cost. While the ensemble expands naturally due to the feeding process along with the changes in the density of the reaction media, the developed approach resizes the ensemble to accelerate the solution via maintaining the simulation volume at a constant level. In parallel, a scaling factor is utilized to prevent the simulation from encountering the effects of radical compartmentalization. The resulting description of the explicit sequence of chains is calculated at a time not much greater than that of deterministic models, which cannot supply this detailed information. Thus, the accelerated KMC methodology provides a solid foundation for optimization of the solvent-borne automotive coatings with the desired population of polymer chains synthesized under industrial conditions. As the next step towards this goal, the kinetic scheme utilized herein will be extended to include complexities such as penultimate copolymerization kinetics, depropagation, and reactions involving the formation and consumption of mid-chain radicals.

5. Supporting Information Further documentation of the results can be found in the supporting information.


Page 24 of 29

Page 25 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6. Acknowledgments The authors would like to thank Axalta Coating Systems and the Natural Sciences and Engineering Research Council of Canada for financial support.

7. References (1)

Nesvadba, P. Radical Polymerization in Industry. In Encyclopedia of Radicals in Chemistry, Biology and Materials; John Wiley & Sons, Ltd: Chichester, UK, 2012; pp 1–36.

(2)

Löffler, M.; Morschhäuser, R.; Hornung, M. Surfactant-Free Cosmetic, Dermatological and Pharmaceutical Compositions. US Pat. 8062630B2 2004.

(3)

Nguyen-Kim, S.; Hoessel, P.; Mueller, G. Cosmetic Product Comprising at Least One WaterSoluble Copolymer Which Contains (Meth)Acrylamide Units. US Pat. 8747824B2 2014.

(4)

Löffler, M.; Morschhäuser, R. Surfactant-Free Cosmetic, Dermatological and Pharmaceutical Agents. US Pat. 8062630B2 2011.

(5)

D’hooge, D. R.; Van Steenberge, P. H. M.; Reyniers, M.-F.; Marin, G. B. The Strength of Multi-Scale Modeling to Unveil the Complexity of Radical Polymerization. Prog. Polym. Sci. 2016, 58, 59–89.

(6)

Mastan, E.; Li, X.; Zhu, S. Modeling and Theoretical Development in Controlled Radical Polymerization. Prog. Polym. Sci. 2015, 45, 71–101.

(7)

Grady, M. C.; Simonsick, W. J.; Hutchinson, R. A. Studies of Higher Temperature Polymerization of N-Butyl Acryalte. Macromol. Symp. 2002, 182, 149–168.

(8)

Liang, K.; Hutchinson, R. A. Solvent Effects in Semibatch Free Radical Copolymerization of 2-Hydroxyethyl Methacrylate and Styrene at High Temperatures. Macromol. Symp. 2013, 325–326 (1), 203–212.

(9)

Wei, W.; Hutchinson, R. A. High Temperature Semibatch Free Radical Copolymerization of Dodecyl Methacrylate and Styrene. Macromol. Symp. 2008, 261 (1), 64–73.

(10)

Li, D.; Grady, M. C.; Hutchinson, R. A. High-Temperature Semibatch Free Radical Copolymerization of Butyl Methacrylate and Butyl Acrylate. Ind. Eng. Chem. Res. 2005, 44



(8), 2506–2517. (11)

Liang, K.; Rooney, T. R.; Hutchinson, R. A. Solvent Effects on Kinetics of 2-Hydroxyethyl Methacrylate Semibatch Radical Copolymerization. Ind. Eng. Chem. Res. 2014, 53 (18), 7296–7304.

(12)

Cunningham, M. F.; Hutchinson, R. Industrial Applications and Processes. In Handbook of Radical Polymerization; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2002; pp 333–359.

(13)

Mayo, F. R.; Walling, C. Copolymerization. Chem. Rev. 1950, 46 (2), 191–287.

(14)

Ito, K.; Yamashita, Y. Copolymer Composition and Microstructure. J. Polym. Sci. Part A Gen. Pap. 1965, 3 (6), 2165–2187.

(15)

Fukuda, T.; Kubo, K.; Ma, Y.-D. Kinetics of Free Radical Copolymerization. Prog. Polym. Sci. 1992, 17 (5), 875–916.

(16)

Ali Parsa, M.; Kozhan, I.; Wulkow, M.; Hutchinson, R. A. Modeling of Functional Group Distribution in Copolymerization: A Comparison of Deterministic and Stochastic Approaches. Macromol. Theory and Simul. 2014, 23 (3), 207–217.

(17)

Regatte, V. R.; Gao, H.; Konstantinov, I. A.; Arturo, S. G.; Broadbelt, L. J. Design of Copolymers Based on Sequence Distribution for a Targeted Molecular Weight and Conversion. Macromol. Theory and Simul. 2014, 23 (9), 564–574.

(18)

Saeb, M. R.; Mohammadi, Y.; Ahmadi, M.; Khorasani, M. M.; Stadler, F. J. A Monte CarloBased Feeding Policy for Tailoring Microstructure of Copolymer Chains: Reconsidering the Conventional Metallocene Catalyzed Polymerization of α-Olefins. Chem. Eng. J. 2015, 274, 169–180.

(19)

D’hooge, D. R. In Silico Tracking of Individual Species Accelerating Progress in Macromolecular Engineering and Design. Macromol. Rapid Commun. 2018, 1800057, 1800057.

(20)

Kryven, I.; Zhao, Y. R.; McAuley, K. B.; Iedema, P. A Deterministic Model for Positional Gradients in Copolymers. Chem. Eng. Sci. 2018, 177, 491–500.

(21)

Wang, L.; Broadbelt, L. J. Tracking Explicit Chain Sequence in Kinetic Monte Carlo Simulations. Macromol. Theory and Simul. 2011, 20 (1), 54–64. - 26 ACS Paragon Plus Environment

Page 26 of 29

Page 27 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(22)

Hamzehlou, S.; Reyes, Y.; Leiza, J. R. A New Insight into the Formation of Polymer Networks: A Kinetic Monte Carlo Simulation of the Cross-Linking Polymerization of S/DVB. Macromolecules 2013, 46 (22), 9064–9073.

(23)

Marien, Y. W.; Van Steenberge, P. H. M.; Barner-Kowollik, C.; Reyniers, M. F.; Marin, G. B.; D’Hooge, D. R. Kinetic Monte Carlo Modeling Extracts Information on Chain Initiation and Termination from Complete PLP-SEC Traces. Macromolecules 2017, 50 (4), 1371–1385.

(24)

Van Steenberge, P. H. M.; Vandenbergh, J.; Reyniers, M. F.; Junkers, T.; D’hooge, D. R.; Marin, G. B. Kinetic Monte Carlo Generation of Complete Electron Spray Ionization Mass Spectra for Acrylate Macromonomer Synthesis. Macromolecules 2017, 50 (7), 2625–2636.

(25)

Tripathi, A. K.; Sundberg, D. C. A Hybrid Algorithm for Accurate and Efficient Monte Carlo Simulations of Free-Radical Polymerization Reactions. Macromol. Theory and Simul. 2015, 24 (1), 52–64.

(26)

Gao, H.; Oakley, L. H.; Konstantinov, I. A.; Arturo, S. G.; Broadbelt, L. J. Acceleration of Kinetic Monte Carlo Method for the Simulation of Free Radical Copolymerization through Scaling. Ind. Eng. Chem. Res. 2015, 54 (48), 11975–11985.

(27)

Brandão, A. L. T. T.; Soares, J. B. P. P.; Pinto, J. C.; Alberton, A. L. When Polymer Reaction Engineers Play Dice: Applications of Monte Carlo Models in PRE. Macromol. React. Eng. 2015, 9 (3), 141–185.

(28)

Van Steenberge, P. H. M.; D’hooge, D. R.; Reyniers, M.-F.; Marin, G. B. Improved Kinetic Monte Carlo Simulation of Chemical Composition-Chain Length Distributions in Polymerization Processes. Chem. Eng. Sci. 2014, 110, 185–199.

(29)

Chaffey-Millar, H.; Stewart, D.; Chakravarty, M. M. T.; Keller, G.; Barner-Kowollik, C. A Parallelised High Performance Monte Carlo Simulation Approach for Complex Polymerisation Kinetics. Macromol. Theory and Simul. 2007, 16 (6), 575–592.

(30)

Barsotti, R. J.; Lewin, L. A.; Scopazzi, C. Coating Compositions Containing Non-Aqueous Dispersed Polymers Having a High Glass Transition Temperature. US Pat. 5763528 1998.

(31)

Wang, W.; Hutchinson, R. A.; Grady, M. C. Study of Butyl Methacrylate Depropagation Behavior Using Batch Experiments in Combination with Modeling. Ind. Eng. Chem. Res. - 27 ACS Paragon Plus Environment


2009, 48 (10), 4810–4816. (32)

Gillespie, D. T. Stochastic Simulation of Chemical Kinetics. Annu. Rev. Phys. Chem. 2007, 58 (1), 35–55.

(33)

Al-Harthi, M. A. MATLAB Programming of Polymerization Processes Using Monte Carlo Techniques. In Applications of Monte Carlo Method in Science and Engineering; InTech, 2011; pp 841–856.

(34)

D’Hooge, D. R.; Van Steenberge, P. H. M.; Derboven, P.; Reyniers, M. F.; Marin, G. B. Model-Based Design of the Polymer Microstructure: Bridging the Gap between Polymer Chemistry and Engineering. Polym. Chem. 2015, 6 (40), 7081–7096.

(35)

Gao, H.; Broadbelt, L. J.; Konstantinov, I. A.; Arturo, S. G. Acceleration of Kinetic Monte Carlo Simulations of Free Radical Copolymerization: A Hybrid Approach with Scaling. AIChE J. 2017, 63 (9), 4013–4021.

(36)

Neuhaus, E.; Herrmann, T.; Vittorias, I.; Lilge, D.; Mannebach, G.; Gonioukh, A.; Busch, M. Modeling the Polymeric Microstructure of LDPE in Tubular and Autoclave Reactors with a Coupled Deterministic and Stochastic Simulation Approach. Macromol. Theory and Simul. 2014, 23 (7), 415–428.

(37)

Schütte, C.; Wulkow, M. A Hybrid Galerkin-Monte-Carlo Approach to Higher-Dimensional Population Balances in Polymerization Kinetics. Macromol. React. Eng. 2010, 4 (9–10), 562– 577.

(38)

Demirel Özçam, D.; Teymour, F. Chain-by-Chain Monte Carlo Simulation: A Novel Hybrid Method for Modeling Polymerization. Part I. Linear Controlled Radical Polymerization Systems. Macromol. React. Eng. 2017, 11 (1), 1600042.

(39)

Sugita, H. Mathematics of Coin Tossing Tossing. In Probability and Random Number; World Scientific Publishing Co. Pte. Ltd.: Singapore, 2017; pp 1–22.


Page 28 of 29

Page 29 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical Abstract for Table of Contents:

A graphical abstract contrasting the conventional strategy in modeling semi-batch (A) and the time-dependent control volume with adjusted simulation volume (B).

Keywords: radical copolymerization, copolymer composition distribution, Kinetic Monte Carlo simulation, starved-feed semi-batch


Modeling the Distribution of Functional Groups in Semibatch Radical

Recommend Documents