Prediction of structural and thermomechanical properties of polymers from multiscale simulations

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

Received 15th December 2014 , Accepted 20th January 2015

First published on 21st January 2015


Abstract

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.


1 Introduction

The improvement of the desired properties and performance of a polymeric material requires investigating the relationship between the structure at the atomic or molecular length and the key physical and chemical macroscopic properties. Although molecular theories developed in polymer science have made paramount progress in relating the structure of polymer chains to their macroscopic properties, such theoretical methods suffer from a limiting predictive power due to the large variety of chemical constituents and molecular architectures that are at the origin of the viscoelastic properties. Molecular modelling1–3 has become an efficient tool for checking the theoretical predictions of properties, for giving a molecular description of the analysis of experimental work and for designing advanced materials with specific properties. However, many problems at the leading edge of materials science involve collective phenomena that occur over a range of time and length scales that are difficult to capture in atomistic simulations. As a result, the ability to perform computer simulations of materials over length scales that are relevant to experiments represents a grand challenge in computational material science. Different approaches have been developed to model the matter at two or more different length and time scales and participate to the overall strategy of multiscale modelling.4–16

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.


image file: c4ra16417b-f1.tif
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.

2 Developments of the coarse-grained (CG) potentials

2.1 Atomistic molecular dynamics (MD) simulations

The atomistic simulations of the cis-1,4-PB, PDMS and PIB polymers were performed using the all-atom representation of the COMPASS force field47 and the Forcite module from MATERIALS STUDIO Accelrys package.48 The bulk atomistic configurations consisted of 40 chains of 200 monomers. Constant-NPT simulations were carried out at 298 K and 0.1 MPa. The equations of motion were integrated using the Verlet leapfrog algorithm scheme with a timestep of 2 fs. The cutoff radius for the Lennard–Jones interactions was fixed to 14 Å. A typical MD simulation consists of a first period of 10 ns followed by an additional production phase of 50 ns. The dimensions of the simulation cell were chosen in order to avoid self-entanglements. The validation of these atomistic force field was carried out by comparing the polymer melt density and its temperature dependence as well as the end-to-end distance. We have also checked some methodological aspects of the atomistic simulations (see Fig. S2 (ESI)) establishing that the description of the polymer melts is as good as possible with an atomistic model. As an example, for the cis-1,4-PB polymer, the atomistic simulations give an end-to-end distance of 0.74 ± 0.03 Å2 mol g−1 for an experimental value of 0.76 Å2 mol g−1 with a ratio to the mean square end-to-end distance to the mean square radius of gyration equal to 6.1. These atomistic simulations18 have established a quantitative prediction of the thermal expansion coefficient of 6.5 × 10−4 K−1 against an experimental value of 6.7 × 10−4 K−1. Once this atomistic force field was validated, the bottom-up strategy used here consists of designing CG potentials from atomistic molecular dynamics simulations.

2.2 Dissipative particle dynamics (DPD)

The dissipative particle dynamics49–53 method solves the Newtonian equations for particles subject to conservative, dissipative and random forces. Thus,
 
image file: c4ra16417b-t1.tif(1)
where fi is a pairwise-additive force defined as the sum of three contributions
 
image file: c4ra16417b-t2.tif(2)
where fCij, fRij and fDij are the conservative, random and dissipative forces, respectively. The conservative repulsive force fCij derives from a soft interaction potential and is expressed as
 
image file: c4ra16417b-t3.tif(3)
where wc(rij) is the conservative potential, rij is the relative displacement of particles i and j and [r with combining circumflex]ij is the corresponding unit vector. rc is the cutoff radius and is fixed to 4.0 nm. The dissipative and random forces are given by
 
fDij = −γωD(rij)([r with combining circumflex]ij·vij)[r with combining circumflex]ij (4)
 
image file: c4ra16417b-t4.tif(5)
where δt is the time step. vij = vivj is the relative velocity, the terms ωD(rij) and ωR(rij) are dimensionless weighting functions. σ is the amplitude of the noise, θij is a random Gaussian number with zero mean and unit variance. γ and σ are the dissipation strength and noise strength, respectively. Español and Warren54 have shown that the system will sample the canonical ensemble and obey the fluctuation-dissipation theorem if the following conditions are satisfied.
 
image file: c4ra16417b-t5.tif(6)
where kB is Boltzmann's constant, T is the temperature and ωD(rij) = (1 − r/rc).2

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.

2.3 Coarse-graining procedure

From 10 independent atomistic configurations, we built bead–bead intermolecular and intramolecular potential functions calculated by coarse-graining atomistic configurations. This procedure allows a direct mapping between the atomistic and CG systems. The degree of coarse-graining was fixed to 5, indicating that each bead corresponds to 5 monomers. This degree of coarse-graining is then much smaller than the entanglement mass and allows the use of a relatively larger time step. By using a value of 5, no bond crossings are detected during the dynamics of the polymer chains. As a consequence, we expect a good reproduction of the entanglement effects. The time step was fixed to 50 fs and the temperature was kept to 300 K. Three different distributions functions are developed simultaneously: the radial distribution function gbond(ri,i+1) between two consecutive beads in the polymer chain; the radial distribution function gbend(θi,i+1,i+2) between three consecutive beads in the polymer chain; the radial distribution function gnb(ri,j) between beads i and j of different polymer chains and beads i and j of the same molecule separated by more than one bond. From these radial distribution functions (RDF), we use the potential of mean force approach35 to derive the corresponding CG potential functions as
 
w0bond(ri,i+1) = −kBT[thin space (1/6-em)]ln[thin space (1/6-em)]gbond(ri,i+1) (7)
 
w0bend(θi,i+1,i+2) = −kBT[thin space (1/6-em)]ln[thin space (1/6-em)]gbend(θi,i+1,i+2) (8)
 
w0nb(ri,j) = −kBT[thin space (1/6-em)]ln[thin space (1/6-em)]gnb(ri,j) (9)
where w0bond(ri,i+1), w0bend(θi,i+1,i+2), w0nb(ri,j) are the CG potentials corresponding to the bonding, bending and nonbonded interactions developed from the CG atomistic configurations. To improve the statistics, we have adopted the strategy37 consisting of using different ways of grouping the five monomers in the atomistic configurations. The resulting potential functions have been tabulated: the force are then obtained by a cubic spline interpolation.60

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 10[thin space (1/6-em)]000 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

 
image file: c4ra16417b-t6.tif(10)
where gn(r) and wn(r) refer to the bonding, bending and non-bonded distribution and potential functions at the step n, respectively. g(r) is the target RDF. This scheme is reiterated until convergence is obtained.

An iteration consists of a constant-NVT DPD simulation of 10[thin space (1/6-em)]000 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


image file: c4ra16417b-f2.tif
Fig. 2 (a) Nonbonded radial distribution functions calculated from MD and DPD simulations and (b) nonbonded CG potential functions used in the DPD simulations obtained from the CG procedure. (c) Bonding and (e) bending CG distribution functions built from MD and DPD simulations. Corresponding (d) bonding and (f) bending CG potential functions obtained for a degree of coarse-graining of 5 monomers.

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.

3 Discussions

3.1 Prediction of the thermomechanical and structural properties

It is now well-established that the combination of the CG models with DPD is able to equilibrate polymer melts containing high molecular mass chains18,32 with a significant computational reduction (see Fig. S1 and S3 (ESI)). Fig. 3a shows the density of the different polymer melts as a function of time. These trajectories resulting from stable constant-NPT simulations allow to extract time average densities that are given in Table 1 for each polymer. Interestingly, the CG potentials allow a very good reproduction of the melt density with a maximum deviation of 1% from experiments.62 Additionally, the change in the density due to the chemical nature of the monomers is also very well predicted by these mesoscale simulations. Fig. 3b shows the end-to-end distance of the polymer chains versus time. The different trajectories show that the CG models are able to distinguish the different polymers on this structural property even if the chemical details of the polymer are no longer explicitly described. The comparison with the experimental end-to-end distance shows a quantitative prediction with a maximum deviation of 7% for PIB (see Table 1). The ratio of the end-to-end distance to the radius of gyration, given for each polymer in Table 1 is in line with the expected value of 6 for ideal gaussian chains.63 The bottom-up approach developed in the constant-NPT ensemble becomes predictive for the density and reproduces also accurately the end-to-end distance.
image file: c4ra16417b-f3.tif
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.
Table 1 Entanglement molecular mass (Me), melt density(ρ), plateau modulus (G0N), end-to-end distance (Ree/M), ratio to the mean square end-to-end distance to the mean square radius of gyration, thermal expansion coefficient αp calculated from CG DPD simulations. The experimental values for each polymer are taken from ref. 62. The subscripts give the accuracy of the last decimal(s), i.e., 0.788 means 0.78 ± 0.08
  Me (g mol−1) ρ (g cm−3) (G0N) (MPa) Ree2/M2 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.

 
image file: c4ra16417b-t7.tif(11)


image file: c4ra16417b-f4.tif
Fig. 4 Specific volume (1/ρ) of the three polymer melts as a function of temperature calculated using the CG models where ρ is the average polymer melt density. The solid and dotted lines represent to the simulated and experimental values of αp at 300 K.

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

 
image file: c4ra16417b-t8.tif(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 image file: c4ra16417b-t9.tif, 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.


image file: c4ra16417b-f5.tif
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.

3.2 Application of an uniaxial deformation

We now turn attention to the transferability of the CG models on the direct calculation of the modulus G0N of a crosslinked network of PIB. A crosslinked network allows to observe a plateau regime within a reasonable computational effort whereas the same calculation for polymer melts would require very large chains of polymers. We consider then a randomly tetrafunctional crosslinked PIB network where the crosslinking is ensured by an insertion of a few amounts of isoprene units which allows chemical bonding between chains through vulcanization. We have checked that the CG potential is not affected with a degree of coarse-graining of 5. The average molecular weight (Mc = 6379 g mol−1) between crosslink units is approximately the entanglement molecular weight of the PIB melt leading to no trapped entanglement between crosslinks. In a first step, we aim to calculate the stress–strain curve of the network by applying a tension in the z-direction and by maintaining a constant pressure in the x and y-directions. More precisely, uniaxial tensile is applied by stretching the z-dimension of the simulation box under a negative normal pressure Pzz. We also apply a positive pressure Pxx = Pyy = 0.1 MPa in the x and y-directions. We used the anisotropic Berendsen barostat66 for the control of the pressure. The magnitude of the elongation along the z-axis is measured through
 
image file: c4ra16417b-t10.tif(13)
or
 
image file: c4ra16417b-t11.tif(14)
where Lz and Lz,0 are the box dimensions along the z-axis at times t and t = 0, respectively. This property is reported in Fig. 6a at three different values of stress tensor defined as
 
σzz = PisotropiczzPzz (15)
where Pisotropiczz = 0.1 MPa was the value of the normal component in the system when the pressure was isotropically maintained. As Pzz varies from −0.2 to 0, σzz decreases from 0.3 to 0.1 MPa.

image file: c4ra16417b-f6.tif
Fig. 6 Stress–strain curve of uniaxial deformation as a function of (λλ−2) over (a) a large range of deformations and (b) a limited range of deformations obtained from CG simulations, affine and phantom models. (c) Stress–strain curve as a function of ε in order to focus on the region of weak deformations and on the linear behavior (black dotted line) of the fitting curve (red curve).

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 image file: c4ra16417b-t12.tif 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 image file: c4ra16417b-t13.tif 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

 
image file: c4ra16417b-t14.tif(16)
with Nαβ = σαασββ where α, β denote the x, y and z components. The components of the stress tensor are defined by
 
image file: c4ra16417b-t15.tif(17)
where fC,αij is the α component of the conservative force defined in eqn (3) and vαi the α component of the velocity of bead i.

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.


image file: c4ra16417b-f7.tif
Fig. 7 Time evolution of the shear modulus G calculated from the stress relaxation function.

4 Conclusions

One of the challenge of the CG models in reducing the number of degrees of freedom is to incorporate the chemical nature of the molecules into the CG element. This is an essential step toward the development of realistic effective potentials. Additionally, the bottom-up approach addresses fundamental questions concerning the thermodynamic and structural consistencies of the CG models with the underlying atomistic configurations. This consistency is required for a quantitative prediction of thermomechanical, structural and viscoelastic properties of polymeric materials that are often only accessible at the mesoscopic scale.

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.

References

  1. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1989 Search PubMed.
  2. Monte Carlo and Molecular Dynamics Simulations in Polymer Science, ed. K. Binder, Oxford University Press, New York, 1995 Search PubMed.
  3. Simulation Methods for Polymers, ed. M. Kotelyanskii and D. N. Theodorou, Marcel Dekker, Inc., New York, 2005 Search PubMed.
  4. F. Müller-Plathe, ChemPhysChem, 2002, 3, 754–769 CrossRef.
  5. V. A. Harmandis, N. P. Adhikari, N. F. A. van der Vegt and K. Kremer, Macromolecules, 2006, 39, 6708–6719 CrossRef.
  6. C. Peter and K. Kremer, Soft Matter, 2009, 5, 4357–4366 RSC.
  7. H. J. Qian, P. Carbone, X. Chen, H. A. Karimi-Varzaneh, C. C. Liew and F. Müller-Plathe, Macromolecules, 2008, 41, 9919–9929 CrossRef CAS.
  8. H. J. Qian, C. C. Liew and F. Müller-Plathe, Phys. Chem. Chem. Phys., 2009, 11, 1962–1969 RSC.
  9. V. A. Harmandis and K. Kremer, Macromolecules, 2009, 42, 791–802 CrossRef.
  10. D. Fritz, K. Koschke, V. A. Harmandaris, N. F. A. van der Vegt and K. Kremer, Phys. Chem. Chem. Phys., 2011, 13, 10412–10420 RSC.
  11. Y. Li, S. Tang, B. C. Abberton, M. Kröger, C. Burkhart, B. Jiang, G. J. Papakonstantopoulos, M. Poldneff and W. K. Liu, Polymer, 2012, 53, 5935–5952 CrossRef CAS PubMed.
  12. Y. Li, B. C. Abberton, M. Kröger and W. K. Liu, Polymers, 2013, 5, 751–832 CrossRef PubMed.
  13. R. Potestio, S. Fritsch, P. Espanol, R. Delgado-Buscalioni, K. Kremer and R. E. D. Donadio, Phys. Rev. Lett., 2013, 110, 10831 CrossRef.
  14. K. Johnston and V. Harmandaris, Soft Matter, 2013, 9, 6696–6710 RSC.
  15. K. Johnston and V. Harmandaris, Macromolecules, 2013, 46, 5741–5750 CrossRef CAS.
  16. S. Trément, B. Schnell, L. Petitjean, M. Couty and B. Rousseau, J. Chem. Phys., 2014, 140, 134113 CrossRef PubMed.
  17. J. E. Mark, K. Ngai, W. Graessley, L. Mandelkern, E. Samulski, J. Koening and G. Wignall, Physical Properties of Polymers, Cambridge University Press, Cambridge, 3rd edn, 2004 Search PubMed.
  18. G. Maurel, B. Schnell, F. Goujon, M. Couty and P. Malfreyt, J. Chem. Theory Comput., 2012, 8, 4570–4579 CrossRef CAS.
  19. G. S. Grest and K. Kremer, Phys. Rev. A, 1986, 33, 3628–3631 CrossRef CAS.
  20. A. Baumgärter and K. Binder, J. Chem. Phys., 1981, 75, 2994–3005 CrossRef PubMed.
  21. Computational Modelling of Polymers, ed. K. Binder and M. Dekker, New-York, 1992 Search PubMed.
  22. G. S. Grest, Curr. Opin. Colloid Interface Sci., 1997, 2, 271–277 CrossRef CAS.
  23. K. Kremer and K. Binder, Comput. Rep., 1988, 7, 259 CrossRef CAS.
  24. K. Kremer and G. S. Grest, J. Chem. Phys., 1990, 92, 5057–5086 CrossRef CAS PubMed.
  25. G. S. Grest, Adv. Polym. Sci., 1999, 138, 149 CrossRef CAS.
  26. T. Kreer, M. H. Muser, K. Binder and J. Klein, Langmuir, 2001, 17, 7804–7813 CrossRef CAS.
  27. J. T. Padding and W. Briels, J. Chem. Phys., 2001, 115, 2846–2859 CrossRef CAS PubMed.
  28. J. T. Padding and W. J. Briels, J. Chem. Phys., 2002, 117, 925–943 CrossRef CAS PubMed.
  29. R. D. Groot and K. L. Rabone, Biophys. J., 2001, 81, 725–736 CrossRef CAS.
  30. A. Ghoufi and P. Malfreyt, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 051601 CrossRef CAS.
  31. A. Ghoufi and P. Malfreyt, J. Chem. Theory Comput., 2012, 8, 787–791 CrossRef CAS.
  32. T. Spyriouni, C. Tzoumanekas, D. Theodorou, F. Müller-Plathe and G. Milano, Macromolecules, 2007, 40, 3876–3885 CrossRef CAS.
  33. T. Mulder, V. A. Harmandaris, A. V. Lyulin, N. F. A. van der Vegt, K. Kremer and M. A. J. Michels, Macromolecules, 2009, 42, 384–391 CrossRef CAS.
  34. L. Lu, J. F. Dama and G. A. Voth, J. Chem. Phys., 2013, 139, 121906 CrossRef PubMed.
  35. J. G. Kirkwood, J. Chem. Phys., 1935, 3, 300–313 CrossRef CAS PubMed.
  36. D. Reith, M. Pütz and F. Müller-Plathe, J. Comput. Chem., 2003, 24, 1624–1636 CrossRef CAS PubMed.
  37. X. Guerrault, B. Rousseau and J. Farago, J. Chem. Phys., 2004, 121, 6538–6546 CrossRef CAS PubMed.
  38. F. Lahmar and B. Rousseau, Polymer, 2007, 48, 3584 CrossRef CAS PubMed.
  39. F. Lahmar, C. Tzoumanekas, D. N. Theodorou and B. Rousseau, Macromolecules, 2009, 42, 7485–7494 CrossRef CAS.
  40. G. Milano and F. Müller-Plathe, J. Phys. Chem. B, 2005, 109, 18608–18619 CrossRef PubMed.
  41. P. Carbone, F. Negri and F. Müller-Plathe, Macromolecules, 2007, 40, 7044–7055 CrossRef CAS.
  42. P. Carbone, H. Varzaneh and F. Müller-Plathe, J. Chem. Phys., 2008, 128, 064904 CrossRef PubMed.
  43. O. Byutner and G. S. Smith, Macromolecules, 2002, 35, 3769–3771 CrossRef CAS.
  44. A. E. Likhtman, S. K. Sukumaran and J. Ramirez, Macromolecules, 2007, 40, 6748–6757 CrossRef CAS.
  45. W. B. Lee and K. Kremer, Macromolecules, 2009, 42, 6270–6276 CrossRef CAS.
  46. A. E. Likhtman and S. K. Sukumaran, Macromolecules, 2010, 43, 3980–3983 CrossRef CAS.
  47. H. Sun, J. Phys. Chem. B, 1998, 102, 7338–7364 CrossRef CAS.
  48. Accelrys Software Inc., Material Studio Modeling Environment, Release 6.2, Accelrys Software Inc., San Diego, 2007 Search PubMed.
  49. P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett., 1992, 19, 155–160 CrossRef.
  50. P. Malfreyt and D. J. Tildesley, Langmuir, 2000, 16, 4732–4740 CrossRef CAS.
  51. F. Goujon, P. Malfreyt and D. J. Tildesley, Macromolecules, 2009, 42, 4310–4318 CrossRef CAS.
  52. C. Ibergay, P. Malfreyt and D. J. Tildesley, J. Chem. Theory Comput., 2009, 5, 3245–3259 CrossRef CAS.
  53. C. Ibergay, P. Malfreyt and D. J. Tildesley, J. Phys. Chem. B, 2010, 114, 7274–7285 CrossRef CAS PubMed.
  54. P. Español and P. B. Warren, Europhys. Lett., 1995, 30, 191 CrossRef.
  55. R. D. Groot and P. B. Warren, J. Chem. Phys., 1997, 107, 4423–4435 CrossRef CAS PubMed.
  56. I. Pagonabarraga, M. H. J. Hagen and D. Frenkel, Europhys. Lett., 1998, 42, 377–382 CrossRef CAS.
  57. G. Besold, I. Vattulainen, M. Karttunen and J. M. Polson, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 62, R7611–R7614 CrossRef CAS.
  58. I. Vattulainen, M. Karttunen, G. Besold and J. M. Polson, J. Chem. Phys., 2002, 116, 3967–3979 CrossRef CAS PubMed.
  59. T. Soddemann, B. Dünweg and K. Kremer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 046702 CrossRef.
  60. W. H. Press, B. P. Flannery, S. A. Teukolsky and W. T. Vetterling, Numerical Recipes: The art of scientific computing, Cambridge University Press, NY, New-York, 2007 Search PubMed.
  61. F. Goujon, A. Ghoufi, P. Malfreyt and D. J. Tildesley, Soft Matter, 2012, 8, 4635–4644 RSC.
  62. L. J. Fetters, D. J. Lohse, D. Richter, T. A. Witten and A. Zirkel, Macromolecules, 1994, 27, 4639–4647 CrossRef CAS.
  63. M. Rubinstein and R. H. Colby, Polymer Physics, Oxford University Press, Oxford, 2003 Search PubMed.
  64. H. Wang, C. Junghans and K. Kremer, Eur. Phys. J. E, 2009, 28, 221–229 CrossRef PubMed.
  65. S. K. Sukumaran, G. S. Grest, K. Kremer and R. Everaers, J. Polym. Sci., Polym. Phys. Ed., 2005, 43, 917–933 CrossRef CAS.
  66. H. J. C. Berendsen, J. P. M. Postma, W. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690 CrossRef CAS PubMed.
  67. H. M. James and E. Guth, J. Chem. Phys., 1947, 15, 669–684 CrossRef CAS PubMed.
  68. W. Kuhn, J. Polym. Sci., Polym. Phys. Ed., 1946, 1, 380–388 CAS.
  69. L. R. G. Treolar, The Physics of Rubber Elasticity, Clarendon Press, Oxford, 1975 Search PubMed.

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
Click here to see how this site uses Cookies. View our privacy policy here.