Journal Information

Article Information


A Kinetic Monte Carlo Simulation of Individual Site Type of Ethylene and α-Olefins Polymerization


Abstract

The aim of this work is to study Monte Carlo simulation of ethylene (co)polymerization over Ziegler-Natta catalyst as investigated by Chen et al. The results revealed that the Monte Carlo simulation was similar to sum square error (SSE) model to prediction of stage II and III of polymerization. In the case of activation stage (stage I) both model had slightly deviation from experimental results. The modeling results demonstrated that in homopolymerization, SSE was superior to predict polymerization rate in current stage while for copolymerization, Monte Carlo had preferable prediction. The Monte Carlo simulation approved the SSE results to determine role of each site in total polymerization rate and revealed that homopolymerization rate changed from site to site and order of center was different compared to copolymerization. The polymer yield was reduced by addition of hydrogen amount however there was no specific effect on uptake curve which was predicted by Monte Carlo simulation with good accuracy. In the case of copolymerization it was evolved that monomer chain length and monomer concentration influenced the rate of polymerization as rate of polymerization reduced from 1-hexene to 1-octene and increased when monomer concentration proliferate.


Expand AllCollapse All

INTRODUCTION

Polyethylene grades have developed over years as the dominant polymers used in many industries such as automotive, pipes, containers, films and electrical conduits.1,2 Ziegler-Natta catalysts particularly those based on titanium, make the majority of commercial polyethylene.35 It is well known that these catalysts have multiple active centers with difference kinetic features.6,7 Some centers are stable and make a higher contribution to the final polymer. Moreover, each center produces a polymer with different molecular weight; they all follow the Shultz-Flory distribution with PDI=2. In addition, various centers have dissimilar susceptibilities to (co)polymerize different monomers.8 As well as this, their reactions toward impurities such as oxygen or water are entirely discrepant.

The polymerization kinetics of Ziegler-Natta catalysts with simplified mathematical models which quantify overall polymerization rates and average molecular weights have been extensively described in the literature.9,10 Also the mechanism of polymerization with these catalysts has been great deal of researcher interests.

Kissin and his coworkers11 derived reaction rate profile for ethylene slurry polymerization and studied the effect of deviation reactions on the reaction rate. Ethylene polymerization kinetics by moment equation modeling to study the effect of different active centers on homopolymerization kinetics was perused by Kalajahi et al.12 Gemoets and his coworkers13 used polymerization kinetics model based on lumping two catalyst sites to predict ethylene polymerization rates and polymer average molecular weights. Study and characterization of polymerization kinetics in the presence of diffusion phenomena were investigated by Ray et al.14 Alshaiban and Soares15 investigated hydrogen and external donors effect on propylene polymerization kinetics.

To the best of our knowledge few publications delineate methods to estimate polymerization kinetics parameters for each site type on the multisite catalyst.16,17 However they did not estimate the minimum number of active site types required to simultaneously clarify instantaneous polymerization kinetics, cumulative polymer yields as well as molecular weight distribution (MWD) using a fundamental mechanism for coordination polymerization.8,18

Chen and coworkers19 introduced method which identifies the minimum number of active site types required to simultaneously clarify instantaneous polymerization kinetics and MWDs evolution for ethylene and α-olefin polymerizations with multisite catalysts. They quantify apparent site activation, monomer/comonomer propagation and site deactivation rate constants for each site type and estimate the MWDs of polymer populations made on each site type.

It should be noted that the statistical nature of chain growth and chain-terminating reaction makes probability theory, particularly Monte Carlo simulation methods, and strong tools in this terrain.20,21 The Monte Carlo simulation methods are able to determine all distributions as well as their mean values due to storing the whole information of (co)polymerization chains while the reaction proceeds.22 Therewith, Monte Carlo simulation method requires only reaction rate constants to simulate polymerization reactions. In addition, Monte Carlo can remark the composition drift and azeotropic properties of copolymers with good accuracy.23,24

In the case of ethylene polymerization kinetics a few researches has been carried out to consider the contribution of each active center to polymerization and determine the composition of the polymer made by each center. Therefore, the aim of this study is to develop appropriate simulation utilizing the Monte Carlo method to investigate the cumulative yield and instantaneous polymerization rate using kinetic constants reported by Chen et al.19

POLYMERIZATION KINETICS

Validity and practicality of kinetic model proposed for polymerization reactions strongly depend on the conception of all the phenomena involved in catalytic polymerization reactions. In the case of Ziegler-Natta catalysts due to multiple centers and effect of diffusion barriers, the model should be more powerful to properly describe the polymerization kinetics. It is well known that broad molecular weight distribution of products synthesized by Ziegler-Natta catalysts is less influenced by diffusion barriers.11 The kinetic mechanisms are based of Chen et al. report.19

THEORY AND SIMULATION

Monte Carlo algorithm

Monte Carlo simulation methods are based on the use of random numbers to sample the variable space using a probability distribution followed by the selection of an event.25 In present study, the principle for the simulation of ethylene polymerization kinetics was developed on the basis of Gillespies’ algorithm.26 Accordingly, the simulation volume, V, was assumed to be divided homogeneously between the reactants as well as simulation volume consists of 106 catalyst precursor molecules and the number of other components in calculated based on catalyst precursor molecules:

(1)
V = N C A T C A T N a v

(2)
N X = X N a v V

where V is simulation volume, [CAT] concentration of catalyst, NCAT number of catalyst precursor, Nav Avogadro’ number, NX number of X molecules, and [X] the concentration of X species. In addition, because the Monte Carlo simulation method converse with a number of molecules, stochastic reaction rates must replace macroscopic ones. While there are L reactions in the simulation, the probability of incidence of reaction l (Pl) is given by27,28

(3)
P l = a l l = 1 l = L a l

where al, the stochastic reaction rate for reaction l, is defined by

(4)
a = h × c

h represents the number of reactants and c stands for stochastic rate constants and correlates with ordinary reaction rates. Equations 5 and 6 demonstrate first and second order of reaction, respectively:

(5)
c = k

(6)
c = k V N a v

Microscopic elementary oligomerization and copolymerization reactions raised discretely and stochastically through M reaction channels and an event was designated in a given time interval (t, t+dt) from uniformly distributed random number in a unit interval, according to following relationships:

(7)
v = 1 μ 1 P v < r 1 v = 1 μ P v

(8)
d t = 1 v = 1 M R v ln 1 r 2

where r1 and r2 are two random numbers, µ is the number of the selected reaction channel, Pv and Rv are the reaction rate probability and rate of reaction v, respectively, while dt is the time interval between two successive reactions.

The simulation program was written in Matlab. A schematic demonstration of the Monte Carlo method algorithm is shown in Scheme 1.

Scheme1.

Monte Carlo algorithm for simulating (co)polymerization of ethylene.

jkcs-62-191-f010.tif

Simulation assumptions

1. It is supposed that the simulation volume is homogeneous and elementary reactions are not controlled by diffusion.

2. The reaction temperature and ethylene pressure were kept constant during polymerization process.

3. It is assumed that all the living chain in the system have and equal selection probability.

4. The model using in copolymerization is terminal model.

RESULTS AND DISCUSSION

Ethylene homopolymerization

Influence of H2 pressure on the variation of polymerization rate with time is shown in Fig. 1. By increasing the concentration of hydrogen in the reactor, the polymerization rate slightly decreased, but did not affect the shape of the ethylene uptake curve.

As it can be seen from Fig. 1 the Monte Carlo simulation can predict well all three stage of polymerization(i.e. activation, balance between activation and deactivation, and deactivation). A slight deviation was observed from model fitting and experimental part as well as rate of polymerization for monomer consumptionand synthesized polymer first increased and then dropped due to site activation and site deactivation respectively.

Figure1.

Effect of H2 on ethylene polymerization rate.

jkcs-62-191-f001.tif

Fig. 2 and 3 demonstrate that the model fits well the instantaneous ethylene polymerization rates andusing 5 site type noticed by MWD deconvolution. It can be seen that the highest deviation occur at short polymerization time due to fluctuation of temperature and pressure which is slightly more at the beginning of the polymerizations. These phenomena can be ascribed to the increasing of the termination rate constant (kd) values and decrease the propagation rate constant (kp) values, which agrees with the observation that H2 usually decreases polyethylene yields with Ziegler-Natta catalysts.

Figure2.

Predicted and measured instantaneous polymerization rates for ethylene at 0.45 bar H2; a) site I, b) site II, c) site III, d) site IV, e) site V, f) total

jkcs-62-191-f002.tif
Figure3.

Predicted and measured instantaneous polymerization rates for ethylene at 0.7 bar H2; a) site I, b) site II, c) site III, d) site IV, e) site V, f) total.

jkcs-62-191-f003.tif

It has been observed that Monte Carlo predicts kinetic behavior of each site with good accuracy using kinetic rate constants from Table 1. Also great prediction was observed specially after 6 min (i.e. activation step). As we can see, site IV has highest activation and then III, II, V, I respectively. However, plateau curve was higher for site I and V than other sites, meaning that step of between activation and deactivation is high for these two sites adoption of this result with kinetic equations for each site, it can be concluded that site IV and III have the largest contribution for polymer production.

Table1.

Apparent rate constants for ethylene homopolymerization19

PH2[bar] Ka[min-1] Kp[min-1] kd[min-1]
Site 1 0.45 0.8538 0.002800 0.01938
0.70 0.9935 0.002750 0.02307
Site 2 0.45 0.3975 0.01249 0.02532
0.70 0.5087 0.01206 0.02677
Site 3 0.45 0.3165 0.02012 0.01416
0.70 0.3655 0.01987 0.01442
Site 4 0.45 0.2274 0.02255 0.01154
0.70 0.2283 0.01941 0.01172
Site 5 0.45 0.1706 0.007330 0.001
0.70 0.2572 0.005932 0.001

The predicted and measured polymer yields are shown in Fig. 4. By increase in hydrogen pressure polymer yields decreased as a result number molecular weight (Mn) decreased but interestingly increased steadily with polymerization time as well as uptake curve was not influenced.

Figure4.

Predicted and measured polyethylene yields.

jkcs-62-191-f004.tif

The ka values increased for all sites as more H2 is employed, maybe because shorter polyethylene chains produced in the presence of H2 fragment the catalyst particles more effectively, better exposing active sites for polymerization. Introducing of hydrogen to the reaction exceedingly reduces molecular weight and due to large amounts of transfer to hydrogen and deactivation reactions.

It should be considered that balance between the rates for site activation, ethylene propagation, and site deactivation for each site type determines how the relative masses of polymer populations differs with polymerization time.

Ethylene/α-olefin copolymerization

It has been observed that the Monte Carlo simulation method was employed to predict kinetic behavior of each site with using kinetic rate constants from Table 2.

Table2.

Apparent rate constants for ethylene/α-olefin Copolymerization19

α-Olefin type α-Olefin[g] Ka[min-1] K ^ p min 1 kd[min-1]
Site 1 1-hexane 2.0 1.4917 0.01469 0.2435
1-hexane 4.0 2.0519 0.02023 0.2503
1-hexane 6.0 3.0746 0.02500 0.2594
1-octene 2.0 0.8469 0.008074 0.2215
Site 2 1-hexane 2.0 0.8568 0.03014 0.1051
1-hexane 4.0 1.0728 0.03502 0.1460
1-hexane 6.0 1.2292 0.04057 0.1808
1-octene 2.0 0.5673 0.02650 0.1040
Site 3 1-hexane 2.0 2.5940 0.03786 0.05112
1-hexane 4.0 4.7167 0.03713 0.06416
1-hexane 6.0 9.1820 0.03678 0.07962
1-octene 2.0 1.5282 0.03173 0.04090
Site 4 1-hexane 2.0 2.7965 0.02868 0.02750
1-hexane 4.0 4.7854 0.02512 0.02849
1-hexane 6.0 6.5020 0.02259 0.03099
1-octene 2.0 1.1888 0.02475 0.02325
Site 5 1-hexane 2.0 1.5590 0.005954 0.001
1-hexane 4.0 2.4329 0.004834 0.001
1-hexane 6.0 3.5810 0.003512 0.001
1-octene 2.0 1.0853 0.005818 0.001

The model predictions with measured values for the instantaneous ethylene/α-olefin copolymerization rate are compared in Fig. 5(a,b,c,d) Model can predict polymerization times exceeding about 10 min with good accuracy. But are less accurate at the very beginning of copolymerization, likely due to temperature fluctuation instantly after the catalyst is injected in the reactor. The Monte Carlo Simulation had better prediction than model for activation step. As expected, Mn declined when more 1-hexene ispresent in the reactor due to higher transfer rates tocomonomer. Substitution 1-hexene with 1-octene leadsto copolymers with higher Mn, likely because 1-octene is not as reactive (for propagation and chain transfer) as1-hexene. Also, from Fig. 4, it can be observed that Monte Carlo simulation can predict copolymerization behavior and each site contribution with excellent accuracy.

Figure5.

(a) Predicted and measured instantaneous ethylene/α-olefin copolymerization (2g 1-hexene); a) site I, b) site II, c) site III, d) site IV, e) site V, f) total

jkcs-62-191-f005.tif
Figure5.

(continued) (b) Predicted and measured instantaneous ethylene/α-olefin copolymerization (4g 1-hexene); a) site I, b) site II, c) site III, d) site IV, e) site V, f) total

jkcs-62-191-f006.tif
Figure5.

(continued) (c) Predicted and measured instantaneous ethylene/α-olefin copolymerization (6g 1-hexene); a) site I, b) site II, c) site III, d) site IV, e) site V, f) total

jkcs-62-191-f007.tif
Figure5.

(continued) (d) Predicted and measured instantaneous ethylene/α-olefin copolymerization (2g 1-octene); a) site I, b) site II, c) site III, d) site IV, e) site V, f) total

jkcs-62-191-f008.tif

The Ka and kd values for each site type increased atdifferent ratios when more 1-hexene is added to thereactor. Also, the Kp values for the low Mn sites (1 and 2) increased with addition of 1-hexene. However, the Kp values for the high Mn sites (3 to 5) decrease withincreasing 1-hexene concentration. It should be noticed that low Mn sites in Ziegler-Natta catalysts have higher reactivity ratiostoward comonomer incorporation; therefore, it is attended that they will suffer the strongest comonomereffect. Finally, the Ka, Kp, and kd values on five active sites for 1-octene were lower than that for 1-hexene, likely because 1-octene has a weaker comonomer effect than 1-hexene.

Ethylene uptake curves for ethylene/α-olefin copolymerization are illustrated in Fig. 6. As depicted in Fig. 6, type and concentration of α-olefin does not affect ethylene uptake rates significantly after the first 10 min of polymerization. Instantly, after catalyst injection, the rate proliferates slightly faster for copolymerization with higher amount of 1-hexene, however this discrepancyswiftly became negligible after a few minutes, suggesting that is a chemical effect rather than physical effect. Moreover, by addition of 1-octene as comonomer, the ethylene uptake rates at the beginning of the polymerization rate are lower, and augment steadily as the concentration of 1-hexene is raised.

Figure6.

Ethylene uptake curves for ethylene/α-olefin copolymerization; a) 2g 1-hexene, b) 4g 1-hexene, c) 6g 1hexene d) 2g 1-butene.

jkcs-62-191-f009.tif

As we can see from this figure Monte Carlo are able to excellent prediction of activation step rather than SSE modelthis is maybe due to nature of statistical method which consider all reaction parts and phenomena during the reaction. In addition to this, the Mn of thepolymer populations is not substantially influenced by theconcentration of 1-hexene in the reactor, which is aninteresting characteristic of this particular Ziegler-Nattacatalyst because it allows near independent control ofmolecular weight and comonomer molar fraction in the copolymer. In general, according to the results simulation error was less than 0.1 due to the fact that kinetic constant was obtained from deconvolution.

It should be noted that structure complexity of Ziegler-Natta catalyst leads to produce polymer with different structure and final properties. Therefore, kinetic constant is so substantial parameters relating to inherent properties of catalyst. Khorasani et al.29 and Kalajahi et al.30 studied Monte Carlo simulation of tandem polymerization of ethylene and Ziegler Natta catalysts, respectively. They used kinetic constants from the articles for their polymerization. However, in this work due to using deconvoluted kinetic constant a good accuracy for theoretical and experimental results was observed. What’s more, in this work the modeling accuracy was compared with experimental data in polymerization and copolymerization, therewith using this method varing structure can be predict. Moreover, by using SSE method the kinetic and polymer processing can be predict which is very substantial for reducing process costs and developing new grades of polymers.

CONCLUSION

In the current work, a Monte Carlo simulation method was developed to study ethylene polymerization kinetics. According to the result obtained, the Monte Carlo simulation predicted stage II and III of polymerization similar to SSE model. Both model had slight deviation from experimental data for stage I. however in homopolymerization SSE had better quality to predict stage I than Monte Carlo while in copolymerization Monte Carlo had better prediction.

Homopolymerization rate of different active center changed in the order of center IV> center III> center V> center II> center I while in copolymerization, the rate changed in the order of center IV> center V> center III> center II> center I.

Moreover, although hydrogen destroyed catalyst activation centers and consequently reduced catalyst yield, it was no specific effect on the uptake curve. Acting as a main transfer agent, hydrogen also reduced polymer molecular weight.

Monte Carlo also predicted polymer yield with excellent accuracy, by addition of hydrogen amount polymer yield decreased, however there was no specific influence on uptake curve.

Acknowledgements

Publication cost of this paper was supported by the Korean Chemical Society.

References

1. 

A. L. Andrady Plastics and the EnvironmentJohn Wiley & Sons2003

2. 

A. Peacock Handbook of Polyethylene: Structures: Properties, and ApplicationsCRC Press2000

3. 

S. M. Ghafelebashi Zarand S. M. M. Mortazavi Journal of Petroleum Science and Technology2012212

4. 

M. Salami-Kalajahi V. Haddadi-Asl M. Najafi S. M. Ghafelebashi Zarand e-Polymers20080041

5. 

Y. Liu B. Zhang Z. Fu Z. Fan Macromol. Res.20171

6. 

Y. V. Kissin R. I. Mink T. E. Nowlin A. J. Brandolini Top. Catal.1999769 [CrossRef]

7. 

Y. V. Kissin J. Polym. Sci. Part A Polym. Chem.2003411745 [CrossRef]

8. 

Y. V. Kissin J. Polym. Sci. Part A Polym. Chem.199533227 [CrossRef]

9. 

V. B. Skomorokhov V. A. Zakharov V. A. Kirillov Macromol. Chem. Phys.19961971615 [CrossRef]

10. 

A. A. Barabanov V. A. Zakharov V. V. Sukulova J. Organomet. Chem.2015798292 [CrossRef]

11. 

Y. V. Kissin F. M. Mirabella C. C. Meverden J. Polym. Sci. Part A Polym. Chem.2005434351 [CrossRef]

12. 

M. Salami-Kalajahi V. Haddadi-Asl M. Najafi S. M. Ghafelebashi Zarand e-Polymers2008829

13. 

F. Gemoets M. Zhang T. W. Karjala B. W. S. Kolthammer Macromol. React. Eng.20104109 [CrossRef]

14. 

R. A. Hutchinson C. M. Chen W. H. Ray J. Appl. Polym. Sci.1992441389 [CrossRef]

15. 

A. Alshaiban J. B. P. Soares Macromol. React. Eng.20126265 [CrossRef]

16. 

K. Z. Yao B. M. Shaw B. Kou K. B. McAuley D. W. Bacon Polym. React. Eng.200311563 [CrossRef]

17. 

K. Chen S. Mehdiabadi B. Liu J. B. P. Soares Macromol. React. Eng.2016

18. 

Y. V. Kissin R. I. Mink T. E. Nowlin J. Polym. Sci. Part A Polym. Chem.1999374255 [CrossRef]

19. 

K. Chen S. Mehdiabadi B. Liu J. B. P. Soares Macromol. React. Eng.20161

20. 

W. Bruns I. Motoc K. F. O’Driscoll Monte Carlo Applications in Polymer ScienceSpringer Science & Business Media201227

21. 

G. G. Lowry Markov Chains and Monte Carlo Calculations in Polymer ScienceMarcel Dekker1970

22. 

J. Shao W. Tang R. Xia X. Feng P. Chen J. Qian C. Song Macromol. Res.2015231042 [CrossRef]

23. 

D. Beigzadeh Macromol. Theory Simul.200312174 [CrossRef]

24. 

J. He H. Zhang J. Chen Y. Yang Macromolecules1997308010 [CrossRef]

25. 

M. Fluendly Markov Chains and Monte Carlo Calculations in Polymer ScienceMarcel DekkerNew York1970

26. 

D. T. Gillespie J. Comput. Phys.197622403 [CrossRef]

27. 

L. U. O. Zhenghon C. A. O. Zhikai S. U. Yaotang Chin. J. Chem. Eng.200614194 [CrossRef]

28. 

M. Al’Harthi J. B. P. Soares L. C. Simon Macromol. Mater. Eng.2006291993 [CrossRef]

29. 

M. M. Khorasani M. R. Saeb Y. Mohammadi M. Ahmadi Chem. Eng. Sci.2014111211 [CrossRef]

30. 

M. S. Kalajahi M. Najafi V. H. Asl Application of Monte Carlo simulation method to polymerization kinetics over Ziegler-Natta catalystsInt. J. Chem. Kinet.20094145 [CrossRef]