Computational and Theoretical Polymer Science 11 (2001) 89–94 www.elsevier.nl/locate/ctps
Photoinitiated bulk polymerization of furfuryl methacrylate. Stochastic simulation results using the new model CORUB J.A. Corzo E. a,*, N. Davidenko b, R. Alvarez b a
Chemistry and Chemical Engineering Department, Matanzas University, Matanzas, Cuba b Center of Biomaterials, Havana University, Havana, Cuba
Received 3 February 1999; received in revised form 8 December 1999; accepted 13 December 1999
Abstract A new stochastic model named CORUB has been used to do a phenomenological description of the photopolymerization of furfuryl methacrylate, with the following methodology. Statistical analysis of the model’s adequacy shows a good correlation between experimental data and those modelled stochastically. This model was developed using two stochastic variables that represent the random nature of the intermolecular encounters upon which the Collision Theory is sustained. The amounts of random numbers generated simulate time in a new way, based on the standard physical concept, with satisfactory results. Also, calculating the number of times that each individual step of the reaction mechanism takes place, a new possibility for carrying out the sensitivity analysis is presented. In addition to which, the analogy assumed between the frequency of the random number generated and the effective intermolecular collision allows us to model, exclusively and accurately enough, the effect of temperature as well as to estimate the apparent activation energy precisely. The algorithm of CORUB model was coded in Turbo Pascal 6.0, resulting in the CHEMOD-X program, which constitutes a software specifically elaborated as an integrative part of this work, allowing the model to be used. 䉷 2000 Elsevier Science Ltd. All rights reserved. Keywords: Stochastic modelling; Monte Carlo method; Photopolymerization
1. Introduction Recently, the polymerization of acryl-furanic monomers was studied through kinetic simulation for the first time [1]. In this study, a mechanism for the photopolymerization of furfuryl methacrylate (FM) initiated by AIBN was proposed and the energy parameters of the Arrhenius equation were estimated for each individual step. The authors applied suitable professional computer programs that incorporated integration methods of Ordinary Differential Equation Systems (ODES) together with different parameter estimation routines, to make possible a mathematical simulation and sensitivity analysis of the photopolymerization. The simulation of chemical reactions could also be carried out from a stochastic perspective, without considering the representation with ODES. The usual mathematical procedure for the calculations derived from the stochastic models in Chemistry is the universal Monte Carlo method [2–8]. In the present paper, the stochastic modelling of the photopolymerization of FM is drawn up, based on results obtained by Lange and co-workers [1]. The name of the
stochastic model used is CORUB [9–11]; this new model is based on the Collision Theory and takes into account the probability of the effective intermolecular encounter. In the model this probability depends of the relative quantities of molecules of the reactants for each one of the individual steps considered.
2. Methods 2.1. Brief description of the CORUB model The CORUB model carries out the simulation of chemical reactions by employing two discreet random variables denominated “effective collision” (e_co) and “reaction step” (re_s). A clear definition of these variables can be shown through their distribution tables [8] as follows: " # 1 2
1 e_co p1 p2 "
* Corresponding author. Tel.: ⫹53-52-261432; fax: ⫹53-52-253101. E-mail address:
[email protected] (J.A. Corzo E.).
re_s
1089-3156/01/$ - see front matter 䉷 2000 Elsevier Science Ltd. All rights reserved. PII: S1089-315 6(99)00093-8
1
2
3
…
z
q1
q2
q3
… qz
#
2
90
J.A. Corzo E. et al. / Computational and Theoretical Polymer Science 11 (2001) 89–94
Fig. 1. General flow chart for CORUB’s programs where c(i)0 is the initial molar concentration of i, c(i) the instantaneous molar concentration of i, conver the maximum conversion, Avog_const the constant related with the Avogrado number, Ni0 the initial quantity of molecules of i, Ni the instantaneous quantity of molecules of i, and NT0 the initial total number of molecules.
2.2. The stochastic variable e_co The variable e_co provides a new viewpoint for a phenomenological description because the probability of a reaction occurrence is given by p1. Then, given an elementary act of the type: aA ⫹ bB ⫹ cC Products
3
where a, b, c are the stoichiometric coefficients such that
a ⫹ b ⫹ c ⱕ 3: Then: p1
NA a
NB b
NC c
NT0 n
4
where n a ⫹ b ⫹ c is the order of the elementary act of reaction, Ni the number of molecules of the reactant i, and NT0 the initial total number of molecules. The value of e_co is determined by generating a random number (RND). This random number is uniformly distributed in [0,1], and when: 0 ⬍ RND ⱕ p1
5
CORUB assigns the value 1 to the variable e_co and otherwise the value 2. If e_co 1 then the elementary reaction act occurs and the Ni values are changed according to their corresponding stoichiometric coefficients, or else the reaction act does not occur. 2.3. The stochastic variable re_s The distribution table of the variable re_s, given by expression 2, shows that this variable could take z values corresponding to the number of steps in the mechanism. The step i has the probability of occurrence qi and is calculated only in terms of the comparative value of its kinetic constants (ki), as can be seen in expression 6. This position differs from previous considerations based on the use of their reaction rates [12–19]. ki qi X z
6 ki
i1
To decide which reaction step will occur, a new random
J.A. Corzo E. et al. / Computational and Theoretical Polymer Science 11 (2001) 89–94 Table 1 Values of the kinetic constants (Lmol ⫺1 s ⫺1) employed in the simulation Step k
e-2 802.6
e-3 752.7
e-4 6 × 105
e-5 621.8
e-6 3:97 × 10⫺2
Step k
e-7 1:001 × 104
e-8 1 × 104
e-9 7 × 105
e-10 998.8
e-11 1.00
91
constant that CORUB can perform, for the simulated chemical reaction. Thus, M 1 is associated with the minimum temperature necessary to activate the modelling reaction. 2.6. General algorithm and description of the computing program created
the elementary reaction act number g is then assumed (g is a natural number in the range 1 ⱕ g ⱕ z: Then, combining the stochastic simulation of the elementary acts, CORUB can reproduce phenomenologically the chemical reaction mechanisms [9–11]. This point provides an intuitive grasp of the simulation, even in the most complex case.
The CORUB method carries out the stochastic simulation uniquely. An immediate consequence is the necessary coding of the programs that execute the simulation, because programs for this purpose have not yet been developed. For the elaboration of the programs, the Turbo Pascal 6.0 language was used. The general algorithm of the programs created, in order to model the chemical reaction systems, can be followed, step by step, on a general flow chart. This chart appears in Fig. 1. In order to model the photopolymerization of FM initiated by AIBN, taking into account the mechanism that appears below [1], it is necessary to code the Pascal program CHEMOD-X. This program expresses the CORUB’s algorithm for the previously cited reaction.
2.4. The variable time
I ! 2R䡠
number must be generated. In this way, the value of re_s is calculated because when: gX ⫺1
qi ⬍ RND ⱕ
i1
g X
qi
7
i1
e⫺1
R䡠 ⫹ FM ! M1䡠
e⫺2
In the CORUB model, the definition of the time unit is identified with the quantity of RND generated. Then, when this quantity coincides with a multiple of the natural number M, chosen previously, a new time unit is attained. This is another of the CORUB model’s original features.
M1䡠 ⫹ FM ! M1䡠
2.5. The parameter M and the temperature
M1䡠 ⫹ FM ! M2䡠
e⫺5
M2䡠 ⫹ FM ! M1䡠
e⫺6
For homogeneous reactions in the gas phase, the frequency of intermolecular collisions in the volume unit is intimately bound to the temperature of the reacting system [20]. Then the CORUB model, for the isothermal chemical act of reaction, can simulate the influence of the temperature in the reaction rate. This can be made possible because the model’s behaviour is related to the maximum frequency of the intermolecular collisions activated by the unit of volume with the aforementioned parameter M. Therefore, in accordance with the Collision Theory [21,22] it can be demonstrated that [11]: 1 R lZ0 ln
8 M T Ea where T is the absolute temperature, R the universal constant of gases, Ea the activation energy, l the own constant model, and Z0 the total frequency of the bimolecular collisions. It can be added that, there is a very close relationship between the parameter M and the stochastic kinetic constant value (k 0 ) in the model. As such, it can be demonstrated that [11]: k 0 k 00 M
9
k 00
signifies the minimum value of the kinetic rate
where
M1䡠 ⫹ M1䡠 ! P
M2䡠 ⫹ M1䡠 ! P
e⫺3 e⫺4
e⫺7
M2䡠 ⫹ R䡠 ! P
e⫺8
M1䡠 ⫹ R䡠 ! P
e⫺9
FM ⫹ R䡠 ! M2䡠 M1䡠 ⫹ FM ! M1䡠
e ⫺ 10 e ⫺ 11
3. Simulation Modelling the photopolymerization of FM with CORUB was carried out taking into account the best proposed group of kinetic constants (Table 1), with the initial conditions employed by Lange et al. [1]. For these reasons the phenomenological description of the stochastic simulation will be given considering the following terms: c
FM0 6:25 mol L⫺1 ; c
AIBN0 1:5 × 10⫺2 mol L⫺1 ; irradiation time, 1020 s; intensity of absorbed light, 1:91 × 10⫺6 Einstein L⫺1 s⫺1 ; initiation quantum yield, 0.4868; and temperature, 0⬚C.
92
J.A. Corzo E. et al. / Computational and Theoretical Polymer Science 11 (2001) 89–94
Table 2 Stochastic concentrations (mol L ⫺1) given by CHEMOD-X for the chemical species involved as a function of time with M 4416727: The error observed for the concentration of the monomer (FM) appears in the last column t(s)
c
R䡠 × 1010
c
M1䡠 × 108 c
M2䡠 × 104 c(FM) Error in c(FM)
60 120 180 240 300 360 420 480 540 600 660 720 780 840 900 960 1020
1.661 1.661 1.661 1.661 1.661 3.322 3.322 1.661 3.322 3.322 1.661 1.661 1.661 1.661 1.661 4.983 1.661
1.030 2.110 2.757 3.306 4.718 4.784 5.365 5.133 6.495 6.561 7.409 7.791 8.538 8.339 8.488 9.236 8.688
1.536 2.962 4.275 5.487 6.603 7.633 8.574 9.434 10.22 10.94 11.59 12.17 12.70 13.17 13.59 13.97 14.30
6.206 6.162 6.118 6.074 6.030 5.986 5.942 5.898 5.854 5.810 5.766 5.722 5.678 5.634 5.590 5.546 5.502
0.007 0.012 0.017 0.022 0.025 0.028 0.030 0.032 0.033 0.034 0.035 0.035 0.036 0.036 0.035 0.035 0.035
3.1. Determining an appropriate M value When using CORUB to reproduce a particular reaction mechanism, it is very important to know which value of M will make it possible. As such, the first step of the stochastic simulation consists of determining the minimum rate of polymerization that CORUB can model with CHEMODX. With this objective, using M 20000; an auxiliary simulation is carried out. The value obtained, 1:6612 × 10⫺10 mol L⫺1 s⫺1 ; allows us to calculate an M that could reproduce the polymerization rate observed
7:3364 × 10⫺4 mol L⫺1 s⫺1 [1]. Then, making M 4416727 and considering the same initial conditions again, CHEMODTable 3 Calculated concentrations (mol l ⫺1) by simulation at 0⬚C of all the involved species in the proposed reaction mechanism as polymerization time (t in seconds) elapses [1] t (s)
c(R 䡠 ) × 10 10
c(M 䡠1) × 10 8
c(M 䡠2) × 10 4
c(FM)
60 120 180 240 300 360 420 480 540 600 660 720 780 840 900 960 1020
1.652 1.654 1.658 1.662 1.668 1.675 1.684 1.693 1.703 1.714 1.726 1.739 1.752 1.765 1.780 1.794 1.809
0.732 1.430 2.104 2.744 3.344 3.896 4.400 4.852 5.255 5.609 5.919 6.187 6.419 6.617 6.786 6.930 7.051
1.112 2.205 3.261 4.265 5.205 6.072 6.862 7.573 8.205 8.763 9.250 9.672 10.04 10.35 10.61 10.84 11.03
6.247 6.239 6.226 6.208 6.185 6.158 6.127 6.093 6.056 6.017 5.975 5.932 5.887 5.842 5.795 5.747 5.699
Table 4 Results obtained from the adequacy analysis of the model using the Kolmogorov–Smirnov test (MNS, minimum level of significance in order to refuse H0) Chemical species MNS
FM 0.4540
M1䡠 0.1120
M2䡠 0.1120
R䡠 0.0463
X is executed once more. The outputs obtained for this stochastic simulation are shown in Table 2, where the error related to reference data (Table 3) [1], for each one of the monomer concentration values, are reported in the last column. 4. Results and discussion 4.1. Confirmation of the outputs obtained by the direct task In Table 3, the concentrations given by Lange [1] for all the species that participate in the mechanism are assembled. As can be seen, in Table 2, the proposed model reproduces the values for the concentrations of all the species involved in the mechanism of the studied reaction. Furthermore, it is observed that CHEMOD-X gives the higher values of concentration for the radical species M2: and the smaller ones for the radical R䡠 : It is very important to realize that the stochastic concentration of FM has a maximum error of 3.5%. This fact shows the good level of accuracy of the simulation. In order to confirm the previous paragraph, for all the chemical species involved, the Kolmogorov–Smirnov test with a significance level / 0:01 was employed. From the results shown in Table 4, the adequacy of the reference data’s description can be seen because in all the cases, the fixed / is less than the corresponding value of the minimum significance level, consequently significant differences between the behaviour of both series of concentration values do not exist; that is to say, they are equivalent. This statistical analysis was carried out with the STATGRAPHICS 6.1 software and its results show the suitability of the fit of the CORUB model. 4.2. The temperature simulation Table 5 reproduces the average values observed for the Table 5 Polymerization rate (Rp) observed for FM and calculated value of M. c
FM0 6:25 mol L⫺1 ; c
AIBN0 1:5 × 10⫺2 mol L⫺1 time of irradiation: 1020 s T (⬚C)
Rp × 10⫺4 (mol L ⫺1 s ⫺1)
M
0 10 20 30 40
7.46 9.33 9.74 11.50 17.86
4490865 5616591 5863408 6922915 10751588
J.A. Corzo E. et al. / Computational and Theoretical Polymer Science 11 (2001) 89–94
93
Fig. 2. Relation between the parameter M and the temperature according to CORUB.
photopolymerization rate as a function of temperature [1]. In this table, a third column is used to show the calculated values of M. This calculation was done taking into account the minimum polymerization rate that CHEMOD-X performs with its current coding. The values that appear in Table 5 allow us to check Eq. (8), now applied to the photopolymerization reaction. This dependence between the reciprocal of the absolute temperature and the logarithm of the inverse of M are appreciated in Fig. 2. The correlation analysis effected in Fig. 2 allows us to formulate the following equation: 1 1 0:01177 ⫹ 5:32482 × 10⫺4 ln
10 T M for Eq. (10), the correlation coefficient has the value 0.9383. In order to verify the statistical level of significance of r, it was checked that it differs significantly from zero with / Table 6 Counting result of the occurrence for each step in the stochastic simulation of the studied photopolymerization. c
FM0 6:25 mol L⫺1 ; c
AIBN0 1:5 × 10⫺2 mol L⫺1 ; time of irradiation, 1020 s; and M 4416727 Step
Occurrences
% of occurrences
e-2 e-3 e-4 e-5 e-6 e-7 e-8 e-9 e-10 e-11
9125362 1690247391 13239 1399521238 1393909744 3486825 10433 42 6496135 2251092
0.2026 37.5189 2.94 × 10 ⫺4 31.0655 30.9410 0.0774 2.32 × 10 ⫺4 9.32 × 10 ⫺7 0.1442 0.0500
0:05; where the fitness is significant. Then, it can calculate the value of the apparent activation energy for the polymerization by CORUB, which is equal to 3.7 kcal mol ⫺1, very close to the observed value (3.4 kcal mol ⫺1) [1]. Eq. (10) can also simulate the temperature value with an error less than 2.5% [11].
4.3. Other features of the simulation When CORUB is employed we can gain a much better insight into the chemical process of the reaction. It can count how many times each reaction step of the mechanism is selected and which of these selected cases correspond with a successful chemical reaction. In this way, CORUB can perform a sensitivity analysis with an alternative technique for the photopolymerization already taken into account, and we can also estimate the time that photoinitiation needs to start. The former, to the author’s knowledge, has never been done previously.
4.3.1. Estimation of photoinitiation time The CHEMOD-X program can determine the quantity of random numbers that were generated before the first occurrence of the photoinitiation step. This allows us to calculate the delay before the process begins. Up to now, these kinetic simulation methods have not been employed to carry out such a task. Then, by considering the same initial conditions, the computing program was executed ten times to obtain the average delay in photoinitiation. The result obtained was 8:93 × 10⫺5 s: Therefore, this supports the fact that it is practically instantaneous.
94
J.A. Corzo E. et al. / Computational and Theoretical Polymer Science 11 (2001) 89–94
4.3.2. An alternative approach to the sensitivity analysis The CHEMOD-X program counts the number of times that each one of the reaction steps are executed. The outputs for stochastic photopolymerization modelling are shown in Table 6. In this table, the percentage of occurrence of each step of the studied mechanism is presented. These results agree perfectly with the analysis of sensitivity [1], and the steps e-4, e-8, e-9 and e-11 are shown to be less sensitive since they have the lowest percentages of occurrence. It can be added, that there are three fundamental steps: the acrylic macroradical propagation, the intermolecular degradative transfer and the re-initiation process (e-3, e-5 and e-6, respectively). This new approach to the sensitivity analysis is a very simple demonstration of CORUB’s possibilities.
5. Remarks As it can be appreciated, to simulate a mechanism of reaction with CORUB, admitting a simple expound in general, it is only necessary to combine the stochastic simulation of each individual step. The parameter M makes it possible to model both time and temperature with notable success. Consequently, the parameter previously mentioned allows the estimation of both the starting time of photoinitiation and the energy of activation associated with the global reaction, with sufficient accuracy. Besides this, a new alternative approach to the sensitivity analysis has been presented in a simpler form.
References [1] Jurgen L, Rieumont J, Davidenko N, Sastre R. Polymer 1998;39(12):2537. [2] Metropolis N, Ulam SM. J Am Statist Assoc 1949;44(247):335. [3] Householder AS, Forsythe GE, Germond HH. Monte Carlo method, Applied Mathematics Series, vol. 12. National Bureau of Standards, 1951. [4] Manork JJ. In: Mattson JS, Mark HB, McDonald HC, editors. The application of Monte Carlo method to chemical kinetics, New York: Marcel Dekker, 1973. [5] Beech G. Fortran IV in chemistry: an introduction to computer assisted methods, vol. 98. London: Wiley, 1975. [6] Niederreiter H. Bull Am Math Soc 1978;84(6):957. [7] Bruns W, Motoc I, O’Driscoll KF. Monte Carlo applications in polymer science. Berlin: Springer, 1981. [8] Sobol IM. Me´todo de Monte Carlo 2nd. Moscow: Mir, 1983. [9] Corzo JA, Alvarez R. Rev Cub Quı´m 1998;10:1–2. [10] Corzo JA, Alvarez R, Sua´rez R. Tecnologı´a Quı´mica 1998;18:3. [11] Corzo JA. PhD dissertation, Havana University, Havana, Cuba, 1998. [12] Kibby ML. Nature 1969;222:298. [13] Nakanishi T. J Phys Soc Jpn 1972;32:1313. [14] Dixon DA, Shafer RM. J Chem Educ 1973;50:648. [15] Moebs WDC. Math Biosci 1974;22:113. [16] Bunker DL, Garret B, Kleindienst T, Long GS. III Combust Flame 1974;23:373. [17] Weber J, Celardin F. Chimia 1976;30:236. [18] Lu J, Zhang H, Yang Y. Makromol Chem Theory Simul 1993;2:747– 60. [19] He J, Zhang H, Yang Y. Macromol Chem Theory Simul 1995;4:811– 9. [20] Halliday D, Resnick R. Physics, T I. Havana, Cuba: Cuban Book Institute, 1971 (chap. XXIII). [21] Glasstone S. Treatise of chemical physics. Havana, Cuba: Cuban Book Institute, 1968 (chap. XIII). [22] Guerasimov Y, et al. Course of physical chemistry, T II. Moscow: Mir, 1980 (chap. X).