Molecular Model for HNBR with Tunable Cross-Link Density

We introduce a chemically inspired, all-atom model of hydrogenated nitrile butadiene rubber (HNBR) and assess its performance by computing the mass density and glass-transition temperature as a function of cross-link density in the structure. Our HNBR structures are created by a procedure that mimics the real process used to produce HNBR, that is, saturation of the carbon−carbon double bonds in NBR, either by hydrogenation or by cross-linking. The atomic interactions are described by the all-atom “Optimized Potentials for Liquid Simulations” (OPLS-AA). In this paper, first, we assess the use of OPLS-AA in our models, especially using NBR bulk properties, and second, we evaluate the validity of the proposed model for HNBR by investigating mass density and glass transition as a function of the tunable cross-link density. Experimental densities are reproduced within 3% for both elastomers, and qualitatively correct trends in the glass-transition temperature as a function of monomer composition and cross-link density are obtained. ■ INTRODUCTION Since its first synthesis in 1930, nitrile butadiene rubber (NBR) has become one of the most widely used elastomers.2−4 Its areas of application are vast and include the aeronautic, automotive, and oil industries,2,4−6 in which its resistance to fuels makes it a suitable material for seals and grommets. NBR is a statistical copolymer consisting of two monomers, acrylonitrile (ACN) and butadiene. Depending on the ratio of these two monomers in the final structure, as well as on the various additives that can be introduced during synthesis, the properties of NBR, such as fuel resistance, tensile strength, heat aging, resilience, and low-temperature flexibility, may be tuned to the particular application. In applications that require elastomer components to operate in relatively harsh environments, for example, in oil wells, where temperatures and pressures can be as high as 500 K and 200 MPa, respectively, hydrogenated NBR (HNBR) is often the preferred material. The hydrogenation process involves using catalysts to react NBR with hydrogen, thereby saturating the carbon−carbon double bonds present in the NBR backbone. The resulting HNBR has enhanced properties compared to those of NBR: resistance to heat, blistering, and amine corrosion are improved, as is elongation to fracture, whereas other key properties that make NBR such an important elastomer in industrial application, such as oil resistance and low-temperature flexibility, are preserved. Failure of critical elastomeric components can have very serious safety and economic consequences. Performing controlled experiments under service conditions to assess performance and elucidate detailed failure mechanisms is often difficult, if not impossible. Conversely, computer simulations enable materials to be studied by varying the parameters of the model in a controlled and systematic manner. Often, these parameters are inaccessible experimentally and therefore simulations can provide information and insight that are complementary to those obtainable from experiments. However, despite the widespread use of NBR and HNBR, there are very few examples of atomistic simulations of these materials: Kucukpinar and Doruker and Song et al. used the all-atom COMPASS force field to compute various properties of NBR chains up to a length of 164 monomers, including the solubilities and diffusivities of small gas molecules. In this article, we present a new chemically inspired all-atom description of HNBR, that is, one that mimics the synthesis of HNBR via saturation of the carbon−carbon double bonds present in the backbone of NBR. An initial NBR structure is first generated by growing an NBR chain one bond at a time in a periodic simulation cell using a Monte Carlo approach. It is then equilibrated with molecular dynamics (MD), using the all-atom “Optimized Potentials for Liquid Simulations” (OPLSAA) force field to describe the interatomic interactions. The carbon−carbon double bonds in the NBR structure are then saturated by either cross-linking or hydrogenation, resulting in the generation of HNBR. Whether a particular double bond is cross-linked or hydrogenated is determined by two factors: (i) a user-defined parameter that sets the target density of the crosslinks in the structure and (ii) whether or not the local atomic environment of the bond permits a cross-link or, in other words, whether or not there is another carbon−carbon double bond available within a user-defined radius, with which to crossReceived: August 3, 2016 Revised: October 13, 2016 Published: November 14, 2016 Article


■ INTRODUCTION
Since its first synthesis in 1930, 1 nitrile butadiene rubber (NBR) has become one of the most widely used elastomers. 2−4 Its areas of application are vast and include the aeronautic, automotive, and oil industries, 2,4−6 in which its resistance to fuels 3,7 makes it a suitable material for seals and grommets.
NBR is a statistical copolymer consisting of two monomers, acrylonitrile (ACN) and butadiene. Depending on the ratio of these two monomers in the final structure, as well as on the various additives that can be introduced during synthesis, the properties of NBR, such as fuel resistance, tensile strength, heat aging, resilience, and low-temperature flexibility, may be tuned to the particular application. 4,7 In applications that require elastomer components to operate in relatively harsh environments, for example, in oil wells, where temperatures and pressures can be as high as 500 K and 200 MPa, respectively, 3,8,9 hydrogenated NBR (HNBR) is often the preferred material. The hydrogenation process involves using catalysts to react NBR with hydrogen, thereby saturating the carbon−carbon double bonds present in the NBR backbone. 2,5 The resulting HNBR has enhanced properties compared to those of NBR: resistance to heat, blistering, and amine corrosion are improved, as is elongation to fracture, whereas other key properties that make NBR such an important elastomer in industrial application, such as oil resistance and low-temperature flexibility, are preserved. 2,3,10 Failure of critical elastomeric components can have very serious safety 11 and economic consequences. Performing controlled experiments under service conditions to assess performance and elucidate detailed failure mechanisms is often difficult, if not impossible. Conversely, computer simulations enable materials to be studied by varying the parameters of the model in a controlled and systematic manner. Often, these parameters are inaccessible experimentally and therefore simulations can provide information and insight that are complementary to those obtainable from experiments. However, despite the widespread use of NBR and HNBR, there are very few examples of atomistic simulations of these materials: Kucukpinar and Doruker 12 and Song et al. 13,14 used the all-atom COMPASS force field 15,16 to compute various properties of NBR chains up to a length of 164 monomers, including the solubilities and diffusivities of small gas molecules.
In this article, we present a new chemically inspired all-atom description of HNBR, that is, one that mimics the synthesis of HNBR via saturation of the carbon−carbon double bonds present in the backbone of NBR. An initial NBR structure is first generated by growing an NBR chain one bond at a time in a periodic simulation cell using a Monte Carlo approach. 17 It is then equilibrated with molecular dynamics (MD), 18 using the all-atom "Optimized Potentials for Liquid Simulations" (OPLS-AA) force field 19 to describe the interatomic interactions. The carbon−carbon double bonds in the NBR structure are then saturated by either cross-linking or hydrogenation, resulting in the generation of HNBR. Whether a particular double bond is cross-linked or hydrogenated is determined by two factors: (i) a user-defined parameter that sets the target density of the crosslinks in the structure and (ii) whether or not the local atomic environment of the bond permits a cross-link or, in other words, whether or not there is another carbon−carbon double bond available within a user-defined radius, with which to cross-link. The resulting HNBR structures are then also equilibrated with MD and the OPLS-AA force field. To the best of our knowledge, this is the first study of NBR and HNBR with the OPLS-AA force field and the first model of HNBR to be derived by direct and tunable saturation of NBR with a controllable cross-link density.
To test our models, we have calculated bulk properties, namely, mass density and glass-transition temperature, T g , for NBR and HNBR. We have compared our results to experimental results, where comparable experimental data were available in the literature. It is worth noting that although glass-transition temperature is one of the most important parameters for describing elastomers, 20,21 atomistic computer simulations struggle to reproduce it quantitatively, 22 as the time scales over which long-range coordinated motion of polymer chains occur is usually far beyond those accessible. An additional complication is that the measured value of T g is known to depend on the time scale of the experiment 20,23 and even different experimental methods give different results. 24 In light of these difficulties, we have (i) paid particular attention to the length and time-scale issues associated with calculating T g and (ii) limited our investigations to trends with respect to other crucial parameters, such as chemical composition (for NBR) and density of cross-links (for HNBR). In spite of the length and time-scale limitations of all-atom MD simulations, models such as those presented in this article are necessary for investigating the properties and trends that are dependent on the chemistry of the polymer, such as diffusivity and solubility of gas molecules, which are difficult to capture using coarsegrained approaches in which much of the detailed chemistry is lost.
The rest of this article is structured as follows: in the Method section, we describe our NBR and HNBR models in detail together with the computational procedures developed and used in this work; in the Result section, we present and discuss the results; and finally, we present our Conclusion.

■ METHODS
Structure of NBR and HNBR. NBR is a statistical copolymer consisting of two monomers: ACN and all isomers of butadiene (Figure 1), among which the trans configuration occurs most commonly. For example, in NBR, with a 33% monomer fraction of ACN, trans-butadiene accounts for 90% of the total butadiene content of the material. 3 Typical polymer chains consist of a sequence of ACN and butadiene, the ordering of which is determined by their reactivity ratios, 3 and typical chains are around 400−20 000 monomers in length. 25 The ratio of the constituent monomers determines the specific volume as well as the glass-transition temperature and other key properties. NBR is often used for oil and gas applications, and the ACN content is normally around 40% to obtain the desired oil resistance. Experiments conducted with un-cross-linked NBR (at 41% ACN) report a mass density of 1000 kg m −3 and a glass-transition temperature of 243 K. 25−27 HNBR is obtained from NBR by saturation of the double bonds present in the structure. The saturation level is usually close to 100%. 28 Hence, the key chemical feature that distinguishes HNBR from NBR is the absence of sp 2 -hybridized carbon atoms in the backbone of the polymer. Experimental studies of un-cross-linked HNBR (at 39% ACN) measure a density of 950 kg m −3 and glass-transition temperature of 248 K. 29,30 Model for NBR. To generate an initial structure for NBR, a 1000 monomer long polymer chain with a specified ACN/ butadiene ratio, using reactivity ratios of 0.02 and 0.28 for ACN and butadiene, respectively, is packed into a periodic simulation cell at a target mass density of 1000 kg m −3 to match the experimental density. This is done using a Monte Carlo process (as implemented in the Amorphous Cell module of BIOVIA Materials Studio 6.1 31 ), in which the polymer is grown one bond at a time, with fixed bond lengths and bond angles but variable dihedral angles for each newly added bond. The change in the potential energy, as calculated with the COMPASS force field, 15,16 is used to generate a self-avoiding torsionally biased random walk. 17 This procedure is repeated several times to generate an ensemble of initial configurations that samples the conformational phase space. As the method is inherently stochastic, each initial configuration generated will be distinct in the sense that they will be inaccessible to each other over the time scales achievable in atomistic MD simulations. The initial configurations are then equilibrated using MD and the OPLS-AA force field, as described later. Fillers as well as chemical cross-links and vinyl-or cis-butadiene (see Figure 1) are currently not included in our NBR model.
Model for HNBR. Our HNBR model is inspired by the actual process by which HNBR is synthesized from NBR, namely, saturation of the carbon−carbon double bonds in the backbone of NBR. Starting from an equilibrated NBR (39% ACN) structure, obtained as described above, all double bonds in the NBR are saturated either by hydrogenation or crosslinking. The process is controlled by two parameters, the crosslink probability, P CL , and the cross-link radius, R CL , and operates as follows. Each carbon−carbon double bond in the original NBR structure is chosen either to be cross-linked (with probability P CL ) or hydrogenated (with probability 1 − P CL ). In the former case, if there is no other carbon−carbon double bond available for cross-linking within radius R CL , then the double bond is simply saturated with hydrogen atoms. Otherwise, if there are one or more sites available for crosslinking within R CL , then the nearest site is the one that is chosen for cross-linking with. When two carbon atoms are cross-linked, the remaining two carbon atoms in the two participating double bonds are automatically saturated with hydrogen atoms to avoid relatively high energy cyclobutane-like The Journal of Physical Chemistry B Article configurations. Figure 2 illustrates the algorithm in the simple case of two unconnected butadiene monomers.
Once the process is completed for all double bonds in the polymer, the result is a fully saturated HNBR structure that has a cross-link fraction, f CL , defined as the fraction of double bonds of the original NBR structure that has been cross-linked as opposed to hydrogenated. In this work, we investigate the full range of values of P CL from 0% (i.e., hydrogenating every NBR double bond) to 100% (i.e., attempting to cross-link every double bond in the NBR structure). R CL is set to 5 Å, a value chosen to ensure that the resulting cross-link bonds are not overstretched. It is also worth noting that the saturation process is performed without any intermediate relaxation of the structure during the procedure; only at the end of the crosslinking process is the resulting HNBR structure equilibrated, as will be described later. The effect of permitting relaxation steps between cross-linking events, that is, a dynamical cross-linking approach, was also investigated. It was found that this did not have any effect on the mass density or distribution of lengths of the cross-linking bonds as compared to that in the static case. As a result, we adopted the static approach, as it is significantly more computationally efficient.
OPLS-AA Force Field and Equilibration. The interatomic interactions in our NBR and HNBR initial configurations are described using the OPLS-AA potential of Jorgensen et al. 19 for both bonded (bond stretching, bending, and twisting) and nonbonded (repulsion, dispersion, and Coulomb) interactions. In the OPLS force field, bonds and angles are described with harmonic potentials, whereas the dihedral terms are expressed as fourth-order Fourier expansions. For the repulsion and dispersion part, the 12-6 Lennard-Jones potential is used. The analytical form of the potential energy together with the parameters used can be found in the Supporting Information. It is worth noting that the functional form of the OPLS-AA force field makes the sensitivity analysis on the missing parameters rather straightforward because, in contrast to that for some of the other force fields that are commonly used for polymer simulations, such as PCFF+, there are no cross-terms between the different types of terms in the force field.
The bonded terms in OPLS-AA have been parameterized using ab initio simulations of small molecules, for example, butadiene and acetonitrile. Therefore, for a polymer composed of monomers based on these molecules, some parameters, such as those corresponding to the dihedral terms involving atoms at the bonding sites between chemically distinct monomers, may be missing from the force field. Among all of our simulations on different compositions and structures, at most 5% of the dihedral parameters were missing for any given system. To determine the impact of these parameters on observable quantities, such as the mass density, we performed a sensitivity analysis. For example, setting the missing C(sp 2 )−C(sp 3 )− C(sp 3 )−C(sp 2 ) dihedral parameters to be the same as the C(sp 3 )−C(sp 3 )−C(sp 3 )−C(sp 3 ) dihedral parameters versus setting them to zero resulted in a difference in the mass density of less than 1%, calculated as an average over eight independently generated and equilibrated configurations, which is within the standard deviation of each set of eight simulations. In the light of these tests, and for clarity and consistency, we chose to set all missing dihedral parameters to zero in all simulations presented in this article.
After the initial NBR and HNBR structures are generated, as described earlier, they are equilibrated broadly following the approach of Kucukpinar and Doruker, 12 with a few modifications; in particular, we used an equilibration procedure that is overall 10 times longer. The procedure is composed of the following six steps: (1) an energy minimization, (2) gradual compression from 0.1 to 500 MPa (to prevent the formation of large voids in the polymer matrix) and decompression back to 0.1 MPa at 298 K over 0.4 ns, (3) a second energy minimization, (4) annealing from 298 to 800 K and back to 298 K at 0.1 MPa over 1.2 ns, (5) a second compression/ decompression cycle as in step (2), and (6) a simulated anneal at 300 K at constant volume for 0.5 ns. A time step of 1 fs was used for all dynamical simulations, for which the energy drift was always less than 1% per 20 ns. The Nose−Hoover method The Journal of Physical Chemistry B Article was used for pressure and temperature control. 32 All dynamical simulations were performed with the Large-scale Atomic/ Molecular Massively Parallel Simulator (LAMMPS) code 18 at the Imperial College London high-performance computing facility. 33 Density and Glass Transition. Once the equilibrated structures were obtained, simulations were run to determine the mass density and glass-transition temperature for NBR as a function of ACN monomer fraction and for HNBR with a 39% ACN monomer fraction as a function of the density of crosslinks. In all calculations, eight independently generated (H)NBR structures were used to sample different parts of the configurational phase space. Each configuration is dynamically evolved, and snapshots were taken every 100 time steps (i.e., every 100 fs). The final numerical results reported are the average of all snapshots across all different structures evolved from the independently generated initial configurations.
To compute the mass density, we ran simulations with a constant number of particles, constant temperature (298 K), and constant pressure (0.1 MPa) for between 2 and 3 ns. The volume of the cubic computational cell was allowed to change isotropically. (Tests allowing anisotropic volume changes did not change the results quantitatively; however, they decreased the efficiency of the runs, requiring more computational time per calculation.) The glass-transition temperature, T g , was identified as the temperature at which there is a change in the slope of the volume of the computational cell as a function of temperature at constant pressure (Supporting Information). Starting at 500 K, the system was cooled by 15 K over 0.5 ns at a constant rate before an anneal for 2 ns to record the average volume at the lower temperature. This procedure was repeated in a stepwise manner until the temperature reached 150 K, generating 24 volume averages at 15 K intervals. The pressure of the system was maintained at 0.1 MPa. To determine T g , the 24 data points were divided into two subsets, [1, n) and [n, 24], for any n ∈ [2, 23], and a linear fit was computed on all possible combinations of subsets, that is for all possible n. The value of T g was calculated as the intersection of the two lines that gave the smallest combined coefficient of determination, R 2 (n), between all couples, following Simperler et al. 34 Heating the computational cell rather than cooling it showed no hysteresis in T g . The error bars on T g are the standard deviations on the computed quantities.

■ RESULTS AND DISCUSSION
Validation of the Force Field and Results for the Mass Density. Figure 3 shows the computed density of NBR as a function of ACN content compared with the experimental data. It is known experimentally that as the ACN content increases the hardness of NBR at room temperature increases. 35 For an ACN content higher than 50%, NBR ceases to be useful, and there is no experimental data available from 50 to 100% ACN NBR, that is, poly(acrylonitrile). Each computed point in Figure 3 is obtained from the average of eight independently generated configurations. The experimentally observed trend of increasing density as a function of ACN at 298 K and 0.1 MPa is correctly reproduced by the simulations. The underestimation of the mass density as a function of the ACN content is related to the transferability of the OPLS-AA force field. The parameters for poly(butadiene) and poly-(acrylonitrile) are taken from the parameters developed for bulk propene/ethane and acetonitrile/ethane, respectively. The results show that these parameters are exceptionally transferable to the case of polymerized butadiene and that the transferability to polymerized ACN, although not quite as exceptional, is still very good: for example, even for the extreme case of polyacrylonitrile (100% ACN), the discrepancy with the experimental density is only 5%.
For HNBR, we first consider the fully hydrogenated ( f CL = 0) configurations. By hydrogenating all double bonds in the NBR structure, we replace them with less stiff single bonds, and the added hydrogen atoms introduce more steric repulsion. As a result, we expect the fully hydrogenated configurations to have a lower density than that of NBR (for the same ACN content). This expectation is borne out by industrial data for HNBR before peroxide curing, that is, with no cross-links in the structure. Table 1 presents a comparison between these industrial data and our computed density. The agreement between the prediction of our model and the industrial value is within 3%, and moreover, the expected trend in density between NBR and HNBR is correct. By increasing the cross-link fraction in HNBR, we replace weak van der Waals forces between unbonded carbon atoms with shorter, permanent covalent single bonds. Therefore, we expect the density of HNBR to increase with the cross-link fraction. This expectation is borne out by the computations, as seen in Figure 4, in which mass density increases by up to 5%. There is a paucity of experimental data in the public domain to compare to, and where there is data available, it is for HNBR containing fillers, so no direct comparisons may be made.  The Journal of Physical Chemistry B

Article
Glass-Transition Temperature. The computed glasstransition temperature (T g ) of NBR as a function of ACN content is shown in Figure 5, where each point is the average T g obtained from eight independently generated configurations. Whereas there is a systematic overestimation of T g by the model, the trend for the experimental data is correctly described.
Furthermore, in Figure 6 we show T g as a function of crosslink fraction for HNBR (with 39% ACN content). The simulations confirm the expected trend that T g increases with increasing cross-link fraction, which is a result of the constraints of intrachain covalent bonds raising the temperature at which coordinated molecular motions can occur. 36 For HNBR no experimental results have been found in the public domain for T g as a function of cross-link density. However, a value of T g was stated in a supplier data sheet for 100% saturated HNBR (34% ACN) before peroxide curing 30 (i.e., comparable to the first point, f CL = 0, of the plot). Table 1 shows the comparison between the T g value obtained in this work and that from the data sheet. However, we have no information about the technique used to obtain the T g in the industrial data sheet and therefore the comparison must be treated with caution. Furthermore, even though we used Monte Carlo procedures to generate eight independent NBR starting structures, no improvement in the agreement with experimental data as well as in the error bars associated with the computed values was observed. To more closely investigate this limitation, we computed T g for 80 independently generated structures; the result can be seen in Figure 7. Given the spread of values obtained, even the use of 80 structures of NBR does not result in a better distribution of glass-transition temperatures.
In experiments, longer time scales compared to those for simulations can be investigated, normally on the order of seconds, minutes, hours, or even days. The typical cooling rate for experiments is 10 −1 −10 0 K s −1 , whereas for simulations is on the order of 10 10 −10 11 K s −1 ; hence, experiments allow for conformational changes that are unlikely to be effectively captured by MD. 21 This is reflected by the observation that T g changes by 3−5 K if the cooling rate changes by 1 order of magnitude. 20,37 Despite the limitations, our models successfully reproduce the expected trends for glass transition, as can be seen in Figures 5 and 6.    The Journal of Physical Chemistry B Article Further Remarks. As Materials Studio does not support the OPLS force field, either natively or through hard-coded parameters, the COMPASS force field was used to generate the initial NBR structures before switching to OPLS for the dynamical simulations. A Protein Data Bank file 38 is saved from the Materials Studio framework and read with a Python script that assigns the new parameters according to the OPLS-AA force field and writes a data file compatible with the LAMMPS engine. Switching between force fields during a simulation requires care, as not only can they be parametrized differently but they can also have different functional forms for the potential energy. In our case, however, the dihedral angles are the main degrees of freedom during the Monte Carlo process for generating the initial configurations (the bond lengths and angles are treated as fixed, as they are considerably stiffer). In both OPLS and COMPASS, the dihedral parameters are derived from ab initio simulations. Regardless of the particular density functional theory method used in the parametrization, qualitatively the same torsional energy landscape is expected to be obtained, for example, similar energetic barriers should be found for a twist of the butadiene configuration from a cis to a trans one. Therefore, we expect to generate similar structures by the Monte Carlo process, independently of the precise force field used.
The reason for using COMPASS and Materials Studio for the initial structure generation and OPLS-AA and LAMMPS for the MD simulations was to exploit the strengths of the two approaches while minimizing their weaknesses. On the one hand, Materials Studio 31,39 has a well-optimized Monte Carlo amorphous polymer builder that is ideal for generating the initial structure, whereas LAMMPS does not. On the other hand, the efficiency and scalability 40 of the LAMMPS package made it the preferred choice as the MD engine for this work. A further consideration was the open-source nature of LAMMPS and OPLS-AA, which enabled testing the sensitivity of properties on specific force-field parameters, the introduction of custom features to assist in data analysis, and the development of models that are generally more transparent and widely accessible. 19 The 1000 monomer long NBR structure is generated as a combination of butadiene and ACN in a ratio that can be specified. This chain length was chosen for computational efficiency, and although it is toward the lower end of the true range of chain lengths of NBR of 400−20 000, 25 it is nonetheless longer than the 164 monomer chain used by Kucukpinar and Doruker. 12 Furthermore, although the initial structure does not contain any vinyl-or cis-butadiene monomers (see Figure 1), this assumption is well supported by the relative ratios of the three butadiene monomers in real (H)NBR. 3 For the study of NBR, chemical cross-links are ignored; therefore, the only cross-linking is physical, that is, resulting from entanglement. Before being used in industrial applications, both NBR and HNBR are typically cross-linked and treated with additives such as carbon black, clay, and plasticizers, to improve the mechanical properties. 41,42 Fillers are usually nanoparticles with dimensions on the order of 100 nm and hence are difficult to simulate in a fully atomistic approach due to time-and length-scale limitations. 42,43 As a result, fillers are omitted in this work; a model able to take fillers into consideration would require coarse-graining of the system and would hence lead to a significant loss of chemical specificity.

■ CONCLUSIONS
In this paper, we have presented a novel all-atom description of HNBR that is inspired by the synthesis of HNBR from NBR. A Monte Carlo approach is used to stochastically generate and pack the initial NBR structure into a periodic box. The HNBR structure is then generated by opening the unsaturated NBR double bonds and either cross-linking or hydrogenating them, depending on a user-defined cross-link probability and the local atomic environment. The atomic interactions of both elastomers are described by the OPLS-AA force field. Our models were tested by computing the mass density and glasstransition temperature and comparing the results with experimental measurements, where available. The investigation on NBR focused on the trends with respect to the ratio between its two constituent monomers, whereas for HNBR the ACN percentage was kept constant at 39% and the trends as a function of the cross-link density of the system were studied.
The glass-transition temperature, T g , was found to increase as a function of ACN content for NBR and as a function of the cross-link density for HNBR. Whereas the qualitative agreement with the experimental data is good, the quantitative agreement is limited by the statistical uncertainties that highlight the sampling problems inherent in calculating T g using atomistic MD. Even by increasing the number of independently generated starting configurations by a factor of 10 (from 8 to 80), the spread of T g across the ensemble of configurations remains large. In addition to the neglect of fillers, a number of other factors make the direct quantitative comparison of our simulations of T g to experimental data problematic, including the use of periodic boundary conditions (which introduce spurious correlations in the model) and the single-chain approximation (which results in an underestimation of the free volume as compared to that of the physical system).
For the mass density, we find that it increases for NBR as a function of the ACN content and it increases as a function of the degree of cross-linking for HNBR. We also find that NBR has a 5% higher mass density than that of the uncross-linked ( f CL = 0) HNBR of the same ACN content (39%). The mass density of HNBR (with a fixed ACN content of 39%) can be systematically increased by up to 5% by tuning the degree of cross-linking. These results are in very good agreement with the experimental results, with only small differences of the order of a few percent.
In spite of their limitations, fully atomistic models are necessary for studying properties that depend on the chemistry of the polymer, such as diffusivity and solubility of gas molecules. These phenomena are difficult to describe accurately with coarse-grained models, in which much of the chemical specificity is lost. The models presented in this paper open up the possibility of investigating such properties for two of the most widely used elastomeric materials, namely NBR and HNBR.
Data underlying this article can be found on figshare.com at https://dx.doi.org/10.6084/m9.figshare.4026375, and used under the Creative Commons Attribution license. Copies of industry data sheets and other references that do not have permanent digital object identifiers are available on request from the corresponding author.
The Journal of Physical Chemistry B Article ■ ASSOCIATED CONTENT