V. M.
Nazarychev
a,
A. V.
Lyulin
b,
S. V.
Larin
a,
A. A.
Gurtovenko
ac,
J. M.
Kenny
ad and
S. V.
Lyulin
*ac
aInstitute of Macromolecular Compounds, Russian Academy of Sciences, Bol'shoi pr. 31 (V.O.), St. Petersburg, 199004, Russia. E-mail: s.v.lyulin@gmail.com
bTheory of Polymers and Soft Matter Group, Technische Universiteit Eindhoven, PO Box 513, 5600 MB Eindhoven, The Netherlands
cFaculty of Physics, St. Petersburg State University, Ul'yanovskaya str. 1, Petrodvorets, St. Petersburg, 198504, Russia
dMaterials Science and Technology Centre, University of Perugia, Loc. Pentima, 4, 05100 Terni, Italy
First published on 16th March 2016
The results of atomistic molecular-dynamics simulations of mechanical properties of heterocyclic polymer subjected to uniaxial deformation are reported. A new amorphous thermoplastic polyimide R–BAPO with a repeat unit consisting of dianhydride 1,3-bis-(3′,4,-dicarboxyphenoxy)diphenyl (dianhydride R) and diamine 4,4′-bis-(4′′-aminophenoxy)diphenyloxide (diamine BAPO) was chosen for the simulations. Our primary goal was to establish the impact of various factors (sample preparation method, molecular mass, and cooling and deformation rates) on the elasticity modulus. In particular, we found that the elasticity modulus was only slightly affected by the degree of equilibration, the molecular mass and the size of the simulation box. This is most likely due to the fact that the main contribution to the elasticity modulus is from processes on scales smaller than the entanglement length. Essentially, our simulations reproduce the logarithmic dependence of the elasticity modulus on cooling and deformation rates, which is normally observed in experiments. With the use of the temperature dependence analysis of the elasticity modulus we determined the flow temperature of R–BAPO to be 580 K in line with the experimental data available. Furthermore, we found that the application of high external pressure to the polymer sample during uniaxial deformation can improve the mechanical properties of the polyimide. Overall, the results of our simulations clearly demonstrate that atomistic molecular-dynamics simulations represent a powerful and accurate tool for studying the mechanical properties of heterocyclic polymers and can therefore be useful for the virtual design of new materials, thereby supporting cost-effective synthesis and experimental research.
However, due to the complex chemical structure of the diamine and dianhydride fragments of PI repeat units, the relationship between the chemical composition of the thermoplastic polyimides and their mechanical properties remains unclear.13 Furthermore, due to their complex chemical structure, the thermoplastic polyimides may show relatively low local translational and rotational mobilities of polymer chain fragments, in contrast to those observed for commodity polymers, such as polyethylene or polystyrene. Overall, these features can lead to specific changes in the relaxation patterns of different fragments of repeat units of thermoplastic polyimides subjected to an external deformation field. They can also significantly affect the correlation between their mechanical properties and temperature, pressure and other external factors, as compared to the corresponding correlations seen for polymers of simpler chemical structure.
Therefore, the synthesis of new polyimides requires large-scale experimental research, including an analysis of the effects of temperature and external pressure. Unfortunately, such experimental studies are quite resource-intensive.1,2 That is why computer simulation techniques used along with fully-atomistic computational models, are regarded as efficient and much cheaper methods for the comparative study of mechanical properties of new thermoplastic polyimides.
As far as the computer simulations of mechanical properties of polymers are concerned, it is still not clear how to correctly account for the influence of the thermal history (cooling rate14–16) and the load application rate (deformation rate17–29). Our study will focus exactly on these two issues, namely the sensitivity of the elasticity modulus to the cooling and deformation rates employed in simulations. This will allow us to determine the optimal values of the cooling and deformation rates to be used in simulations in order to appropriately model the thermomechanics of thermoplastic polyimides.
The next section of the paper focuses on the polymer materials in question and the simulation details, as well as on the uniaxial deformation strategy. The effects of the degree of equilibration, cooling and deformation rates, the molecular mass and the size of a periodic simulation cell on the elasticity modulus are studied in the subsequent sections of the paper. Finally, the influence of the temperature and external pressure on mechanical properties of the thermoplastic polyimides is reported and discussed.
Comparing three previously considered thermoplastic polyimides (R–BAPS, R–BAPB and R–BAPO),33–36 the choice of PI R–BAPO for studying the mechanical properties can mainly be justified by the following two reasons. First, in contrast to the partially crystallizing PI R–BAPB, the polyimide R–BAPO is fully amorphous, which allows us to rule out completely the effects related to the crystallization processes and investigate the mechanical properties of a fully amorphous polymer sample. Second, PI R–BAPO is characterized by a negligible anisotropy in the distribution of partial charges. This is in drastic contrast to PI R–BAPS, whose repeat unit contains a highly polar sulfone group that creates additional computational issues.
The polymer samples were then cooled down to temperature of 290 K. In order to examine the influence of the simulated cooling rate γc on the mechanical properties, a cooling procedure was applied with five different values of γc, spanning the range from 1.5 × 1010 K min−1 to 1.5 × 1014 K min−1.
In order to reduce the required simulation time, we do not take into account the electrostatic interactions (EI)35 in the present study. Such an assumption is justified due to the rather small contribution of EI energy of PI R–BAPO to the total potential energy (∼25%),36 electrostatic interactions do modify the mechanical properties of PI R–BAPO but do not change the qualitative conclusions of the present study. Similar simulations without explicit EI are widely used in studies that dealt with polyimide mechanical properties;38–40 the possible influence of EI will be discussed in our future publication.
After the cooling step, each of the prepared 11 samples has been deformed with a stretching rate γd that was varied in the range from 10−1 nm ps−1 to 10−5 nm ps−1. For a periodic simulation cell of length of ∼5.8 nm, the corresponding deformation rate amounted to 1.8 × (106–1010) s−1. The simulations have been performed using the molecular-dynamics (MD) package Gromacs 4.5.641,42 and the force field Gromos53a5.43 All simulations required ∼1.6 million processor hours and were performed using 64 cores on Lomonosov (Intel Xeon X5570 CPUs) and Chebyshev (Intel Xeon E5472 CPUs) supercomputers at Moscow State University.
The uniaxial deformation changes the periodic cell size at a constant rate, along the positive direction of one of the reference axes – X, Y or Z –44 so that the isotropic Berendsen barostat46 with a time constant τp = 0.5 ps was replaced with the anisotropic Berendsen barostat with τp = 1 ps. In the direction of the applied deformation the samples remained incompressible, i.e. the compressibility of the system in this direction was set to zero. In the transverse direction the system compressibility was set to 4.5 × 10−10 Pa−1. Therefore, upon stretching a simulation cell elongates in the direction of deformation and compresses in the directions transversal to deformation in response to the external pressure (1 bar), see Fig. 2.
Fig. 2 PI R–BAPO sample configuration before (left) and after (right) uniaxial deformation along the X axis. The periodic cell is highlighted in blue. |
During deformation the values of the pressure tensor Pi, i = x, y, z, and the simulation cell size Li in the stretching direction were saved each 1 ps. The characteristics obtained were converted to the dependence of the stress σ on the relative strain ε as47
σ = −Pi, |
ε = (Li − L0i)/L0i, | (1) |
σ = Eε. | (2) |
Elastic properties can also be analyzed with the transversal dimensions fixed during the deformation.48 The value of the modulus determined in this way combines the Young's modulus E and the polymer bulk modulus, and differs from the elasticity modulus calculated using the present approach.
The yield peak σy (the onset of plastic deformation in the material) and the strain-hardening modulus Gh (the slope of the linear dependence σ(ε) in the post-yield deformation regime) were calculated as49–52
σ = σy + Gh(λ2 − λ−1). | (3) |
Such behaviour of Young's modulus dependence on the degree of equilibration can be due to the preliminary annealing, where some equilibration of the samples in the glassy state has been done. This equilibration led to the rather weak dependence on the sample thermal history. Thus it might be sufficient for the accurate calculation of the mechanical properties in the present molecular-dynamics simulations.
Our findings also show that the modulus of elasticity for well-equilibrated polymer samples (Fig. S2 and Table S1, ESI†) hardly depends on their molecular mass (Fig. S3 and S4, ESI†). The dependence of mechanical properties on the size of a periodic simulation cell was examined as well, and no sufficient influence was observed (Table S2 and Fig. S5, ESI†). Such insensitivity of the elastic modulus is due to the fact that the modulus is mainly determined by the processes on the scales that do not exceed the entanglement length.
After switching off the LINCS algorithm, the samples obtained by cooling with five different rates γc from 1.5 × 1010 K min−1 to 1.5 × 1014 K min−1 were deformed at a stretching rate γd = 1.8 × 108 s−1, see Fig. 3a. The results show that the extracted values of the modulus E depend logarithmically on the cooling rate γc, see Fig. 3b.
Fig. 3 (a) Stress–strain dependences for different values of γc and (b) cooling-rate dependence of the elastic modulus in simulations; the logarithmic fit is shown by the dotted line. |
The logarithmic relationship between the elasticity modulus and the cooling rate is known from experimental studies on polyheteroarylenes.53 A similar dependence was also reported for a simple two-dimensional model of an amorphous glass.14 To the best of our knowledge, such a logarithmic dependence has never been observed in atomistic computer simulations of thermoplastic polyimides.
Such mechanical behavior of the PI samples with different thermal history is likely to be determined by both differences in the density at room temperature,49 see Fig. 4a, and different values of the typical intramolecular energy of excluded-volume interactions, see Fig. 4b.
Fig. 4 (a) Strain dependence of the PI density and (b) energy of excluded-volume interactions (EL-J) at various cooling rates γc. |
As the cooling rate γc decreases, the density of a sample in the glassy state systematically increases, leading to an increase in the number of contacts between polymer chains and a drop in the free volume. In the strain-hardening regime the dependences ρ(ε) that are calculated for various values of γc reach some constant density, see Fig. 4a. A similar behavior is observed for the strain dependence of the energy of excluded volume interactions, see Fig. 4b, where the saturation becomes clearly visible for deformations that exceed 20%. This is a manifestation of mechanical rejuvenation or (at least in part) of the removal of the thermal history upon deformation.54,55
The results obtained for PI R–BAPO confirm the logarithmic dependence E(γd) observed previously in coarse-grained simulations of polyethylene and polypropylene21,56 and also in experimental studies of various polymers and composite polymer-based materials.57,58
It is noteworthy that for the samples with very high deformation rates that exceed γd = 1.8 × 1010 s−1 (∼100 m s−1) and approach acoustic velocities, the dependence E(γd) clearly deviates from the logarithmic behavior, see Fig. 5b. Coarse-grained simulations of polyethylene demonstrated that the elasticity modulus also deviated from the logarithmic dependence21 at even higher rates (as high as 1011 s−1).
Some possible explanations for the observed non-logarithmic dependence of the elasticity modulus on the deformation rate can be linked to the presence of additional sub-Tg relaxation processes that are activated in polymer glasses under fast deformation. At high stretching rates the polymer internal structure undergoes substantial changes that are manifested in a rather rapid density drop, even in the linear viscoelastic regime (relative deformation of a few percent), see Fig. 6.
Thus, we can conclude that a correct choice of deformation-rate values in atomistic computer simulations is important for determining the mechanical properties of thermoplastic polyimides and for a subsequent comparison with experimental data. In the following the values of the deformation were taken below 1.8 × 109 s−1, i.e. well within the logarithmic regime of E(γd) dependence.
Fig. 7 (a) The stress in PI R–BAPO samples as a function of the engineering strain ε at different temperatures. The inset shows the linear viscoelastic regime. (b) The stress in PI R–BAPO samples as a function of the strain εt at different temperatures. Solid black lines indicate the linear fits in the strain-hardening regime as highlighted by the vertical dotted lines, see eqn (3). |
The initial part of the σ(ε) dependence at different temperatures was used for calculating the temperature dependence E(T) of the elasticity modulus, see Fig. 8.
Fig. 8 The elasticity modulus as a function of temperature. A dashed grey line corresponds to the linear fit. Dashed arrows indicate Tbulk (E equals zero at Tbulk) and Tg of R–BAPO PI as estimated previously.36 |
The analysis of the data in Fig. 8 shows that the modulus drops almost linearly with temperature. Similar dependences E(T) were also found in other simulation studies,21,59 and they were also qualitatively confirmed by experimental data.60,61
When temperature approaches Tg from the low-temperature region, the value of E decreases by ∼45% as compared to the E value at room temperature. Extrapolation of the dependence E(T) to zero gives the flow temperature of Tbulk ∼ 580 K; this corresponds to the temperature at which the highly elastic response of the material to the applied load changes to a viscous flow. One can therefore conclude that the temperature T = 600 K that was selected for the equilibration stage corresponds to highly-mobile PI melt conditions.
The temperature range of the viscoelasticity for the majority of polymers is determined by their chemical structure. For thermoplastic PIs the highly-elastic state can be maintained in a rather broad temperature range from room temperature to temperatures that exceed by ∼100 K the glass transition temperature Tg = 470 K.36
We also analyzed the density–strain dependences at different temperatures, see Fig. 9.
During deformation, the density decreases until a T-dependent constant value is reached at different strain values. The combined influence of temperature and deformation rate makes the sample plastic at lower values of ε when temperature increases. The state of viscous flow defines the material reaction to external pressure, and the polymer liquid becomes practically incompressible under such conditions. Note that we did not observe any significant increase in the density and the magnitude of (negative) excluded volume (LJ). The deformation above Tg has to be done with caution, due to the large fluctuations of these characteristics. To obtain reliable results, the number of simulated samples should be large enough to provide good statistics. The results obtained are in line with the simulation findings for mechanical properties of commodity polymers.44,59
With temperature the yield peak and the strain-hardening modulus decrease,62 see Fig. 10.
Fig. 10 Temperature dependence of the strain-hardening modulus. The inset shows the temperature dependence of the yield peak. The error bars are smaller than the symbol size. |
The analysis of the results presented in Fig. 10 shows that both the yield peak and the strain-hardening modulus decrease almost linearly as temperature increases. In the case of the yield peak dependence such behavior is in good agreement with other computer simulations.44 The theory of rubbery elasticity predicts that the strain-hardening modulus (Gh) should increase with temperature,
Gh = kBTρ/Me | (4) |
The modulus increases linearly with external pressure P from 2.7 GPa (P = 0) to ∼3.7 GPa (P = 500 MPa), see Fig. 12.
The results obtained are qualitatively supported by experimental data for commodity polymers.63,64 To the best of our knowledge, the influence of the external pressure on mechanical properties of polyimides has not been reported yet.
The values of the yield peak σy and the strain-hardening modulus Gh also depend on the amount of the applied external pressure, see Fig. 13. With the increase of pressure both σy and Gh linearly increase. This correlation was established for the first time in simulations of atactic polystyrene under external pressure.49 Our findings in Fig. 13 are qualitatively supported by experimental data for other polymers.28,63,64 It is known that the ratio Gh/σy may act as an empirical relation that describes the strength of the polymer material49,50 (the so-called Considere criterion), which is often related to the neck formation as a result of the deformation when Gh/σy < 1/3. As Gh/σy exceeds 1/3, the material strength is supposed to increase. The dependence Gh/σy on external pressure at T = 290 K clearly distinguishes two regimes: Gh/σy is approximately equal to 1/3 under an external pressure of ≤200 MPa, and is above 1/3 under an external pressure of >200 MPa. Such behavior is likely to be associated with the substantial increase in the material strength under certain values of external pressure.
Summarizing, the application of external pressure has a significant impact on practically all investigated mechanical properties, leading to a considerable increase in the elasticity modulus, the yield peak and the strain-hardening modulus. The strength of thermoplastic PI R–BAPO can be controlled by changing the external pressure; this property can be used to create new high-performance materials.
The dependence of the elasticity modulus on the deformation in the acoustic regime (∼1010 s−1) deviates from the logarithmic dependence observed at much lower deformation rates (in both simulations and experiments). These findings most likely can be related to the thermally-activated nature of the polymer glass deformation due to the possible additional contributions of fast sub-Tg relaxation processes that are activated at high deformation rates.
Extrapolation of the T-dependence of the elasticity modulus to the zero modulus value was used to determine the temperature of transition from the highly elastic state to the viscous-flow state, Tbulk ∼ 580 K. Simulations showed the linear dependence of the yield peak and the strain hardening modulus on external pressure in agreement with previously obtained data.49 We also found an almost linear dependence of these characteristics on temperature. Our computational results also showed an increase in higher strength performance of R–BAPO with temperature (for temperatures below the glass transition temperature). Applying high external pressure (200 MPa and above) at a temperature of 290 K may lead to a further improvement of PI strength.
Overall, our study clearly demonstrates the ability of the state-of-the-art atomic-scale MD simulations to appropriately describe the mechanical properties of heterocyclic polymers. The developed protocols will be used in the nearest future for a comparative study of mechanical properties of composite materials based on these thermoplastic polyimides, which are reinforced with carbon nanofillers of various shapes and sizes.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6sm00230g |
This journal is © The Royal Society of Chemistry 2016 |