Journal of Non-Crystalline Solids 352 (2006) 5008–5012 www.elsevier.com/locate/jnoncrysol
Simulation of polymer glasses: From segmental dynamics to bulk mechanics Alexey V. Lyulin *, M.A.J. Michels Group Polymer Physics, Eindhoven Polymer Laboratories and Dutch Polymer Institute, Technische Universiteit Eindhoven, P.O. Box 513 5600 MB Eindhoven, The Netherlands Available online 21 August 2006
Abstract Molecular dynamics (MD) simulations of the low-molecular-weight glass-former, isopropylbenzene (iPB), atactic-polystyrene (PS) and bisphenol A polycarbonate (PC) reveal the polymer-specific effects of the glassy segmental orientational dynamics. We also simulate the mechanical responce of these two typical amorphous polymers, brittle PS and tough PC, under uniaxial extension. The Young moduli, yield stresses, strain-softening and strain-hardening phenomena of the simulated polymers are numerically close to the experimental data. The local mobility of the segments in the deformation direction is increased drastically (but differently) for both polymers beyond the yield point. In the non-deformed case three relaxational processes have been found for all simulated glass-formers: high-frequency transient ballistic process, b process corresponding to the motions within the cage, and final cage relaxation. For all three glass-formers the slowing down of the cage relaxation (a relaxation) follows mode-coupling theory above Tg, with non-universal exponents, and continues below Tg as an activated relaxational process. For iPB, PS and PC at high temperature the a process is merged with the b process. The b process can be described by an activation law both above and below Tg. 2006 Elsevier B.V. All rights reserved. PACS: 62.20.Fe; 83.10.Mj; 64.70.Pf Keywords: Glass transition; Mechanical, stress relaxation; Molecular dynamics; Polymers and organics
1. Introduction Understanding the mechanisms responsible for the tremendous slowing down of the local mobility when approaching the glass transition is one of the most important challenges in the modern soft condensed matter physics both for low-molecular-weight (LMW) and polymeric materials. Our recent molecular dynamics (MD) simulation studies [1–3] on the atomistic model of an atactic-polystyrene (PS) melt reveal the details of segmental polymer dynamics in the glassy state. The simulations show that above Tg the characteristic translational diffusion time for this polymer follows the same universal power law as the characteristic time of dynamic scattering, in good agree*
Corresponding author. Tel.: +91 44 52197750; fax: +91 44 52197792. E-mail address:
[email protected] (A.V. Lyulin).
0022-3093/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jnoncrysol.2006.02.156
ment with the predictions of the mode-coupling theory (MCT) [4]. At the same time, the character of mobility is qualitatively changed below Tg: relaxational segmental dynamics is described by one and the same activation law. This is completely different from the divergence of relaxation times predicted by the MCT. In order to investigate the possible polymer-connectivity-specific effects of the glassy dynamics we have simulated the orientational relaxation of the LMW analogue of PS, isopropylbenzene (iPB), in a melt at an isotropic atmospheric pressure and in a wide temperature range, both above and below the corresponding Tg. The characteristic relaxational times of iPB are within the time window usually covered by the MD simulation, even at very low temperature. This fact allows us to separate different relaxational processes with rather high precision. We also report here the results of MD simulations for a typical
A.V. Lyulin, M.A.J. Michels / Journal of Non-Crystalline Solids 352 (2006) 5008–5012
brittle polymer, atactic-polystyrene, and a typical ductile polymer, bisphenol A polycarbonate. We have realistically reproduced the experimentally observed [5] strain-softening and strain-hardening phenomena for these polymers, in spite of the inevitable difference between experimental and simulated deformation rates. The MD simulations of the orientational segmental mobility of these polymers in the isotropic state have been performed as well, and the temperature-dependent relaxation times have been analyzed. 2. Algorithm and simulation details United-atom models are used in the present study; the models, the force fields and details of the cooling procedure are described in detail by Lyulin et al. [1–3]. In short, the NPT MD simulation has been performed for a model PS melt consisting of 1–8 polymer chains of Np = 80–160 monomers each. sixty four chains of Np = 10 monomers (about one entanglement length, 64 · 10 system) and eight chains of Np = 80 monomers (about eight entanglement lengths, 8 · 80 system) have been used for the PC melt. For the present calculations well-equilibrated initial structures have been taken from our previous studies. In order to check the quality of the equilibration we have also started to use the new effective end-bridging Monte Carlo algorithms [6], this study is in progress now. The simulated values of Tg for both polymers (375 K for PS and 433 K for PC) were reported in our previous study [7]. Uniaxial deformation with constant deformation velocity L_ ¼ ˚ /ps was applied along one (X, Y or Z) of the axes 0:005 A to five independent sets of the relaxed isotropic PS and PC samples. All simulations in the non-deformed (isotropic) state have been performed by using the NVE ensemble. A low-molecular-weight analogue of polystyrene, isopropylbenzene (iPB), C9H12, has been also simulated for N = 80 and N = 640 molecules, in a wide range of temperatures both above and below its Tg (the glass transition was simulated by the step-wise cooling from the initial liquid
5009
state at T = 360 K with steps of 30 K, and equilibration of about 4 ns between two steps). The simulated value of Tg was 175 ± 15 K. 3. Results: Uniaxial deformation The uniaxial deformation was applied to the relaxed polymer samples of PS and PC at room temperature (T = 300 K), well below their corresponding Tg. After the application of the deformation, the initial elastic regime for both polymers is clearly seen for extensions up to 3%. The elastic regime is followed by the onset of plastic flow in the yield point, Fig. 1. At T = 300 K the simulated Young moduli and values of the yield stress and yield strain are comparable to the known experimental values at room temperature, in spite of the enormous difference in deformation rate. After the yield point pronounced strain-softening is clearly seen in Fig. 1 for simulated PS. Such a softening is not present for PC. For both polymers strain hardening starts to develop at larger extensions (above deformations of about 20%) with the value of the strain hardening modulus GR = 12.0 ± 0.5 MPa for PS (almost twice lower than the simulated value for PC). These results are in fair agreement with available experimental data [8]. The segmental translational mobility (mean-square translational displacement, MSTD) has been studied for both polymers both in the isotropic state and under the uniaxial deformation, Fig. 2. The trivial contribution to the MSTD due to the convective motion of the monomers under deformation was taken into account and is not shown in Fig. 2. For both polymers it is clearly seen that the deformation does not influence significantly the translational mobility below the yield point (up to t 200 ps). Changes occur in the post-yield response: for the PS glass the uniaxial mechanical deformation leads to an acceleration of the parallel diffusion by more than one order of magnitude. For PC the observed acceleration is weaker. Is it possible to relate these differences to the orientational segmental mobility of these two mechanically different
Fig. 1. Stress-strain behavior of a non-entangled PS melt (left) and two specimens of PC of different molecular weight (right) at T = 300 K. Strainsoftening is clearly seen for a PS sample, and is clearly absent for PC. The straight lines are the least-square fits of the data in the strain-hardening regime. The strain-hardening moduli were calculated from these fits for both polymers. Results are averaged over all independent runs. Typical fluctuations are 10 Mpa.
5010
A.V. Lyulin, M.A.J. Michels / Journal of Non-Crystalline Solids 352 (2006) 5008–5012
polymers already in the isotropic, non-deformed state? If so, what are the polymer specific effects? In particular, what is the effect of the hindered motions of the side phenyl groups for PS local dynamics? The remaining part of the paper is devoted to the detailed study of the local orientational mobility in undeformed iPB, PS and PC. 4. Results: Local orientational mobility
2
P 2 ðtÞ ¼ h3=2ðbð0ÞbðtÞÞ 1=2i;
where b is the unit vector directed along the main-chain bond or along the phenyl side-group (in PS or iPB) and brackets denote the averaging over positions and over time. Kohlrausch–Williams–Watts (KWW) [9] stretched exponential functions have been used to fit the final parts of the ACFs at different temperatures, P 2 ðtÞ ¼ expððt=s2 Þb Þ;
To have insight into possible common mechanisms underlying the segmental mobility of these mechanically extremely different polymers, and to reveal the possible polymer-specific effects, the local orientational mobility in undeformed iPB, PS and PC has been studied with the help of Legendre polynomials of the second-order (autocorrelation functions, ACFs)
ð1Þ
ð2Þ
where s2 is the characteristic orientational relaxation time and b is the parameter effectively taking into account the non-exponential nature of the relaxation process. Each ACF is also analyzed using the CONTIN [10] method Z þ1 P 2 ðtÞ ¼ F ðln sÞ expðt=sÞdðln sÞ; ð3Þ 1
where F(ln s) is a normalized distribution function of relaxation times. This function is shown for iPB in Fig. 3. Three relaxational processes can be identified at each temperature. The very fast transient process (times less than 1 ps) is due to the ballistic motion of monomers, is purely determined by the thermal velocity of the united atoms at each prescribed temperature, and has very weak temperature dependence (Fig. 3, right). On the other side of the spectrum of the relaxation times is the so-called a process. The temperature dependence of the a-relaxation times above Tg is clearly non-Arrhenius, Fig. 3, and is fitted by the Vogel–Fulcher–Tamman (VFT) empirical equation, s1 expðEa =k B ðT T 0 ÞÞ; Fig. 2. The fluctuating stochastic contributions to the parallel MSTD (symbols, averaged over backbone monomers for PS and over all monomers for PC) together with the corresponding dependencies in the isotropic case (lines), 8 chains by 80 monomers, T = 300 K. For PS the onset of the Rouse regime in the direction of the deformation starts much earlier as compared to the isotropic case. For PC the translational diffusion just after the yield point is changed less by the deformation. The error bars are smaller than the symbols size.
ð4Þ
with Ea = 1 ± 0.1 kJ/mol and T0 = 150 ± 10 K < Tg = 175 K. The temperature dependence of the relaxation times of this process is changed to an activated one below Tg, and a simple activation law with Ea = 2.7 ± 0.1 kJ/mol is used to fit the simulation data. This change in the character of the relaxational behavior from the VFT-like behavior above Tg to the activated behavior below Tg is qualitatively
Fig. 3. Normalized distribution function Eq. (3) of the relaxation times for the P2(t) ACFs for the backbone bond of iPB at different temperatures (left). Three relaxation processes are clearly seen at all temperatures. The temperature dependence of the characteristic relaxation times is also shown (right), for the backbone (closed symbols) and the side-group vectors (open symbols): times for the b and transient processes are calculated from the positions of the corresponding maxima in the CONTIN distribution functions; times for the a and activated processes are calculated from CONTIN (inverse triangles) and from a KWW fit (squares). The VFT fit, Eq. (4), for the a process above Tg is also shown. Straight dashed lines represent an activation law fit for the b- and activated relaxation processes. For all data sets the error bars are smaller than the symbols size.
A.V. Lyulin, M.A.J. Michels / Journal of Non-Crystalline Solids 352 (2006) 5008–5012
similar to that observed for PS in [1–3]. The MCT fit [4], s (T Tc)c, is also used to fit the T-dependence of the a-relaxation times above Tg (the MCT-like fit is very close to the fit produced by Eq. (4) and is not shown in Fig. 3 for clarity) and gives Tc 181 K (above the simulated value of Tg, qualitatively in agreement with the predictions of MCT). The value of the critical exponent for iPB is c 1 which is significantly lower than the critical exponent for PS (c 2.8) calculated in [7]. This difference reflects the well-known molecule-specific aspects of MCT, even for chemically related structures. The temperature dependence of the characteristic times of the b process can also be extracted from the CONTIN analysis. Also this process is described by a simple activation law, both below and above Tg, with Ea = 6.8 ± 0.2 kJ/mol. From Fig. 3 it is seen that at high temperature the a and b processes are merged; they start to split at temperatures above Tg. How does this picture change for specific polymers? The analysis explained above for iPB has been implemented here also for PS and PC, Fig. 4. As in the case of iPB, three relaxational processes can also be identified at each of the
5011
simulated temperatures: the very fast transient process, the b process, and, finally, the a process. At high temperature the b process is merged with the a process for both polymers; they start to split at temperatures well above the corresponding Tg. This splitting starts already at rather high temperature for PS, Fig. 5, left panel. MCT fits for the a process hold for both polymers above Tc > Tg (see Fig. 5), with values of the critical exponents much larger than for iPB (c 2 for PC, c 3 for PS and c 1 for iPB). The VFT fits, Eq. (4), also hold (they are not shown in Fig. 5 for clarity), but with the values of T0 much lower than the corresponding Tg for both polymers. 5. Discussion Of course at temperatures below Tg both the low-molecular and polymer systems are in the non-ergodic state and the calculated relaxation times produced by the fit (2) are outside the simulated time window. Nevertheless, simulations for 80 ns allow for almost 30% of the relaxation of the P2 autocorrelation function for PS at T = 370 K. To
Fig. 4. Normalized distribution functions of the relaxation times for the P2(t) ACF for the backbone bond of PC (left) and PS (right) at different temperatures, calculated using the CONTIN method. Three relaxation processes are clearly seen for both polymers at temperatures close to the corresponding Tg. At high temperature (for example, at T = 600 K for PS) because of the merging of the a and b processes only two processes can be unambiguously identified. For all data sets the error bars are smaller than the symbols size.
Fig. 5. Temperature dependence of three relaxation processes for the backbone (closed symbols) and the side-group vectors (open symbols for PS) of PS (left, 1 · 80 melt) and PC (all bond vectors are averaged) (right, 64 · 10 melt). Times for the a process and its continuation below Tg are calculated from a KWW fit and also using CONTIN (grey circles). The MCT fit for the a process above Tsg is shown for both polymers. Straight dashed lines represent an activation law fit for the b- and activated relaxation processes. For all data sets the error bars are smaller than the symbols size.
5012
A.V. Lyulin, M.A.J. Michels / Journal of Non-Crystalline Solids 352 (2006) 5008–5012
check the accuracy of the fit (2) two different initial configurations produced by the cooling from the high temperature melt for PS sample have been used for the production runs at T = 270 K Tg. The calculated ACFs were not that different and the corresponding relaxation times were averaged. The CONTIN analysis used in the present paper for both polymers is also quite precise because (even for temperatures slightly below Tg) the simulated window is close to the position of the maximum in the time dependence of the spectral density at this temperature. This fact allows us to speak about the reliability of our results. For even lower temperatures the relaxations times are 2 orders of magnitude above the simulated window. Still, even at these temperatures there exists some relaxation of the correlation functions. Not surprisingly the absolute values of the ballistic time for polymers are very close to those for the low-molecularweight compound. As in the case of iPB, their temperature dependence for PS and PC is also very weak (Fig. 5). For both polymers the temperature dependence of the relaxation times of the b process is very close to an activation law, but the values of the activation energy are much larger as compared to iPB (Ea = 20 ± 2 kJ/mol for PS and Ea = 28 ± 2 kJ/mol for PC, as compared to Ea = 6.8 ± 0.2 kJ/mol for iPB). As in the case of iPB the a relaxation is continued below Tg but the character of this relaxation is changed from the MCT-like behavior to an activated behavior. Again, qualitatively this picture is similar to that for the low-molecularweight glass, but the activation energy is one order of magnitude larger: Ea = 24 ± 2 kJ/mol for PS, Ea = 23 ± 2 kJ/mol for PC as compared to Ea = 2.7 ± 0.1 kJ/mol for iPB. 6. Conclusion The local orientational segmental mobility for the lowmolecular-weight glass-former, isopropylbenzene, atacticpolystyrene and bisphenol A polycarbonate has been simulated in the non-deformed isotropic state. The molecular dynamics simulations for the uniaxial deformation of the chemically realistic models of these two polymers clearly
have been performed as well, reproducing the experimentally observed strain softening and strain-hardening, in spite of the difference between computational and experimental deformation rates. For all simulated non-deformed compounds the analysis of the distribution functions of the orientational relaxation times for the second-order ACFs (Legendre polynomials) reveals the existence of three different relaxation processes: an initial ballistic process at very short times, followed by the b relaxation (motions within the cage) and, finally, the relaxation of the cage itself. This final relaxation, which is a collective a-relaxation process above Tg, is changed into an activated process below Tg, with very different activation energies for the polymerand low-molecular-weight glasses. The b relaxation is an activated process which continues below Tg, and merges with the a process at high temperature. This picture is universal for both main-chain and side-group (phenyl rings for PS) orientational motions in the isotropic case. The predictions of the MCT are non-universal, and depend on the specific chemical structure; in particular value of the critical exponent is c 1 for iPB and increases significantly for polymers. The effect of the uniaxial deformation on the local orientational mobility of the two polymers is the subject of our future investigation.
References [1] A.V. Lyulin, M.A.J. Michels, Macromolecules 35 (2002) 1463. [2] A.V. Lyulin, N.K. Balabaev, M.A.J. Michels, Macromolecules 35 (2002) 9595. [3] A.V. Lyulin, N.K. Balabaev, M.A.J. Michels, Macromolecules 36 (2003) 8574. [4] W. Go¨tze, J. Phys. Condens. Matter. 11 (1999) A1. [5] E.T.J. Klompen, T.A.P. Engels, L.E. Govaert, H.E.H. Meijer, Macromolecules 38 (2005) 6997. [6] M. Doxastakis, V.G. Mavrantzas, D.N. Theodorou, J. Chem. Phys. 115 (2001) 11352. [7] A.V. Lyulin, B. Vorselaars, M.A. Mazo, N.K. Balabaev, M.A.J. Michels, Europhys. Lett. 71 (2005) 618. [8] H.G.H. van Melick, L.E. Govaert, H.E.H. Meijer, Polymer 44 (2003) 2493. [9] G. Williams, D.C. Watts, Trans. Faraday Soc. 66 (1970) 80. [10] S. Provencher, Comput. Phys. Commun. 27 (1982) 213.