Letter Cite This: ACS Macro Lett. 2017, 6, 1414−1419
pubs.acs.org/macroletters
Kinetic Monte Carlo Simulation for Quantification of the Gel Point of Polymer Networks Rui Wang,† Tzyy-Shyang Lin,† Jeremiah A. Johnson,‡ and Bradley D. Olsen*,† †
Department of Chemical Engineering and ‡Department of Chemistry, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139, United States S Supporting Information *
ABSTRACT: Accurate prediction of the gel point for real polymer networks is a long-standing challenge in polymer chemistry and physics that is extremely important for applications of gels and elastomers. Here, kinetic Monte Carlo simulation is applied to simultaneously describe network topology and growth kinetics. By accounting for topological defects in the polymer networks, the simulation can quantitatively predict experimental gel point measurements without any fitting parameters. Gel point suppression becomes more severe as the primary loop fraction in the networks increases. A topological homomorphism theory mapping defects onto effective junctions is developed to qualitatively explain the origins of this effect, which accurately captures the gel point suppression in the low loop limit where cooperative effects between topological defects are small.
P
olymer networks are widely used for applications varying from commodity materials like elastomers and superadsorbers1−5 to functional materials used in tissue engineering,6,7 drug delivery,8,9 soft electronics,10,11 and smart actuators.12,13 Gelation, which characterizes the transition from a viscous fluid to an elastic solid due to progressive linking of polymer chains, is one of the most important problems in the study of polymer networks.1−4 Although determining the transition point from sol to network, known as the gel point, is crucial for the application of these materials, quantitative prediction of the gel point for real polymer networks remains an outstanding challenge.14−16 Most of our fundamental knowledge about gelation is based on the classical Flory−Stockmayer (F−S) theory that assumes a loop-free treelike structure;17,18 however, real polymer networks inevitably possess cyclic defects: loops formed by intrinsic intramolecular reactions.19,20 These wasted internal connections suppress gelation, and their absence in network models is thought to be the essential difference between the Flory−Stockmayer prediction and experimental measurements.1,2,21 Since it is difficult to characterize the structure of cross-linked networks,22 a reliable theory or simulation that enables calculation of the gel point, verified by experiments, would play a central role in understanding the formation−structure−property relationships of polymer networks.23−26 © XXXX American Chemical Society
Despite many efforts to modify the F−S theory or to use alternative ways to compute the same results,27−34 the kinetics of loop formation and the topology of loop structures in polymer networks have not been properly accounted for in the theoretical prediction of the gel point. Existing gelation theories often rely on presumed network configurations17,18,27,28 and approximate estimations of loop densities.29−34 The topological constraint due to the saturation of the junction functionality, accounted for in the theoretical description of the loop structures,46 is missing in existing gelation theories,29−34 which leads to overestimation of the density of higher-order loops. Furthermore, existing theories usually assume that loops are isolated in the polymer networks. The cooperative effect due to adjacent loops in the formation of defective connections is also ignored,29−34 which results in a discrepancy between theoretical prediction and experimental measurement of the gel point as the loop fraction increases. Due to the lack of proper treatment of loops, existing theories and simulations often involve fitting parameters to reconcile gel point predictions with experimental data.31−34 Therefore, it is desirable to Received: August 8, 2017 Accepted: November 20, 2017
1414
DOI: 10.1021/acsmacrolett.7b00586 ACS Macro Lett. 2017, 6, 1414−1419
Letter
ACS Macro Letters
other hand, if this distance is infinite, intermolecular reaction occurs which results in a change in both the number and size distribution of clusters. The simulation algorithm can be generalized to more complicated polymer networks with disperse chain length, other chain statistics, and junction functionality. Monte Carlo simulations clearly capture gel point suppression due to an increased density of topological defects. The gel point pc in our simulations is the critical conversion for the transition from finite size clusters to an infinite network, which can be determined from the cluster size distribution ns(p) recorded in the simulation. The weight-average cluster size DPw at conversion p is defined as
systematically account for the impact of loops on the growth and percolation of polymer networks to quantify the correlation between cyclic defects and the gel point. Recently, we reported “network disassembly spectrometry” (NDS), which provides the first experimental approach to directly count the number of primary loops in polymer networks.35−38 Motivated by these new experimental data, we modified the kinetic graph theory39,40 based on previous work of Stepto41 and other researchers42,43 to quantify the cyclic topology of polymer networks. Without any fitting parameters, the primary loop fraction predicted by the theory shows good agreement with the NDS experimental measurements for endlinked polymer networks prepared over a wide range of polymer concentrations and chain lengths.39 We previously applied kinetic Monte Carlo simulation to quantify the cyclic topology of polymer networks with different junction functionalities,44 matching direct measurements of primary loops in trifunctional and tetrafunctional polymer networks with no adjustable parameters.35−37 The kinetic graph theory and kinetic Monte Carlo simulations show a nonmonotonic dependence of higher-order loops on the primary loop fraction and an odd−even alternation of the loop densities on the junction functionality, which results from the topological constraint of the network and has not been properly accounted for in the theoretical prediction of the gel point. In this work, kinetic Monte Carlo simulations are applied to simultaneously quantify the topological defects and growth kinetics of polymer networks during the gelation process. Molecular simulation has been applied to study the loop size and its distribution in polymer networks.26,45,46,60 Molecular simulation also provides an important tool to study gelation,47−51 since it does not need to invoke assumptions that have not been validated in existing gelation theories. The kinetic Monte Carlo simulation algorithm previously developed for counting loops35−37,39,44 is modified to track the cluster size distribution during the process of network growth (see Supporting Information Section I for the detailed description of the MC algorithm). The polymer networks considered are formed via end-linking bifunctional polymers (A2) and tetrafunctional junctions (B4). It has been shown by Nishi et al. that the end-linking process is reaction-limited over a wide concentration range,52 where diffusion of the reactants and adjustment of the chain configuration are much faster than the reaction rate. This reaction-limited assumption allows the spatial location of polymer segments and junctions to be neglected, enabling a large increase of the possible system size (1.5 × 104 ∼ 5 × 104 polymers are included in the simulation). The primary loop fraction predicted in the previous theoretical and simulation work 39 based on the reaction-limited assumption shows good agreement with the NDS experimental measurements without any fitting parameters, which supports the validity of this assumption. In the diffusion-limited regime, modifications to the topological simulation53,54 or the use of real space simulations26,45,46,60 would be required. The pure topological perspective allows us to simplify the simulation algorithm by only tracking the topological distance between reactive groups. The structures of loops and the size of clusters are determined by the topological distances between different molecules, which are updated at each reaction step. Two reactive groups belong to the same cluster if their topological distance is finite. If the topological distance between the selected pair of reactive groups is finite, intramolecular reaction occurs, and the loop densities will be updated. On the
DPw =
∑ s2ns(p)/∑ sns(p) s
s
(1)
where the sum runs over all cluster sizes s. As shown in Figure 1a, DPw increases sharply at a certain p. The cluster growth is
Figure 1. (a) Plot of weight-average cluster size DPw versus conversion p at 6 different levels of primary loop fraction at full conversion of the network (p = 1.0). DPw is averaged over 20 independent simulations. (b) Representative reduced weight-average cluster size DPrw versus p for different primary loop fractions. The maximum of DPrw is taken as the gel point for a given system size.
investigated for networks of varying loop density, which is defined as the number of polymer chains contained in primary loops divided by the total number of chains at full conversion of the network. Due to the one-to-one correspondence between the primary loop fraction and the dimensionless variable cR3 (c is polymer concentration and R is the root-mean-square end-toend distance of polymer chains),39 the loop density in the network can be controlled by adjusting cR3 in the simulation. Loop-free networks are prepared by forbidding the looping 1415
DOI: 10.1021/acsmacrolett.7b00586 ACS Macro Lett. 2017, 6, 1414−1419
Letter
ACS Macro Letters reaction during the simulation. Figure 1a shows that the sharp increase of DPw is delayed to later conversion as the primary loop fraction in the network increases. Since the divergence of DPw at the gel point cannot be realized in a finite size simulation due to the restriction on size of the largest cluster, the reduced weight-average cluster size DPrw is employed to aid accurate determination of the gel point.39,41,43 The reduced weight-average cluster size is defined as DP rw =
∑ s2ns(p)/∑ sns(p) s*
s*
(2)
with the sum running over all clusters excluding the largest one. The maximum of DPrw as shown in Figure 1b is taken as the gel point. Intuitively, the sudden drop in DPrw occurs as the second largest cluster chemically bonds with the largest one, corresponding to the point at which percolation is achieved. Extrapolation of the maximum in DPrw for finite systems enables the gel point to be estimated for a thermodynamically large system. The scaling relation pc − pmax (L) ∼ L−σ between the p at which a maximum of DPrw in the finite size system is observed (denoted by pmax(L)) and the system size (reflected by the total number of polymers L) is used to extrapolate to an infinitely large system.43 The parameter σ is the critical exponent that characterizes the divergence of the largest cluster near the gel point.14 However, the value pc is insensitive to the value of σ over a wide range (0.3 ≤ σ ≤ 0.7, see Figure S1 in Supporting Information Section II). Therefore, σ = 1/2, the value obtained in the Flory−Stockmayer theory,1,14 is chosen to avoid challenges with nonlinear fitting and to reduce the number of parameters in the fit. The gel points for infinite systems are quantified from the intercepts of these linear fits. Figure 2 shows good linear relationships between pmax(L) (averaged over 20 independent simulations) and L−1/2 for a series of primary loop fractions in the network. Comparisons between the simulations and experimental results of gel point suppression show good agreement with no adjustable parameters. Figure 3 shows pc as a function of the dimensionless variable cR3, which characterizes the ratio between the intramolecular distance and intermolecular distance. The pervaded volume R3 = N3/2b3 can be expressed
Figure 3. Gel point pc versus the dimensionless variable cR3 (c is polymer concentration and R is the root-mean-square end-to-end distance of polymer chains) predicted by the simulation in comparison with four sets of experimental data by Wile,55,56 by Gordon and Scantlebury,57 and by Cail and Stepto24 (denoted by Exp1 and Exp2 for networks prepared by polyurethanes with 29 and 33 Kuhn monomers, respectively). The inset shows the dependence of pc on the primary loop fraction. The error bars on the pc represent 95% confidence intervals obtained by the linear fitting (as shown in Figure 2) over 20 independent simulations.
in terms of the number of Kuhn monomers (N) and the Kuhn length (b) that are either experimentally measurable or reported in the literature. As shown in Figure 3, the gel point predicted by the simulation is in good agreement with four independent sets of experimental results. Gelation experiments by Wile55,56 and by Gordon and Scantlebury,57 respectively, on the polycondensation of adipic acid and pentaerythritol are compared to the simulation using the intrinsic parameter of the Kuhn monomer (b = 0.247 nm) reported in the literature.56,58 Experimental data for end-linking of polyurethanes (PUs) of two different lengths (N = 29 and 33) by Cail and Stepto24 are also compared to the simulation using the Kuhn length b = 0.308 nm reported for PU.33 Without any fitting parameters, the simulation can quantitatively describe gel point suppression in networks synthesized by different types of polymers, which highlights the essential role of loops in gel point suppression. Other effects, such as packing,26 entanglement,59 and compositional fluctuations,42 may affect the looping probability and the resulting topology and percolation of the network. However, the results here indicate that the gel point for reaction-limited systems can be well-quantified using a simple model based on a purely topological perspective. The quantitative agreement between the simulation and experiments also suggests that packing effects in the vicinity of the gel point are not essential for quantifying the gel point, perhaps because only a small fraction of reactive groups are consumed in the Ginzburg zone due to its narrow width compared to the mean-field region. The deviation between the simulation and experiment increases as cR3 decreases, particularly in the regime cR3 < 1. This observation is consistent with the increase of the width of the Ginzburg zone, which scales as N−1/3 (or R−1/2).14 Gel point suppression is directly correlated with the density of primary loops in the network, as shown in the inset of Figure 3. pc increases linearly with primary loop fraction when the primary loop fraction is small (
