&KHPLFDO3K\VLFV/HWWHUV  ²Contents lists available at ScienceDirect

Chemical Physics Lettersjournal homepage: www.elsevier.com/locate/cplett

Research paper

A multiscale model for the synthesis of thermosetting resins: From theaddition reaction to cross-linked network formation

7⁎

Jing Lia, , Sakamoto Jumpeia, Hiroki Waizumib, Yutaka Oyaa, Yue Huangc, Naoki Kishimotob,
Tomonaga Okabea,caDepartment of Aerospace Engineering, 6-6, Aramaki Aza-Aoba, Tohoku University, Sendai 980-8579, Japan
Graduate School of Science and Faculty of Science, 6-3, Aramaki Aza-Aoba, Tohoku University, Sendai 980-8578, Japanc
Department of Materials Science and Engineering, University of Washington, BOX 352120, Seattle, WA 98195, Unites Statesb

H I GH L IG H T S

model simulated the inﬂuence of conditions on thermosetting resins.
• AWemulti-scale
the eﬀect on the structure and physical properties.
• We simulatedon the synthesis of the phenol formaldehyde resin.
• The focusedreaction model included the addition and condensation reactions.
• The comprehensive• eﬀectiveness of the reported model was validated.

A R T I C LE I N FO

A B S T R A C T

Keywords:
Thermosetting resin
Multiscale model
Polymerization
Cross-linked network
Molecular dynamics

A multi-scale model that can simulate the inﬂuence of the synthetic conditions (e.g., molar ratio) on the structure
and physical properties of thermosetting resins (e.g., phenol resin) through the introduction of a comprehensive
reaction model that includes the addition reaction was developed and reported herein. The eﬀectiveness of themodel was validated by verifying the inﬂuence of the primary synthetic index, i.e., the molar ratio, on thephysical and thermal properties (i.e., density and glass transition temperature) of Resol resins. This model can
aid in developing more eﬃcient and suitable synthetic conditions to obtain products exhibiting the desired
material properties.

1. Introductionin the system, or where the reaction is considered to be a direct condensation reaction, thereby rendering examination of the inﬂuence ofthe synthetic conditions rather challenging. Thus, the eﬀect of themolar ratio and the monomers employed on the prepolymer size cannot
be predicted. Due to the inﬂuence of these factors on estimation of thephysical properties (e.g., density and glass transition temperature [4])of a polymer, the importance of the addition reaction to the cross-linked
structure appears to have been underestimated.
A model which can discuss the inﬂuence of the synthetic conditions
is therefore required. More speciﬁcally, this model should be able toclarify the relationship between the synthetic conditions, the resulting
polymeric structure, and the physical properties. We also expect that
this could reduce the cost of the synthetic process, while allowing thedesign of resins exhibiting desirable properties.
In recent years, very few studies have focused on the synthetic

As a class of synthetic thermosetting polymers, phenol resin exhibitsa range of properties determined by the chemistry of its constituent
monomers and how these monomers combine with each other in athree-dimensional network structure upon heating [1]. Compared tothe continuum model, atomic-level simulations enable us to examine
this cross-linking process and obtain information relating to the microscopic response from the fundamental chemical structure and the
dynamics of each atom, which ultimately determine the macromechanical properties of the polymer. Indeed, several groups [2,3]have employed atomic-level molecular dynamics (MD) simulations to
discuss the physical properties of phenol resin. However, all previous
models have ignored or simpliﬁed the chemical reaction process, where
the linear polymer chain has been directly used as the initial molecule⁎

Corresponding author.
E-mail address: lijing@plum.mech.tohoku.ac.jp (J. Li).

https://doi.org/10.1016/j.cplett.2019.02.012
Received 18 September 2018; Received in revised form 3 February 2019; Accepted 4 February 2019
$YDLODEOHRQOLQH)HEUXDU\
(OVHYLHU%9$OOULJKWVUHVHUYHG

&KHPLFDO3K\VLFV/HWWHUV  ²

J. Li, et al.

intrinsic reaction coordinate calculations were carried out to conﬁrm
the reaction pathway. In addition, the Global Reaction Mapping
Method (GRRM) was employed to ﬁnd transition states by uphillwalking from equilibrium structures by anharmonic downward distortion following algorism [18]. The activation energies were then calculated by taking the energy diﬀerence between the transition state andthe reactants. Such activation energy simulation work has also been
carried out previously [19].

model due to the diﬃculty in representing the chemical reaction that
results in formation of the three-dimensional cross-linked structure [5].
In 2001, the reactive force ﬁeld ReaxFF [6] was developed for MD simulations to enable chemical reactions involving bond formation andbreakage to take place. However, such simulations are particularly time
consuming. The large-scale model is only employed to simulate systems
containing ∼10,000 atoms [7]. This contrasts to the nonreactive force
ﬁeld, which can be applied to billion atoms systems [8]. Another approach is the use of Monte Carlo (MC) simulations, which have been
employed since the 1950s for building polymer networks. One of themain advantages of MC over MD simulations is that the MC simulations
are rapid, even for fully atomistic modeling [9]. Notably, MC models
based on kinetic theory are becoming popular [10,11], as such models
built for polymer networks can be used to represent the various distributions of the ﬁnal product following the polymerization reaction.
However, the chemical reaction hypothesis based on sphere capture islimited, as description of the non-bonded interactions between atoms ischallenging, and these interactions play a crucial role in describing thephysical cross-linking reaction and the resulting conﬁguration. In contrast, the use of MD models with a force ﬁeld based on van der Waals
and electrostatic interactions can solve this problem.
Thus, to address the above issues, we herein report the developmentof a multi-scale model based on quantum chemistry calculations, MC
simulations, and MD simulations to reproduce the polymerization
process of phenol resin. Indeed, the MC/MD model is identical to that
previously presented by our group for the cross-linked epoxy resin[12,13]. However, in this case, we adopt a complex reaction model
based on quantum chemistry combined with the MC/MD model. Morespeciﬁcally, we consider the two reaction steps involved in the polymerization process [14], namely the addition of formaldehyde to theortho- and para-phenol positions, and the subsequent condensation reaction. Furthermore, a range of monomer molar ratios are examined
due to the importance of this variable in determining the network crosslinking structure [15].

2.3. The polymerization process by MC and MD simulations
During the PF polymerization process, the diﬀerent monomers reacts with one another to form the prepolymers, which then react to givethe ﬁnal polymer [15]. To reproduce this process, quantum chemistry
was initially used to calculate the activation energies (the reactionenergy barriers), which was converted to give the reaction probability.
Subsequently, combination of the MC algorithm with the reaction
probability was used to reproduce this random event. Finally, the force
ﬁeld theory in the MD simulation was applied to ﬁnd the minimum
total potential energy state of the system.
In the MC/MD simulation employed to reproduce the cross-linking
reaction of the phenol resin system, the dynamics (relaxation) calculation and the recombination of covalent bonds that are associated withthe reaction are repeatedly calculated as follows. Firstly, all pairs ofreaction sites are determined between a phenol and a formaldehyde
molecule within 4 Å of one another, which is achieved using the radial
distribution functions. Secondly, MC calculations are used to determine
whether a reaction occurs for each pair or not. In this stage, the random
number (P), which ranges from 0 to 1, compares to the Arrhenius typereaction probability (K):E
K = A exp ⎛− a ⎞,⎝ RT ⎠(1)

where R is the gas constant, T is the local temperature obtained fromthe reacting molecules, Ea is the activation energy obtained from the
GRRM calculations, and A is an acceleration factor (1018) employed tobridge the time gap between an actual reaction (hour scale) and the MD
simulation (picosecond scale). In the subsequent step, bond breakage
and bond creation are considered to occur simultaneously if the random
number is greater than the reaction probability (i.e., P > K). Finally,geometrical optimization and relaxation calculations for the system are
performed over 2 ps. By repeating the above steps under iso-thermal(298 K) and iso-pressure (1 atm) conditions, we can obtain phenol resin
systems with their corresponding curing conversions (α), which aredeﬁned by:

2. Materials and methods
2.1. Formation of a Resol resin
Although the highly cross-linked phenol formaldehyde (PF) resin iscommercially valuable, once the highly intractable network has been
formed, it cannot be molded due to the irreversible changes involved in
its formation. The preparation process is therefore carried out in twostages, namely preparation of the prepolymer, and subsequent growth
polymerization. More speciﬁcally, the desired prepolymers are ﬁrstly
obtained in an alkaline medium through the reaction of formaldehydewith the ortho- and para-positions of the phenol substrate, as shown in
Fig. 1(a). Subsequently, during growth polymerization, the cross-linked
network structure is formed through condensation of the prepolymers
and the elimination of water, as indicated in Fig. 1(b).

α=freactftotal,(2)where ftotal is the number of the phenol units in molecule before reaction, and freact is that of the reacted phenol molecular during the crosslinking process. It should be noted that the water molecule formed
during the condensation reaction (see Fig. 2) is omitted for simplicity.
Following the cross-linking simulation, we evaluated the thermal
properties (density and glass-transition temperature) for the systemwith each curing conversion. The density (ρ) is deﬁned by the total
weight of the particle to the cell volume, which is obtained by relaxation dynamics under the NPT (P = 1 atm and T = 298 K) ensemble. The
glass-transition temperature (Tg) is obtained by ﬁnding the intersection
temperature of two straight lines on volume vs temperature plots inwhich each data point is a result of the relaxation dynamics under the
NPT ensemble where T is increased from 220 to 580 K while maintaining P at 1 atm.
As mentioned above, our MD simulation procedure is divided into
two stages, namely the curing dynamics for obtaining the polymerization structure, and the relaxation dynamics to determine the thermal
properties. In the case of curing dynamics, the COMPASS II force ﬁeld

2.2. Quantum chemistry simulations for the activation energy
We initially studied the reactions at the ortho- and para- positions ofphenol using the transition state theory. The activation energy of thesereactions was elucidated using the Molecular Orbital method, with the
HF/3-21G level of theory being employed for the transition state. The
solvation environment was mimicked by implicit water molecules (i.e.,the solvent) using the polarizable continuum model (PCM) [16] in aquantum chemistry package program, Gaussian 09 [17]. Detailed information regarding the method employed is summarized in Table 1.
The geometric structures for formaldehyde (H2CO), phenoxide ion(C6H5O−), and the reaction products (C6H4O(CH2OH)−) were optimized, and frequency calculations were used to conﬁrm all structures
were at their true minimums. The kinetics of the reaction were thenrevealed by optimizing the transition state of the reaction, where

&KHPLFDO3K\VLFV/HWWHUV  ²

J. Li, et al.

Fig. 1. (a) Chemical structures of the phenol and formaldehyde substrates/monomers employed for preparation of the phenol resin. The reactive positions of phenol
and formaldehyde are shaded in red. (b) During the polymerization process, prepolymer formation takes place through addition reactions between the monomers.
The various prepolymers then react via a condensation reaction to eliminate water and form a cross-linked network structure. (For interpretation of the references tocolor in this ﬁgure legend, the reader is referred to the web version of this article.)[21] was applied and the Materials Studio 8.0 was utilized. In contrast,
to calculate the thermal properties, the TEAM force ﬁeld was employed
along with the LAMMPS software package. To validate the eﬀectivenessof this quantum mechanics/MC/MD model, the inﬂuence of the primary synthetic index, i.e., the molar ratio, on the mechanical properties(i.e., the Resol glass transition temperature (Tg)) was veriﬁed. Morespeciﬁcally, formaldehyde/phenol molar ratios of 1.2, 1.4, 1.5, 1.6, 2,
2.5, and 3 were examined, which are the typical molar ratios employedin the synthesis of Resol resin [20], and the atom numbers of phenol
was kept as constant, 2400. Fig. 2 shows snapshots of the three

Table 1
Quantum chemistry computational conditions employed based on
the Molecular Orbital method for determination of the activationenergy.
Calculation
Level of Theory
Spin
Environment

GRRM11 with Gaussian09
HF/3-21G
SingletAlkaline, Charge is −1

Fig. 2. Example snapshots of each equilibrium structure in the initial conditions with molar ratios of 1.2, 2, and 3, and atom numbers ranging from 3560 to 5000. The
bulk size is on the nanometer scale.


&KHPLFDO3K\VLFV/HWWHUV  ²

J. Li, et al.

Fig. 3. The activation energies and heats of formation for the addition reactions as determined by GRRM calculations. (a) In the reaction at the ortho-position, theactivation energy is 55.05 kcal/mol and the heat of formation is 9.55 kcal/mol. (b) In the reaction at the para-position, the activation energy is 50.20 kcal/mol andthe heat of formation is 13.22 kcal/mol.

observed using this multi-scale model, and the variation in the atomic
system was reported, as shown in Fig. 4. The monomers were randomly
distributed in the initial atomic system (Fig. 4(a)), where a formaldehyde:phenol molar ratio of 2 was employed along with a bulk sizeof a 4.018 nm3 cubic cell with 4200 atoms. Thus, the addition reactionsin the para- and ortho-positions tended to occur following growth of thepolymerization degree, and in this state (see Fig. 4(b)), small molecules
were the predominant species. Subsequently (Fig. 4(c)), the prepolymers reacted with one another to give linear and branched polymers. When the polymerization reaction was allowed to proceed until
close to completion, large-scale molecular and network structures wereformed. These results indicate that the polymerization process involving both addition and condensation reactions is successfully reproduced in our quantum mechanics/MC/MD multiscale model.
The room temperature (298 K) density results simulated under different synthetic conditions (i.e., formaldehyde to phenol molar ratios
ranging from 1.2 to 2.5) are shown in Fig. 5. Importantly, our results
agree with commercial values ranging from 1.15 to 1.3 g/cm3 [23] andexperimentally measured values reported previously [24]. In contrast,
previously reported simulated density results that do not consider the
addition reaction were < 1.15 [2,3], while the simulated density results
where the addition reaction was included are closer to the typical PF
density data, i.e., 1.24 g/cm3 [25].
Diﬀerent Resol resins were then characterized based on the freephenol content following preparation, and the simulation result isshown in Fig. 6 along with a comparison with previously reported

equilibrium structures in the initial conditions with diﬀerent molar
ratios and atom numbers.
3. Results and discussions
The activation energy and heat of formation of the addition reactions taking place at the phenolic ortho- and para- positions under alkaline conditions were predicted using the GRRM tool, as shown in
Fig. 3(a) and (b). More speciﬁcally, the activation energy for the orthoposition is 55.05 kcal/mol, while that of the para- position is50.20 kcal/mol, indicating that the probability of a reaction occurringat these two positions is similar, because the activation energy can beinterpreted as the magnitude of the energy barrier to the reaction. The
activation energies of the other types of prepolymers, such as o,o-, o,p-,and o,o,p-, have been reported in Ref. [22], and the simulation results
suggested that the reaction probability at the ortho-position was signiﬁcantly higher than that at the para-position based on the obtained
activation energies. Indeed, almost twice the amount of energy is required for reaction at the para-position under acidic conditions. These
results can therefore account for the fact that diﬀerent resin structures
are formed in the reaction between phenol with formaldehyde underdiﬀerent pH conditions [15]. More speciﬁcally, Novolac, which has alinear polymeric chain structure, is prepared under acidic conditions,while Resol, which has a network polymer structure, is synthesizedunder alkaline conditions.
For the synthesis of phenol resin, the polymerization processes were

Fig. 4. Constructed molecular structure of the 4.018 nm3 cubic cell with 4200 atoms in various reaction conversion stages (formaldehyde to phenol molar ratio = 2):(a) Initial condition prior to reaction; (b) Addition reaction with 37% curing conversion; (c) Formation of linear and branched polymer chains following condensation
with 63% curing conversion; (d) and (e) Polymerization continued until the reaction approaches completion to give large-scale molecular and network structures.
&KHPLFDO3K\VLFV/HWWHUV  ²

J. Li, et al.

420
400

1.30
1.20

T [K]

Typical value [25]Com. value range [23]1.85E+04360
1.35E+04
8.50E+03

320
300

1.50

2.00
F/P Molar ratio

2.50

22.0

3.00

5.5
0.0
1.4

1.8
F/P Molar ratio

2.2

2.0
2.5
P/F Molar ratio

3.0

3.5

Tg∞ = 434 K, λ = 0.41, as described in the literature [26]. Indeed, ourestimation of Tg agrees well with the results of the Perez model and withtypical experimental results at 423 K [23].
It should also be noted that molecular size can inﬂuence the value of
Tg. Thus, the simulation results for Tg and for the largest molecular size
(L.M.S) under diﬀerent molar ratios were compared, as presented inFig. 8. The obtained results suggest that the L.M.S exhibits a positive
correlation with the Tg, where larger molecules increase the value of Tg.
These results conform with the literature explanation [27] that theformation of large molecules can increase the number of strong intermolecular forces. In contrast, the presence of small molecules that candiﬀuse within the cross-linked system can reduce Tg. Such trends havebeen studied previously [4,28].

11.0

1.0

1.5

Fig. 8. Simulation results for the Tg and the largest molecular size (L.M.S) at
diﬀerent molar ratios. These results conform with the literature explanation[27] that the formation of large molecules can increase the number of strong
intermolecular forces. In contrast, the presence of small molecules that candiﬀuse within the cross-linked system can reduce Tg.

This study
Ref.[20]

16.5

3.50E+03
1.0

Fig. 5. The density simulation was performed under diﬀerent synthetic conditions (i.e., formaldehyde to phenol molar ratios of 1.2–2.5). Literature values
for the commercial [23] and simulated results are also shown for comparison.

Free phenol (%)

2.35E+04

380

340

1.10
1.00
1.00

2.6

Fig. 6. Diﬀerent resol resins were characterized based on their free phenol
contents following preparation, and the simulation results are compared withpreviously reported experimental data [20]. The variation trend of the freephenol with the increasing molar ratio agree well with the experimental data.

4. Conclusions
We herein reported the development of an atomic-level model thatcan simulate the inﬂuence of synthetic conditions (i.e., the molar ratio
and pH) on the physical and thermal properties of the phenol resin,
phenol formaldehyde (PF). This model clariﬁed the relationship between the synthetic conditions employed, the structure formation, andthe physical and thermal properties of the resin by considering acomprehensive reaction model that includes both the addition andcondensation reactions. The eﬀectiveness of the reported model wasthen validated by verifying the inﬂuence of the primary synthetic index(i.e., the molar ratio) on the physical and thermal properties of the
Resol resins, including the glass transition temperature, density. In
addition, the computing cost was reduced through the adoption of amulti-scale model, which combined quantum chemistry calculations,Monte Carlo simulations, and molecular dynamics simulations. Thismodel would be expected to contribute to reducing costs in the synthesis of PF resins by the removal of trial and error, and should also aid indetermining suitable synthetic conditions for obtaining the desired
material properties.

experimental data [20]. As shown, the free phenol content tended todecrease with an increasing molar ratio in both cases, and from both theexperimental and simulation results, it was apparent that phenol existed in the system following completion of the polymerization reaction.
These results also suggest that the synthetic model should consider the
addition reaction, where the initial reactants are the monomers.
The Tg values were then predicted for the various curing conversions, and the results are presented in Fig. 7. More speciﬁcally, the freevolume approach and two ﬁtting curves were used to calculate Tg,
where we used the least squares error method to determine the two
ﬁtting curves, and the Perez model was ﬁt using Tg0 = 268 K, and450

Glass Transition, Tg [K]2.85E+04
Tg
L.M.S

Mass [g/mol

Density [g/cm^3]

1.40

363
275

Declaration of interests

188

The authors declare that they have no known competing ﬁnancial
interests or personal relationships that could have appeared to inﬂuence the work reported in this paper.

This study
Ref.[26]

100
0.0

0.2

0.4

0.6

0.8

1.0

Acknowledgments

Conversion

We thank Dr. Yamashita for his fruitful discussions. We also acknowledge the support of the Cross-ministerial Strategic Innovation
Promotion Program, Tohoku University Center for Gender Equality
Promotion (TUMUG) and the vitally important encouragement andFig. 7. Tg was predicted with variation in the curing conversion. The Perez
model was ﬁt using Tg0 = 268 K, and Tg∞ = 434 K, λ = 0.41, as described in theliterature [26]. The estimation of Tg agrees well with the results of the Perez
model and with typical experimental results at 423 K [25].


&KHPLFDO3K\VLFV/HWWHUV  ²

J. Li, et al.

support of the University of Washington-Tohoku University: Academic
Open Space (UW-TU:AOS). N. K and T. O. express gratitude for a research grant from the Institute for Quantum Chemical Exploration(IQCE).

and molecular simulation, Polymer 54 (2013) 4660–4668.
[13] T. Okabe, Y. Oya, K. Tanabe, G. Kikugawa, K. Yoshioka, Molecular dynamics simulation of crosslinked epoxy resins: curing and mechanical properties, Eur. Polym.
J. 80 (2016) 78–88.
[14] P.J. Flory, Principles of Polymer Chemistry, Cornell University Press, New York,
1953.
[15] G. Odian, Principles of Polymerization, John Wiley & Sons, New Jersey, 2004.
[16] B. Mennucci, Polarizable continuum model, Wiley Interdis. Rev.: Comput. Mol. Sci.
2 (2012) 386–404.
[17] M.J. Frisch, et al.,, Gaussian 09, Revision B.01, Gaussian Inc, Wallingford, CT, 2009.
[18] K. Ohno, S. Maeda, A scaled hypersphere search method for the topography ofreaction pathways on the potential energy surface, Chem. Phys. Lett. 384 (2004)
277–282.
[19] T. Tokiwa, N. Kishimoto, K. Ohno, Automated exploration of isomerization anddissociation pathways of ethylene sulﬁde cation by the global reaction route
mapping method, Chem. Phys. Lett. 641 (2015) 97–103.
[20] L.B. Manfredi, O. de la Osa, N. Galego Fernandez, A. Vazquez, Structure–properties
relationship for resols with diﬀerent formaldehyde phenol molar ratio, Polymer 40(1999) 3867–3875.
[21] H. Sun, Z. Jin, C. Yang, R.L.C. Akkermans, S.H. Robertson, N.A. Spenley, S. Miller,
S.M. Todd, COMPASS II: Extended coverage for polymer and drug-like molecule
databases, J. Mol. Model. 22 (2016) 47.
[22] J. Sakamoto, H. Waizumi, J. Li, Y. Oya, N. Kishimoto, T. Okabe, to be submitted for
publication.
[23] https://ciechgroup.com/en/products/zywice-i-tworzywa/zywice/phenolformaldehyde-resins/.
[24] R.A. Haupt, T.J. Sellers, Characterizations of phenol-formaldehyde Resol resins,Ind. Eng. Chem. Res. 33 (1994) 693–697.
[25] C.P.R. Nair, Advances in addition-cure phenolic resins, Prog. Polym. Sci. 29 (2004)
401–498.
[26] J.M. Pérez, F. Rodríguez, M.V. Alonso, M. Oliet, Time–temperature–transformation
cure diagrams of phenol–formaldehyde and lignin–phenol–formaldehyde novolac
resins, J. Appl. Polym. Sci. 119 (2011) 2275–2282.
[27] J.M.G. Cowie, V. Arrighi, Polymers: Chemistry and Physics of Modern Materials,
CRC Press, Florida, 2007.
[28] E.R. Weeks, J.C. Crocker, A.C. Levitt, A. Schoﬁeld, D.A. Weitz, Three-dimensional
direct imaging of structural relaxation near the colloidal glass transition, Science
287 (2000) 627–631.

References[1] K. Rajan, Materials informatics, Mater. Today 9 (2005) 38–45.
[2] J.D. Monk, E.W. Bucholz, T. Boghozian, S. Deshpande, J. Schieber,
C.W. Bauschlicher Jr., J.W. Lawson, Computational and experimental study ofphenolic resins: thermal– mechanical properties and the role of hydrogen bonding,Macromolecules 48 (2015) 7670–7680.
[3] Y. Shudo, A. Izumi, K. Hagita, T. Nakao, M. Shibayama, Large-scale molecular
dynamics simulation of crosslinked phenolic resins using pseudo-reaction model,
Polymer 103 (2016) 261–276.
[4] R.J. Larsen, C.F. Zukoski, Eﬀect of particle size on the glass transition, Phys. Rev. 83(2011) 051504.
[5] C. Li, A. Strachan, Molecular scale simulations on thermoset polymers: a review, J.
Polym. Sci. B Polym. Phys. 53 (2015) 103–122.
[6] A.C.T. van Duin, S. Dasgupta, F. Lorant, W.A. Goddard, ReaxFF: a reactive force
ﬁeld for hydrocarbons, J. Phys. Chem. A 105 (2001) 9396–9409.
[7] Mo Zheng, Xiaoxia Li, Fengguang Nie, Li Guo, Investigation of model scale eﬀects
on coal pyrolysis using ReaxFF MD simulation, Mol. Simulat. 43 (13–16) (2017)1081–1088.
[8] F.F. Abraham, R. Walkup, H. Gao, M. Duchaineau, T.D. De La Rubia, M. Seager,
Simulating materials failure by using up to one billion atoms and the world's fastest
computer: brittle fracture, Proc. Natl. Acad. Sci. 99 (2002) 5777–5782.
[9] P.V.K. Pant, D.N. Theodorou, Variable connectivity method for the atomistic Monte
Carlo simulation of polydisperse polymer melts, Macromolecules 28 (1995)
7224–7234.
[10] J. He, H. Zhang, J. Chen, Y. Yang, Monte Carlo simulation of kinetics and chain
length distributions in living free-radical polymerization, Macromolecules 30(1997) 8010–8018.
[11] K. Binder (Ed.), Monte Carlo and Molecular Dynamics Simulations in Polymer
Science, Oxford University Press, Oxford, 1995.
[12] T. Okabe, T. Takehara, K. Inose, N. Hirano, M. Nishikawa, T. Uehara, Curing reaction of epoxy resin composed of mixed base resin and curing agent: experiments