Gaëtan Maurelab,
Florent Goujona,
Benoit Schnellb and
Patrice Malfreyt*a
aClermont Université, Université Blaise Pascal, Institut de Chimie de Clermont-Ferrand, ICCF, CNRS, UMR 6296, BP 10448, F-63000 Clermont-Ferrand, France. E-mail: Patrice.Malfreyt@univ-bpclermont.fr; Fax: +33 473 40 53 28; Tel: +33 473407204
bManufacture Française de Pneumatiques MICHELIN, Centre de Ladoux, 23, place des Carmes, 63040 Clermont-Ferrand, France
First published on 21st January 2015
We report mesoscale simulations of polymer melts and crosslinked polymer networks by using realistic coarse-grained (CG) models that are developed from atomistic simulations of polymer melts. We apply this multiscale strategy to different polymers in order to predict quantitatively some structural and thermomechanical properties such as the melt density, the end-to-end distance, the entanglement mass and the plateau modulus. The temperature dependence of the CG models is investigated through the calculation of the melt specific volumes at different temperatures and the calculation of the isothermal compressibility gives some insight into the pressure transferability of the CG models. We also show that the CG models can be applied successfully to high molecular weight chains. We test the performance of the CG models by calculating directly the plateau modulus of a crosslinked PIB network from mesoscopic simulations under a tensile stress. We compare the value of the plateau modulus with that calculated from the autocorrelation of the stress tensor during equilibrium simulations.
Indeed, a powerful investigation of the relationship between molecular structure and mechanical properties of polymer melts requires statistical averages on well-equilibrated melts of very long chains. For a polymer chain of length N, the longest relaxation of an entangled polymer2,17 scales as N3.4 with a computational cost growing as N.4 Atomistic simulations require solving the equations of motion on a timescale of 10−15 s and a length of Angstroms whereas collective relaxations can occur on a scale of micrometers and a time that can exceed 1 ms (see Fig. S1 (ESI†) for a comparison between atomistic and mesoscopic simulations for the calculation of the end-to-end vector autocorrelation function). To simulate materials on larger time scales, there is no alternative but to simplify the model. These simplified/coarse-grained (CG) models allow for longer length and time scales than atomistic models that are unable to reach the complete relaxation of the polymer melts7,8,18 with a reasonable computational effort. A number of CG models19–28 have been applied in the past to reproduce universal (scaling laws) properties of polymer melts but the quantitative prediction of properties as a function of the polymer chemistry remains much more challenging. The key problem is how to simplify the model without loosing the essential physics and being able to relate the parameters of the CG model to molecules of specific chemistry. Top-down and bottom-up parameterization schemes represent then two different alternatives to derive CG models. Top-down approaches derive parameters from macroscopic properties (compressibility, diffusion)29–31 and requires the use of a well-developed experimental database. Bottom-up approaches use the configurations at the atomistic level to develop interaction forces and parameters for mesoscopic model.5,9,16,18,32–34 It represents an attractive alternative for designing new polymeric materials from the structure–properties relationship.
Within the bottom-up approach, the Iterative Boltzmann Inversion (IBI) method is an iterative method for potential inversion from distribution functions using the potential of mean force approach.35 CG intramolecular and intermolecular potential functions can be then deduced from appropriate distribution functions of atomistic configurations.36–39 This method has been successfully applied to a variety of polymer melts formed by vinyl polymer chains,40 dendrimers,41 polystyrene7,8,32 and polyethylene chains.37–39 These coarse-grained models have been implemented either in Molecular Dynamics7,32,40–42 or in Dissipative Particle Dynamics (DPD)8,37–39 methods. There are some major advantages to use the DPD method: (i) in addition to the conservative force, the dissipative and random forces are short-ranged and pairwise additive so that the hydrodynamic interactions are preserved; (ii) the random and dissipative forces are independent of the conservative forces and are coupled to each other through the fluctuation–dissipation theorem to act as a thermostat. However, the use of IBI method for building CG potentials addresses a number of fundamental issues:7,40 (i) we have no guarantee that the potentials developed from atomistic configurations to match the distributions functions are able to reproduce all the properties of the original system; (ii) the way of grouping the atoms into the CG element and the degree of coarse-graining are not unique; (iii) the dynamics were often found to be significantly faster than the atomistic simulations.
In this paper, we propose to combine the IBI approach with the DPD method to reproduce thermomechanical and structural properties of different polymer melts (see Fig. 1). In a first step, we develop CG potentials for three polymers (cis-1,4-PB, PDMS and PIB) to calculate the density, end-to-end distance, entanglement mass and plateau modulus of melts. We examine the temperature and pressure transferabilities throughout the calculation of the coefficients of thermal expansion and isothermal compressibility, respectively. We also investigate the transferability to polymer chains of higher molecular weight. In a second step, we apply these potentials to different methodologies for the calculation of the plateau modulus of a crosslinked PIB network. We opt for a crosslinked network in order to avoid too large system-sizes that are expected to show a plateau in polymer melts.43–46 The plateau modulus is then calculated from non-equilibrium simulations by means of stress–strain curve and by equilibrium simulations using the stress tensor elements.
![]() | ||
Fig. 1 cis-1,4-Polybutadiene (cis-1,4-PB), poly(dimethylsiloxane) (PDMS), polyisobutylene (PIB) chemical structures along with the corresponding beads formed by five monomers. |
The outline of this paper is organized as follows. In Sec. II, we present the DPD method and the coarse graining procedure. Sec. III presents the results of this work concerning the calculation of the thermomechanical and structural properties of PIB, PB and PDMS melts and of the shear modulus of a crosslinked network of PIB. The main conclusions of this work are summarized in Sec. IV.
![]() | (1) |
![]() | (2) |
![]() | (3) |
fDij = −γωD(rij)(![]() ![]() | (4) |
![]() | (5) |
![]() | (6) |
The conservative potential wc(rij) sums intramolecular and intermolecular interactions. The equations of motion were integrated using a modified version of the velocity-Verlet (DPD-VV) algorithm.55 The development of integration schemes in DPD simulations is an area of active research.8,56–59 We have compared the modified velocity-Verlet scheme55 with that developed by Pagonabarraga et al.56 using a self-consistent technique to calculate forces and velocities (DPD-SC). With the values of time step used here, we did not detect any significant deviations between the target and the calculated temperatures calculated by the DPD-VV and DPD-SC algorithms.
w0bond(ri,i+1) = −kBT![]() ![]() | (7) |
w0bend(θi,i+1,i+2) = −kBT![]() ![]() | (8) |
w0nb(ri,j) = −kBT![]() ![]() | (9) |
From these w0bond(ri,i+1), w0bend(θi,i+1,i+2), w0nb(ri,j) coarse-grained potentials, we develop constant-DPD simulations of polymer melts formed by 40 polymer chains of 40 beads. A first DPD simulation of 10000 steps was performed with the initial w0bond(ri,i+1), w0bend(θi,i+1,i+2) and w0nb(ri,j) potentials. From these mesoscopic simulations, we re-calculate the corresponding radial distribution functions gbondn(ri,i+1), gbendn(θi,i+1,i+2) and gnbn(ri,j), respectively. We observe that the target RDFs in eqn (7), (8) and (9) are not accurately reproduced. By calculating the ratio between these intermediary RDFs and the target RDFs according to the Iterative Boltzmann Inversion process, we obtain new potential functions as
![]() | (10) |
An iteration consists of a constant-NVT DPD simulation of 10000 steps and 5 iterations were required to obtain a good convergence. Additional constant-NPT DPD simulations were performed to modify these intermolecular potentials by using an attractive linear function18 in order to make the pressure in line with a pressure of 0.1 MPa. Fig. 2a shows that the intramolecular target RDFs obtained from atomistic simulations are very well reproduced by the CG models after the iterative process for the cis-1,4-PB, PIB and PDMS polymers. The resulting CG potential functions used in the DPD simulations are given in part b of Fig. 2. We observe some significant changes in the distribution functions between the different polymers indicating that the presence of the methyl groups and their resulting excluded-volume in the PDMS and PIB chains are well-reproduced by oscillations in the distribution functions in spite of a coarse-graining degree of five. Fig. 2c and e show the bonding and bending distribution functions and the corresponding CG potentials are represented in Fig. 2d and f, respectively. It is worth noting that the CG potential for the bond stretching resembles a simple harmonic potential that has been widely applied to model the interactions between consecutive beads36,50,51,61 at the mesoscale. The deviation from this simple harmonic potential are polymer-microstructures dependent. From Fig. 2f, we observe that the CG bending potentials tend to disfavor the small angles in line with the bending potentials used for flexible and semi-flexible polymer chains.46
Since the dissipative forces are present in the simulations, another issues concern the choice of the friction coefficient7,38,39 that has been shown to significantly impact on the dynamical properties of the polymer chains. It is well-known that the dynamical properties of polymer melts can be significantly affected by the choice of these adjustable parameters8,38 rc and γ. Based on the work38 of Lahmar and Rousseau who have studied the influence of rc and γ on the global and local dynamics of a polymer melt, we checked that γ = 300 kg mol−1 ns−1 was a good choice for the observation of long-time dynamical processes. Indeed, by using this value of γ, we have shown in a previous work18 that the diffusion coefficient of the cis-1,4-PB is slightly faster than that calculated from atomistic simulations but it remains on the same order of magnitude.
![]() | ||
Fig. 3 (a) Simulated and experimental polymer melt densities as a function of time for cis-1,4-PB, PIB and PDMS polymers; (b) mean square end-to-end distances as a function of time. The experimental data, shown in dotted lines in (a) and (b) are taken from ref. 62. |
Me (g mol−1) | ρ (g cm−3) | (G0N) (MPa) | Ree2/M (Å2 mol g−1) | Ree2/RG2 | αp (10−4 K−1) | ||||||
---|---|---|---|---|---|---|---|---|---|---|---|
CG | Exp. | CG | Exp. | CG | Exp. | CG | Exp. | CG | Exp. | ||
cis-1,4-PB | 2500 | 2347 | 0.9109 | 0.900 | 0.731 | 0.76 | 0.788 | 0.76 | 6.1 | 6.51 | 6.7 |
PIB | 6238 | 5686 | 0.9255 | 0.918 | 0.301 | 0.32 | 0.615 | 0.57 | 6.2 | 5.21 | 5.5 |
PDMS | 10588 | 9613 | 0.9735 | 0.970 | 0.181 | 0.20 | 0.423 | 0.42 | 6.2 | 5.61 | 9.0 |
We now focus on the transferability of the CG potentials to the polymer chain length, temperature and pressure. Whereas the CG potentials have been developed from simulations of polymer chains of 200 monomers (40 beads), we use these potentials to simulate larger polymer chains of 800 monomers corresponds to a length of 160 beads. The polymer melt density and the end-to-end distances are given in Table S1 (ESI†). The deviations from experiments are not larger than those calculated with smaller chain lengths. This is an interesting result that shows the possibility of developing the CG potentials from shorter polymer chains. We now turn attention to the temperature transferability of the CG potentials by calculating the specific volume of the polymer melt in the [200–500] temperature range above the glass transition temperature. The simulated specific volumes result from constant-NPT simulations and their evolution with the temperature are represented in Fig. 4 at different temperatures. The linear temperature-dependence of the specific volume on the temperature allows to determine the thermal expansion coefficient defined by eqn (11). The value of the thermal expansion, given in Table 1, is calculated from the slope of the different lines at 300 K. We also report for comparison the experimental17 and simulated values of αp at 300 K by dotted and solid lines in Fig. 4, respectively. We observe a good agreement between experiments and simulations for cis-1,4-PB and PIB polymers with a deviation of less of 5%. The deviation is more significant with PDMS leading to a deviation of 35%. We can conclude however that the order of magnitude for αp is reproduced with the CG models developed in this work. A better agreement between experimental and calculated values of αp for PDMS is possible but would require some refinements of the CG potential. The temperature dependence of the polymer melt density is then well-reproduced by the CG models in the temperature range above the glass transition temperature.
![]() | (11) |
To complete the analysis, we investigate the pressure dependence of these CG potentials through the calculation of the isothermal compressibility κT from the volume fluctuations in the constant-NPT ensemble as
![]() | (12) |
The CG models give values of κT at 300 K equal to 1090 × 10−15, 130 × 10−15, 120 × 10−15 bar−1 for cis-1,4-PB, PIB and PDMS polymers, respectively. The corresponding experimental values17 for cis-1,4-PB, PIB and PDMS polymers are 7.2 × 10−15, 4.8 × 10−15, 11 × 10−15 bar−1, respectively. This comparison establishes the inability of the CG potentials to reproduce any pressure dependence leading to deviations from experiments of two orders of magnitude as already observed for other CG models.32,64 The procedure applied here allows to reproduce some thermodynamic properties at a target pressure but excludes any transferability to other pressures due to the deficiency of reproducing the incompressibility of the polymer material at equilibrium.
The explanation of this poor reproduction of the incompressibility comes from various origins: (1) the degree of coarse-graining, (2) the correction of the pressure during the development of the CG potentials. First, concerning the degree of coarse-graining, we develop CG potentials for a level of coarse-graining of 4 leading to harder potentials. The results are presented in Table S1 (ESI†). Fig. S4 (ESI†) shows that the trajectory of the simulated polymer density shows weaker fluctuations of the density than with λ = 5. This clearly establishes a relationship between the softness of the potential and the compressibility property. Second, the modification of the CG potentials to reproduce the pressure of the atomistic configurations introduces a pressure dependence in these CG potentials.32,64
We propose here to analyze the CG configurations of the different polymer melts by applying the primitive path analysis (PPA65) methodology, thoroughly described in ref. 18. Fig. 5 shows the distributions of the molecular weight (Me) between entanglements. The average entanglement mass calculated from this distribution is reported in Table 1 along with the experimental corresponding property. We establish here a very good performance of the prediction with a maximum deviation of 10% from experiments. Additionally, the differences between the entanglement molecular weight of each polymer are also very well reproduced: the ratios MPDMSe/MPIBe and MPDMSe/Mcis-1,4-PBe are equal to 1.7 and 4.2 for CG simulations against 1.7 and 4.1 for experiments, respectively. For PIB, an average entanglement mass of 6238 g mol−1 corresponds to an entanglement length of 9 beads (45 monomers) and to the presence of approximately 17 entanglement segments per chain formed by 160 beads. As a result, we model here polymer chain lengths on the order of 20 entanglement lengths Ne whereas the longest chains that can be simulated by atomistic models within a reasonable computational cost represent roughly 2Ne. Since the plateau modulus G0N is related to Me through the relationship62 , it is then possible to estimate the plateau modulus G0N of linear polymer melts. However, this plateau modulus could be calculated directly from CG simulations without making any assumption in any theoretical model but would require very long chains. The values of G0N are reported in Table 1 and the agreement with experiments is quantitative with a maximum deviation of 10% for PDMS and 4% for cis-1,4-PB.
![]() | ||
Fig. 5 Distributions of the entanglement molecular weight for cis-1,4-PB, PIB and PDMS polymers. The dotted curves represent the fits obtained from a Poisson distribution. |
![]() | (13) |
![]() | (14) |
σzz = Pisotropiczz − Pzz | (15) |
Fig. 6a shows the stress–strain curve in a wide range of deformations from the elastic regime to the high elongation regime. We also plot the stress–strain curves deduced from the classical theories, the phantom network model67 and the affine deformation model68,69 of the rubber elasticity. We observe that the curve deviates from the phantom behavior for an elongation ratio of εz = 0.6 leading to an hardening of the network in the regime of strong deformations. We also observe that the stress–strain curve of the crosslinked PIB network exhibits a linear portion in the range of weak deformations and matches very well the phantom model. The agreement between the simulated and theoretical curves of the phantom model is not surprising since our network corresponds to a weakly constrained system due to the absence of additional topological constraints (no entanglement) between two crosslinks whereas the affine model describes a strongly constrained network. For strong deformations (Fig. 6a), the stress–strain curve deviates from both the affine and phantom theories as expected for the hardening of the network. In the linear regime (see Fig. S6 (ESI†)), the stress–strain curve allows to determine the Young's modulus E for uniaxial deformation where E is defined as the slope of where (ε = λ − 1) → 0. We find that E = 0.57 MPa. Since G = E/3 in the elastic regime for incompressible systems, i.e., with a Poisson's ratio of 0.5 (see Fig. 6b and S5 (ESI)†), we obtain a numerical value of the plateau modulus G0N of 0.19 MPa in line with that calculated from
MPa.
It should be noticed that the CG models are able here to reproduce the Poisson's coefficient characterizing incompressible systems whereas these same models were unable to reproduce the incompressibility of polymer melts (see Section 3.1). It means that adding topological constraints improves the description of the compressibility behavior. It seems that decreasing the interpenetration between the beads leads to a better reproduction of the compressibility in line with the impact of harder potentials on the density fluctuations of pure polymer melts. In the case of a crosslinked network of polymer, we can conclude that the CG models developed here perform very well in the reproduction of the Young's modulus and plateau modulus through the calculation of the stress–strain curve.
Another alternative for a direct calculation of the plateau modulus from equilibrium simulations is to consider the shear relaxation modulus G(t) from the autocorrelation of the stress tensor elements44 by running very long equilibrium simulations. The shear relaxation modulus can be calculated from the autocorrelation of the stress tensor elements44 as
![]() | (16) |
![]() | (17) |
Fig. 7 reports the calculated stress relaxation function of the crosslinked PIB network. As expected, the G(t) curve exhibits a plateau region. Very interestingly, the value of the plateau modulus fits very well with the values predicted by the nonequilibrium CG simulations and the phantom theory. We show here in the case a crosslinked network that the two different methodologies converge to the same value of the plateau modulus establishing that the combination of the CG models with the DPD method as an efficient tool for determining this property.
We have developed CG models of different polymers (cis-1,4-PB, PIB and PDMS) by using the iterative Boltzmann inversion method. Nonbonded and intramolecular (bonding and bending) potentials have been built from distribution functions of the atomistic configurations. We took the route of developing the CG models in the constant-NPT ensemble in order to reproduce the pressure of the atomistic simulations. These models have been incorporated in the DPD method. We did not aim to focus on the dynamical properties of the polymer melt since we have shown recently that the dynamics of our CG models18 is slightly faster than that of the atomistic models.
Here, we have established that the (IBI + DPD) combined approach can be applied to different polymers and different properties. We have shown that the CG models are able to reproduce accurately the polymer melt density, the end-to-end distance, the entanglement mass of different polymer melts. From theoretical models, it is then possible to estimate the plateau modulus with a reasonable computational effort without running simulations of very long chains. The temperature and pressure transferabilities of these potentials have been discussed. It results a relatively good reproduction of the thermal expansion coefficient whereas the calculation of the isothermal compressibility shows significant deviations of two orders of magnitude from experiments. We have proposed some assumptions that could explain this behavior requiring however further investigations.
In order to perform a direct calculation of the plateau modulus with these CG models within a reasonable system-size, we have opted for a crosslinked network of PIB. We have shown that the two ways of calculating the plateau modulus either by applying a deformation either by calculating the autocorrelation of the stress tensor elements lead to the value expected by the classical theories of linear elasticity. This work underlines the performance of the bottom-up approach for quantitative predictions of thermomechanical properties of polymer with a reasonable computational cost.
Footnote |
† Electronic supplementary information (ESI) available: The atomistic and CG models are compared on their ability of simulating well-equilibrated melts of long polymer chains. The results with a smaller degree of coarse-graining and different polymer chain lengths are given for comparison. The methodology corresponding to the simulation of the elongation of the crosslinked polymer network is described. See DOI: 10.1039/c4ra16417b |
This journal is © The Royal Society of Chemistry 2015 |