- Email: [email protected]
Molecular dynamics calculation to clarify the relationship between structure and mechanical properties of polymer crystals: the case of orthorhombic polyethylene
Computational and Theoretical Polymer Science 11 (2001) 357±374
www.elsevier.com/locate/ctps
Molecular dynamics calculation to clarify the relationship between structure and mechanical properties of polymer crystals: the case of orthorhombic polyethylene K. Tashiro* Department of Macromolecular Science, Graduate School of Science, Osaka University, 1-1 Machikaneyama-cho, Toyonaka, Osaka 560-0043, Japan Received 21 September 2000; revised 20 December 2000; accepted 26 December 2000
Abstract The molecular dynamics (MD) technique was used to calculate the temperature dependence of the structure, molecular motion, and mechanical property of the orthorhombic polyethylene (PE) crystal. The potential functional parameters reported by Karasawa et al. (J Phys Chem, 95 (1991) 2260) were re®ned further so that the vibrational frequencies of infrared and Raman bands, measured by us at ultra-low temperatures for the normal and fully deuterated PE, could be reproduced well. The ¯ip-¯op motion around the chain axis and the torsional motion of the skeletal chains were found to start above ca. 350 K and increase the amplitude of these motions progressively. Coupling these two types of chain motion resulted in a steep increase of the thermal vibration parameters or the mean-square-displacements of carbon and hydrogen atoms, corresponding well with the X-ray data. The lattice constants and the related linear thermal expansion coef®cients were also found to be in good agreement with the observed data. The calculated Young's modulus along the chain axis decreased gradually with the increasing temperature: 330 GPa at 0 K to 280 GPa at room temperature. The latter was in good agreement with the value of 280±305 GPa evaluated from the Raman measurement of the longitudinal acoustic mode. Young's modulus was found to relate intimately with the chain contraction caused by the skeletal torsional motion. Only 0.3% contraction of the chain resulted in the reduction of the modulus by ca. 35%. A similar behavior was also seen in the trigonal polyoxymethylene and nylon 6 a forms. q 2001 Elsevier Science Ltd. All rights reserved. Keywords: Molecular dynamics; Polyethylene; Crystal structure; Molecular motion; Young's modulus; Potential functions
1. Introduction In a series of papers [1±3] we have discussed the intimate relationship between the structure and mechanical properties of polymer crystals from the molecular level by way of theoretical calculations and related experiments. At the present stage we may state that, as long as the structure of a polymer substance is known reasonably well, Young's moduli in various directions of the polymer crystal can be predicted with fairly high probability. In these evaluations, however, we have not explicitly taken into consideration the effects of the thermal motion of the chains on the mechanical properties of the polymer crystals. The mechanical properties of polymer crystals may be assumed as the limiting values of the bulk samples. The mechanical properties of the bulk samples are a function of temperature. Therefore it is absolutely and essentially important to evaluate the effect of the thermal motion in the theoretical prediction of the mechanical properties of polymer crystals. Such trials were * Tel.: 181-6-6850-5455; fax: 181-6-6850-5455. E-mail address: [email protected] (K. Tashiro).
made by several researchers from both the experimental [4±8] and theoretical [9±16] points of view. The temperature dependence of the elastic constants of polyethylene (PE) crystal was calculated on the basis of the free energy treatment under the quasi-harmonic approximation [13,14], molecular dynamics (MD) method [12,15], and lattice dynamics [16]. On the experimental side, the temperature dependence of Young's modulus along the chain axis (or the crystallite modulus) was measured for PE by the X-ray diffraction method [4±6]. In this case, however, the situation became rather complicated because of the problem of the so-called stress distribution in the bulk sample, giving some uncertainty to the experimental values obtained by the X-ray diffraction method. The temperature dependence of the vibrational frequency of the longitudinal acoustic mode (LAM) may be another useful method, but requires correction to study the effect of the intermolecular interactions [17±20], the mass effect at the end groups of the chains [21], and so on. In the evaluation of the mechanical properties of polymer crystals as a function of temperature, we also need to consider the temperature dependence of the crystal structure.
1089-3156/01/$ - see front matter q 2001 Elsevier Science Ltd. All rights reserved. PII: S 1089-315 6(01)00005-8
358
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Another motivation of this paper is based on our recent study on the temperature dependence of the crystal structure of orthorhombic PE. In a previous paper [22], we carried out the X-ray diffraction measurements of highly oriented PE samples at various temperatures from 2100 to 1908C and could collect numerous X-ray re¯ections at each temperature by using the imaging plate system. By analyzing the obtained sets of re¯ection data, the temperature dependence of the crystal structure was revealed with high reliability. The lattice constants, particularly the a-axial length showed the de¯ection point around 1108C. The setting angle of the planar-zigzag chain measured from the b-axis and the thermal parameters of the constituent atoms C and H also showed similar de¯ection points in the same temperature region. The atomic thermal parameters or the mean-square amplitudes of the atomic vibrations were calculated by Kitagawa and Miyazawa [23] under the quasi-harmonic approximation, and were found to be in good agreement with the above-mentioned X-ray values in the low-temperature region. In the temperature region higher than 108C, however, a large discrepancy was detected between the observed and the calculated thermal parameters. This discrepancy led us to speculate that the large-amplitude thermal motion of the chains become signi®cant in the temperature region above 108C. However, we could not explain the type of thermal motion that occurs in the hightemperature region. This needs to be investigated from a theoretical point of view. In the theoretical calculation of the temperature dependence of the structure, thermal motion, elastic constants and so on, several methods have been presented as mentioned above. One is to treat the free energy of the crystal as a function of temperature [13,14], where the elastic constants are obtained as second derivatives of the free energy with respect to the stress component; and the other method is to introduce the anharmonic terms into the force constants: f f0 1 f 0 DX where f, f0, and f 0 are the force constant and the corresponding harmonic and anharmonic terms, respectively, and DX is the atomic displacement caused by the external perturbation. Phenomenological treatment of the anharmonic effect is to use the mode-GruÈneisen constants, which are de®ned as g 2 dn=d10 =n where n is the vibrational frequency of the mode and 1 is the strain [11,24]. The MD method was also used for the evaluation of the temperature dependence of the elastic constants [9,10]. The MD calculation on the temperature dependence of the lattice and elastic constants of orthorhombic PE crystal was made by Cagin et al. [12,15]. One serious problem in these trials is the reasonableness of the utilized potential function parameters, which could not be said to be so reasonable as to give a satisfactorily good reproduction of the observed vibrational frequencies of the infrared and Raman spectra. Kobayashi et al. re®ned the force constants of the valence force ®eld type so that the calculated vibrational frequencies could be reproduced as reasonably as possible [25,26]. Karasawa et al. presented a set of potential function
parameters (MCXX) [27] and calculated the vibrational frequencies and elastic constants of orthorhombic PE crystal. Palmo et al. re®ned the molecular mechanics potential function for n-alkane chain base on scaled ab initio force ®elds [28]. In spite of these efforts, however, the discrepancy between the calculated and observed vibrational frequencies has not yet been removed satisfactorily. As reported in the present paper, we re®ned the potential function parameters proposed by Karasawa et al. and applied them to the MD calculation of the orthorhombic PE crystal to evaluate quantitatively the thermal behaviors of lattice and elastic constants, particularly Young's modulus along the chain axis. As mentioned above, the MD calculation has already been made by Cagin et al. [12,15], but the total number of atoms included in the crystal lattice was limited to 360±504. These atom numbers may not be large enough to simulate the thermal behavior of the actual crystal, although they check the effect of the model size on the calculated results. In the present study we built the PE lattice model consisting of ca. 5200 atoms in total and calculated the trajectory (the coordinates and velocities of the constituent atoms and the lattice constants) over the long time span of 100 ps under the three-dimensional periodic boundary conditions. On the basis of the trajectory data obtained at the various temperatures, Young's modulus along the chain axis was calculated and compared with the observed data. It took about 50 days to obtain the trajectory ®le at one temperature point because of the limited power of our computer system. It took us as long as one and a half years, inclusive of the time taken for the data analysis, to obtain the results presented in this paper. A similar calculation was also made for trigonal polyoxymethylene (POM). As the details have already been already reported in Ref. [29], only the results will be compared with those of orthorhombic PE; and some common features on the relationship between the structure and mechanical property of polymer crystals will be extracted. 2. Re®nement of potential function parameters On the basis of the vibrational spectroscopic data, the force constants or the potential function parameters concerning the intra and intermolecular interactions were proposed for PE crystals [25±27]. As already mentioned above, Karasawa et al. proposed a set of potential functions (MCXX), which gave a relatively good agreement between the observed and calculated vibrational frequencies for orthorhombic PE crystal. Recently Sun reported the potential function parameters, COMPASS, which gave a good reproduction of the equations of state for many organic substances [30]. By using their proposed parameters we calculated the vibrational frequencies for orthorhombic PE crystals of both the hydrogenous and the deuterated species and compared the results with the observed values, though the comparison was not satisfactory. We then tried to
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
359
Table 1 Comparison between the calculated and observed values for the various parameters Unit cell parameters
Karasawa et al. [27]
Ê) a(A Ê) b (A Ê) c (A Vibrational spectra (cm 21) CH2
7.202 4.795 2.546 Lattice mode Lattice mode Lattice mode Rocking Bending
CD2
Lattice mode Rocking Bending
Young's modulus (GPa)
a
Karasawa et al. 84 115 132 784 798 1484 1485 ± ± ± ± ±
a-axis b-axis c-axis
This study 6.995 4.733 2.531
Observed (4 K) [65] 7.121 4.851 2.548
This study 97.1 108.2 123.3 718.7 735.2 1494.2 1496.7 89.4 513.9 525.5 1089.0 1093.0
Observed (7 K) 79.1 108.5 137.0 721.8 734.7 1460.8 1474.8 74.5 519.1 529.1 1083.8 1094.5
This study 10.9 9.9 329.5
Observed 10.2 (110) a ± See text
Ref. [47] (Young's modulus along 110 direction).
modify the potential function parameters proposed by Karasawa et al. further so that the vibrational spectroscopic data observed for various kinds of samples such as the hydrogenous PE, the deuterated PE, the random copolymers of H and D ethylene monomeric units, and the PE blends of the pure H and D species could be reproduced consistently and satisfactorily [31,32]. The re®nement of the potential function parameters was made by a trial-anderror method. The calculation of the vibrational frequencies was carried out for the energetically minimized structure models by using the commercial software Cerius 2 (Version 3.8, Molecular Simulations Incorporation) on the basis of the Hessian-matrix method. The set of potential function parameters obtained is listed in Appendix A, where the de®nitions of the potential functions and the corresponding parameters are given concretely. The comparison between the observed and calculated values is shown in Table 1, where the lattice constants, the vibrational frequencies, and the crystallite moduli are listed. We have recently collected the observed vibrational data at 7 K [33]. As seen in Table 1, the values calculated by the potential function parameters listed in Appendix A are in good agreement with the observed data. The detailed discussion on the vibrational data will be reported in a separate paper [33]. The discussion about the calculated Young's modulus will be made in a later section of the present paper.
(Version 3.8, Molecular Simulations Incorporation). In an MD cell all the 72 chains were packed, each of which was constructed by 24 CH2 monomeric units. Therefore the basic MD cell contained 5184 atoms in total. The threedimensional periodic boundary conditions were applied: the 24 CH2 chain segments were connected in®nitely along the c-axis by covalent bonding. Fig. 1(a) shows the thus-constructed starting model. By applying the abovementioned potential function parameters, the atomic charges of 20.0533 and 10.0267 electron units were assumed for all the carbon and hydrogen atoms, respectively, which were calculated by the Gasteiger method [34]. The MD calculation was made for the energetically minimized model under the constant NPT condition by the Hoover±Nose method, where the temperature (T ) was ®xed to the predetermined value (0, 100, 200, 300, 350, 400, and 500 K), and the external pressure was assumed to be 0 atm. The time interval of each step was 1 fs and the total time taken for calculation was 100 ps at each temperature. The trajectory data created in the MD calculations were saved at every 0.1 ps and analyzed by using our home-made programs to obtain the various physical quantities described in the following sections.
4. Results and discussion 4.1. Temperature dependence of lattice constants
3. Molecular dynamics calculations The MD calculation was performed by using the Cerius 2
The lattice constants a, b, and c were evaluated by averaging the trajectory data in the time region reaching the constant energy. The averaged angles a , b , and g were
360
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Fig. 3. Comparison of the temperature dependence of the b-axial length between the observed and calculated values. The horizontal bars in this ®gure indicate the temperature regions where the linear thermal expansion coef®cients were calculated.
Fig. 1. Snapshots of PE crystal model viewed along the chain axis at: (a) 0, and (b) 300 K.
Fig. 2. Comparison of the temperature dependence of the a-axial length between the observed and calculated values. The horizontal bars in this ®gure indicate the temperature regions where the linear thermal expansion coef®cients were calculated.
essentially 908. Figs. 2±4 show the temperature dependencies of the calculated a-, b-, and c-axial lengths. The absolute values of the lattice constants a, b, and c were not suf®cient to reproduce the experimental values accurately. This disagreement could probably be because of the unsatisfactory re®nement of the used potential function parameters, which is still being made to obtain a better agreement between the observed and calculated values. At the same time we cannot ignore the conditions of the experiment. It is a well-known fact that the lattice constants of polymers are dependent sensitively on the sample preparation conditions. Therefore it may be safer to compare the calculated and observed cell constants by taking into account some amount of ¯exibility. The thermal expansion itself is found to show almost the same behavior irrespective of the sample except
Fig. 4. Comparison of the temperature dependence of the c-axial length between the observed and calculated values. The horizontal bars in this ®gure indicate the temperature regions where the linear thermal expansion coef®cients were calculated.
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
361
Fig. 5. Time dependence of the setting angle of an arbitrarily chosen chain calculated at the various temperatures.
for the extremely different sample preparation conditions. Therefore, here we will focus our attention mainly on the thermal expansivity. The a-axis shows the most remarkable temperature dependence, which is not linear over the whole temperature range. The linear thermal expansion coef®cient a a was calculated from the slope of the curve in each temperature region and compared with the observed values as shown in Fig. 2. In Ref. [22] the temperature dependence of a- and b-axes was measured, but the accuracy was not very high. The measurement was made again by increasing the sample-to-camera distance for the same ultra-drawn sample reported in the previous paper. Although the aaxial length itself is a little different, the calculated a a is in good agreement with the observed value in a wide
temperature region. The temperature dependence of the baxis is one digit smaller than that of the a-axis. The linear thermal expansion coef®cient a b is compared with the observed value in Fig. 3. The b-axis is contracted in the high-temperature region, and therefore the a b is negative, which is also reproduced relatively well by the calculation. The c-axis along the chain direction contracts by increasing the temperature as shown in Fig. 4, where the linear thermal expansion coef®cient a c is negative and two digits smaller than a a. The agreement between the observed and calculated values is good, ca. 25 £ 10 26 K 21 at 250 K. When the temperature dependencies of the a-, b- and caxes are investigated carefully, the de¯ection points are found around 350 K although the degree of the de¯ection
362
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Fig. 6. Distribution of the setting angle of the chain extracted from Fig. 5.
is not very high. This might correspond to the actually observed phenomenon: the a-axial length, e.g. increased linearly with temperature but changed its slope gently above the room temperature [22]. The degree of this de¯ection is dependent on the sample. The uniaxially oriented sample with draw ratio of ca. 15 times the original length shows the clearly detectable change around 108C, while this change is small for the ultra-drawn sample with the draw ratio of ca. 150 times. As discussed in Ref. [22], one possible origin of this different behavior may be the difference in the crystallite size between these two different sample types; the ultra-drawn sample should show quite a long and wide crystallite size [5,35]. We may say that the in®nitely large crystal of PE shows a small but signi®cant change in the curve of the a-axial length vs temperature around 108C, and that the difference in the degree of the curve de¯ection might come from the effect of the crystallite surface. For the smaller crystallite, the surface of the lamella gives stress on the structure and molecular motion of the inner part of
the crystallite [36,37]. In fact, several researchers including us clari®ed the change in the unit cell parameters and atomic thermal parameters depending on the samples with a different history [22,38]. By taking it into consideration that the three-dimensional periodic boundary conditions are assumed in the present MD calculations, the small but detectable change in the curves shown in Figs. 2±4 may correspond to the behavior of a very large crystal of orthorhombic PE. However, we need to clarify the concrete origin of these de¯ections in the cell parameters, which will be discussed in the following section. 4.2. Temperature dependence of molecular motion The main purpose of the present paper is to clarify the relationship between the temperature-dependent mechanical property and the conformational disordering of chains caused by thermal motion. As for the molecular motion of PE chains, a detailed discussion was made by Noid et al. in a
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
363
Fig. 7. A series of snapshots of PE crystal taken at 400 K in the time region of 9±12 ps. The molecular chains of red and yellow are named as chains A and B in the text, respectively.
364
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
series of papers on MD calculations [39±44]. Therefore, in the present paper, the description of the thermal motion of PE chains will be made brie¯y so that the information on the essential features of this thermal motion can be helpful for the discussion on the temperature dependence of mechanical property. Fig. 1(b) shows the snapshot taken at an arbitrary time of the trajectory data obtained at 300 K. The orientation of the molecular chains viewed along the chain axis was disordered and at the same time the relative height of the neighboring chains also changed. The setting angle, which is de®ned as an angle between the zigzag plane of the chain and the b-axis, was calculated for the arbitrarily chosen chain in the cell and was plotted against time as shown in Fig. 5. The distribution of the setting angles of the chain is shown in Fig. 6. In the temperature region up to 300 K, the setting angle is almost constant, about 438 in average, but the ¯uctuation of the setting angle from the averaged value increases with the increasing temperature. Around 350 K, as seen in Fig. 5, this ¯uctuation becomes much larger and some chains were found to change their setting angle from 143 to 243 and then to 1438, i.e. the ¯ip-¯op motion was
being observed. This motion becomes more remarkable at 400 K. At 500 K, the setting angles of the chains change almost randomly without any correlation among the neighboring chains. The chain motion at 400 K is now seen in more detail. Figs. 7 and 8 illustrate a series of snapshots picked up at various times. Around 10 and 60 ps, the cooperative rotations of the neighboring chains A (red) and B (yellow) can be observed; the naming of the chains is shown in Fig. 9. For example, around 9.5 ps as shown in Fig. 7, the yellow chain B gradually changed its setting angle from 143 to 2438 (or 11378) by rotating around the chain axis. With some time delay the red chain A started to change the setting angle around 10±10.5 ps. At 11±12 ps these two chains became parallel to each other and at the same time they were parallel to the other neighboring chains. That is, the local packing structure similar to the monoclinic form was generated in the orthorhombic crystal lattice. In Fig. 8 an inverse phenomenon was observed. Around 58 ps the chain B (yellow) changed its setting angle and the red chain A changed the orientation after a delay of about 5 ps. Then the local structure returned to the original orthorhombic type.
Fig. 8. A series of snapshots of PE crystal taken at 400 K in the time region of 50±60 ps. The molecular chains of red and yellow are named as chains A and B in the text, respectively. The small ®gure in the lower right side shows the clusters of monoclinic-type structure.
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
365
Fig. 9. Naming of the PE chains discussed in the text.
How about the correlation between the rotational motion and the translational motion of the chain? Fig. 10 shows the changes in the setting angle and the relative height plotted against time for the chains A and B. For example, at 5 ps, the setting angle of chain B changed steeply and the relative height (the z coordinate) changed correspondingly. Around 55 ps, the relative height changed again after the setting angle returned to its original value. That is, the molecular motion of chain B may be described as a helical motion. The period of this ¯ip-¯op motion of the PE chain is estimated to be about 100 ps at 400 K. However, such a motion was not observed clearly in chain A; around 60 ps the relative height was kept almost unchanged even when the setting angle changed drastically. Therefore, we may say that the rotational motion around the chain axis and the translational motion along the chain axis do not have an intimate correlation with each other. At the same time we should notice that, even when the two neighboring chains A and B changed their setting angles cooperatively, the other chains adjacent to these two chains did not show such ¯ip-¯op motion over the time region of at least 100 ps, as seen in Fig. 11. Besides,
the change in the relative height or the translational motion along the chain axis also did not show any correlation as shown in Fig. 12. This way the ¯ip-¯op motion and the translational motion of the zigzag chain seem to occur with very local and small correlations among the neighboring chains. Local correlation or formation of monoclinic-type clusters can be observed at the various parts of the crystal lattice (see the snapshot of the whole lattice taken at 60 ps in Fig. 9). It should be pointed out here that this monoclinic-type or parallel packing structure could be seen even at a relatively low temperature such as 300 K when the MD calculation was made for much smaller models of ca. 10±20 chains. This is considered to come from overly strong correlation between the chains or from the large effect of the lattice boundary. Besides, the chain packing structure of the monoclinic type is energetically more stable, though only slightly, and competes with the orthorhombic type of packing structure [45]. Therefore once the chain packing changes from the orthorhombic to monoclinic structure, the latter is stabilized and is kept unchanged over a long time. Such a situation was
366
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Fig. 10. Time dependence of the setting angle and relative height of the chains A and B at 400 K.
not seen in the present large model of ca. 5200 atoms. The ¯ip-¯op-type rotation of chains occurred not only in the boundary of the repeating unit cell but also at the central parts of the cell or at an arbitrary position of the cell. At 500 K the ¯uctuation of the setting angles of the chains became drastically large and an almost free rotation was attained at various places. Another important aspect is the relation of this rotational motion with the change in the chain conformation. We calculated the time dependence of the C±C torsional angles of the skeletal chains at various temperatures. The torsional angles ¯uctuated around the center of 1808. As shown in Fig. 13, the amplitude of ¯uctuation became larger as the temperature increased, but the averaged torsional angle was almost 1808 even at 500 K. In Fig. 14 the time dependence of the setting angle, torsional angle, and the relative height of the arbitrarily chosen chain is shown at 350 and 400 K. Any time correlation could not be found for the torsional motion of the chains with the large-amplitude ¯ip-¯op motion. To some extent the ¯uctuation of the torsional angle is related with the small ¯uctuation in the rotational motion of the chain around the chain axis (or the ¯uctuation of the setting angle). The thermal parameters or the meansquare displacements from the equilibrated position were calculated for carbon and hydrogen atoms. Fig. 15 gives the calculated results. As the temperature increased, the
mean-square displacements became larger, and increased dramatically above 350 K. The X-ray structure analysis reported in Ref. [22] showed that the thermal parameters increased with the de¯ection point around 300 K. Although the temperature is a little different, a good agreement can be seen between the observed and calculated results. From Figs. 13±15, the origin of the remarkable increment of the mean-square atomic displacements is considered to be the synergistic effect of the intramolecular torsional motion and the ¯ip-¯op-type rotational motion of the chains. According to the results of the MD calculation made by Noid et al. [39±41] for the aggregation models of the chains of ®nite length without the three-dimensional periodic boundary conditions, the generation of gauche bonds could be seen at high temperatures. In our present calculation such largely twisted segments could not be detected even at 400±500 K because of the imposition of the threedimensional periodic boundary condition. That is, the chain segments were in®nitely connected by the covalent bonds along the c-axis and so the chain motion is more or less constrained. Originally the chains in the crystal should be changed to the random coils at a temperature as high as 500 K, because the equilibrium melting temperature of the PE crystal is about 410 K. However, by connecting the chain segments included in the neighboring MD cells by strong covalent bonds, the chain orientations were
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
367
Fig. 11. Time dependence of the setting angles calculated at 400 K for the various chains surrounding chain A as shown in Fig. 9.
maintained as seen in Figs. 1, 7 and 8, and the planar-zigzag chains were conformationally disordered as a result of the torsional ¯uctuation of the successive monomeric sequences. This conformational disordering became more remarkable at 400±500 K. The situation shown in Figs. 7 and 8 is apparently similar to the generation of the hexagonal phase in the temperature region immediately below the melting point. As already clari®ed experimentally, the transition from orthorhombic to hexagonal phase can be observed when the ultra-drawn PE ®ber is tightly wound around some metal holder and heated up 420 K beyond the equilibrium melting point of the orthorhombic PE crystal [35,46]. The hexagonal phase consists of the conformationally disordered chains. The tensile stress (or the hydrostatic pressure) is required to generate the hexagonal phase. This is important, because we must keep the molecular chains straight along the chain axis to allow the chains to experience the conformational disordering without any drastic contraction to the random coils in the molten state. In the temperature region of 400±500 K, the MD calculations gave structures with highly disordered chain conformations with the orientation along the c-axis as the
three-dimensional periodic boundary conditions were imposed. This situation is quite similar to the experimental observation of the generation of the hexagonal phase caused under the tensile stress. It must be noticed, however, that the structural disordering observed at such high-temperature regions does not necessarily mean that the phase transition takes place from the orthorhombic to the hexagonal phase. Here we have pointed out only the possibility of such a structural change accompanying the conformational disordering phenomenon. 4.3. Temperature dependence of Young's modulus along the chain axis The elastic constants of the crystal lattice were calculated at various temperatures by the method of Parrinello and Rahman [9] and Ray [10]. The contribution of the potential energy term (the second derivatives with respect to the strain) is overwhelmingly large compared with the other terms such as the kinetic energy term and the thermal ¯uctuation term [12,15]. Besides, the calculation of the second derivatives of the potential energy has a rather
368
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Fig. 12. Time dependence of the relative height calculated at 400 K for the various chains surrounding chain A as shown in Fig. 9.
large error in the actual calculation, about ^10% of the averaged value. Therefore we calculated only the term of the second derivatives of the potential energy, in particular Young's modulus along the chain axis. Many snapshot structures were extracted from the trajectory data and Young's moduli were calculated. The calculated Young's modulus along the chain axis is plotted against temperature as shown in Fig. 16(a). Young's modulus at 0 K is 330 GPa. This value is in good agreement with the calculated values reported in Refs. [1,2,14,16,22]. Young's modulus at a low temperature was evaluated as 235±254 GPa by the X-ray diffraction method (120 K) [47]. The bulk Young's modulus of the ultra-drawn ®ber was reported to be 288^10 GPa at liquid nitrogen temperature [48]. Young's modulus of the crystalline region at a low temperature should be higher than this. As the temperature is increased, the modulus decreases gradually. At room temperature the calculated value is about 280 GPa. The observed values at room temperature are 280±305 GPa (from the Raman scattering measurement of the LAM of the n-alkane chain [17±20]), 329 GPa (from the coherent inelastic neutron scattering measurement of the dispersion curve of the acoustic phonon mode [49]), and 210±254 GPa (from the X-ray diffraction measurement under the tensile stress [5,6]). As already pointed out, the
X-ray diffraction method contains the problem on stress distribution within the bulk sample. The evaluated modulus is different depending on the samples and the researchers. Young's modulus estimated after the correction of the heterogeneous stress distribution is about 260 GPa [50]. The neutron scattering method contains a rather large experimental error because of the low signal-to-noise ratio of the scattering intensity. Compared with these two methods, the Raman scattering method seems more reliable. The LAM frequency of the n-alkane single crystal is determined by the coupling of intra and intermolecular interactions. The calculated value of 280 GPa is close to the value observed by this LAM method, 280±305 GPa, corrected for the intermolecular interactions. Palmo and Krimm reported the modulus 303 GPa calculated at room temperature [16], where they reduced the value at 0 K (330 GPa) by changing the effective cross-sectional area and length of the chain to the values observed experimentally at room temperature. Our value of 280 GPa is a little lower than theirs. One of the reasons might be the difference in the used potential function parameters. Another reason is an error in the evaluation of the modulus. As mentioned above, the calculation error is estimated to be ca. 10%. Therefore it is a little dif®cult to say de®nitely that our 280 GPa is really different from their value of 303 GPa beyond the estimation error.
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
369
Fig. 13. Distribution of the torsional angle of a skeletal chain calculated at the various temperatures.
Anyway, Young's modulus at room temperature decreases from the value at 0 K by ca. 40^10 GPa. The factors governing the temperature dependence of Young's modulus along the chain axis are now discussed. Fig. 16 shows the temperature dependencies of: (a) Young's modulus; (b) the averaged torsional angle about the C±C bond; (c) the ¯uctuation of the skeletal torsional angle from the averaged value; and (d) the percentage of the chain contraction from the repeating period of the planar-zigzag form (0 K). The averaged torsional angle of the skeletal chain is almost 1808. Therefore the molecular chain takes essentially the planar-zigzag conformation in average at any temperature. However, the atomic displacements from the average positions become larger at higher temperatures, as seen in the increasing torsional ¯uctuation. The deviation of the torsional angle from 1808 means the contraction of the chain axis. Therefore the ¯uctuation of the skeletal torsional angle and the chain contraction have a good correlation with each other.
In Fig. 16 we also notice the intimate relationship between Young's modulus and the chain contraction. Fig. 17 shows the plot of Young's modulus against the chain contraction calculated at various temperatures. It should be emphasized here that Young's modulus decreases dramatically by quite a small degree of chain contraction from the perfectly planar-zigzag form: only 0.3% contraction gives the decrement of the modulus by 35%! This situation was already reported by Tashiro and Tadokoro for the case of nylon 6 a form [51]. Nylon 6 chain takes the planar-zigzag conformation in average, but the repeating period becomes shorter at higher temperatures because of the torsional motion around the bond between the amide group and the alkyl chain segment. Young's modulus along the chain axis decreases by 45% from 312 to 172 GPa by the chain contraction of only 0.5%! This remarkable change of the modulus was experimentally con®rmed by Miyasaka et al. [52] and Nakamae et al. [53].
370
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Fig. 14. Comparison of the time dependence of the setting angle, torsional angle and relative height calculated for chain A in Fig. 9 at 350 and 400 K.
As reported in Ref. [29], we also made the MD calculation on the temperature dependence of the chain conformation and Young's modulus along the chain axis for trigonal POM crystal. The details were described in a previous paper. Fig. 18 shows the conformational disordering of the POM chain. At higher temperatures the helical chain conformation is more disordered and the chain is
Fig. 15. Temperature dependence of the mean square displacements calculated for C and H atoms of orthorhombic PE crystal.
contracted as seen in Fig. 19, where the plot of the calculated Young's modulus and the chain contraction is made against temperature. As the temperature increases Young's modulus decreases gradually, being in good agreement with the data observed using X-ray diffraction measurement [7,54]. The contraction of the chain length by 2% results in the decrease of Young's modulus by about 50%. We now discuss a possibility of the origin of the so-called a relaxation process observed for the temperature dependence of the mechanical property of the bulk PE sample. We know the remarkable dispersion of the bulk Young's modulus in the temperature region higher than room temperature [55]. Several relaxation mechanisms were proposed for this mechanical property change [36,37,56±63]: the chain slippage between the neighboring crystallites; the anharmonic chain motion in the crystallites; and so on. The present MD calculation is considered to give a concrete and quantitative image to the anharmonic chain motion as a possible origin of the a relaxation mechanism. A similar concept was reported by Mori et al. [64], although their model was not so large as ours, having a possibility to be affected by the periodic boundary conditions as discussed in the previous
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
371
Fig. 16. Temperature dependencies of: (a) Young's modulus along the chain axis; (b) averaged torsional angle of the skeletal chain; (c) the ¯uctuation of the torsional angle from the averaged value; and (d) the chain contraction calculated for orthorhombic PE crystal.
section, and the reasonableness of the used potential functions were not checked. Of course, we do not consider that the present calculation could give the correct and ®nal answer to the a relaxation process of PE. In order to clarify the real image of this process, we are now calculating the relaxation time for the thermal motions of ¯ip-¯op, translation, and skeletal twisting of the PE chains, which might be compared with the observed relaxation time. 5. Conclusions
Fig. 17. Correlation between Young's modulus and the contraction of the molecular chain calculated for orthorhombic PE crystal at the various temperatures (see Fig. 16).
In the present paper we have mainly described the results of the MD calculations made for the orthorhombic PE crystal model consisting of 5184 atoms under the threedimensional periodic boundary condition. The cell constants were found to correspond well with the actually observed data. In particular, the linear thermal expansion coef®cients
372
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Fig. 18. Molecular conformations in trigonal POM cystal calculated at the various temperatures.
were in good agreement with the experimental values. The molecular chains were found to show the ¯ip-¯op-type motion at high temperatures, and at the same time the torsional ¯uctuation of the skeletal chains became large. These thermal motions re¯ect on the large thermal parameters of the H and C atoms at high temperature. The thermal parameters show the de¯ection point in the vicinity of 350 K, corresponding relatively well to the X-ray analysis result. Young's modulus along the chain axis was evaluated by the second derivatives method of the potential function. The modulus at room temperature, ca. 280 GPa, is in good agreement with the experimental value 280±290 GPa obtained by the LAM measurement of the n-alkane crystals. Young's modulus was found to closely relate with the chain contraction, which was caused by the thermal ¯uctuation of the skeletal chain conformation. Young's modulus was found to decrease by 35% for quite a small contraction of the planar-zigzag chain by 0.3%. A similar situation was also found in the case of nylon 6 and trigonal POM. This way, Young's modulus along the chain axis is affected sensitively by the slight change in the chain conformation. Therefore, if we plan to produce a new polymer material with high thermal stability of the mechanical properties, we will still have to ®nd a way to ®x the chain length at higher temperatures. Acknowledgements The author wishes to thank Dr Tahir Cagin of Molecular Simulations Incorporation for his kind discussion about the numerical calculation of the elastic constants.
Fig. 19. Temperature dependence of Young's modulus along the chain axis and the c-axial length calculated for trigonal POM crystal. The open circles in the top ®gure are the observed data by the X-ray method. The bottom ®gure shows the correlation between Young's modulus and the chain contraction extracted from the above two ®gures.
functions and the re®ned parameters are listed. The total energy of the model is expressed as a summation of the following potential functions. EBond Stretch exp 2ab R 2 Rb 2 12 ; kb 2Db a2b EAngle Bend Ccos u 2 cos ua 2 =2; ku C sin2 ua ETorsion Vt 1 1 cos 3f=2
Appendix A. Potential functions used in the present MD calculation and the numerical values of the related parameters
EBond-Angle D cos u 2 cos ua R 2 Rb ; kru 2D sin ua
As stated in the text the potential functions proposed by Karasawa et al. [27] were re®ned so that the vibrational spectroscopic data could be reproduced as well as possible. In this appendix the concrete forms of these potential
EOne-Center Angle±Angle Cross-Terms
EBond±Bond krr R1 2 R1b R2 2 R2b
G cos uIJK 2 cos uaIJK cos uIJL 2 cos uaIJL
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
ETwo-Center Angle±Angle Cross-Terms F cos f cos uIJK 2 cos uaIJK cos uJKL 2 cos uaJKL Eelectrostatic Qi Qj = 10 1RIJ EvdW A exp 2BRij 2 C=R6ij DV 6 exp z 1 2 r 2 zr26 = z 2 6 van der Waals parameters for C and H atoms Karasawa (MCXX) Tashiro zC 13.0 13.0 3.8410 3.8410 RVC 0.07918 0.0792 DVC zH 11.2 11.2 3.1665 3.1665 RVH 0.0200 0.0190 DVH Valence force ®eld parameters for C and H atoms Bond stretch C±C
Rb kb Db Rb kb Db
1.4841 884.9940 85.80 1.0765 741.3720 95.10
1.5300 884.9940 85.80 1.0765 580.0000 90.00
ku ua ku ua ku ua
55.6076 119.3933 65.7301 117.7291 84.1810 121.2400
48.0000 119.3933 72.0000 117.7291 84.1810 121.2400
Vt Vt Vt
5.1686 6.1626 5.7070
5.1686 6.1626 5.7070
222.6583 3.1321 234.3195 225.9234 1.3684 254.0185 26.2187
222.6583 3.1321 234.3195 225.9234 1.3684 254.0185 26.2187
One-center angle±angle cross-terms GCC:HC GCC:HH GCH:CC GCH:CH
27.6083 25.3356 25.0824 25.3356
27.6083 25.3356 25.0824 25.3356
Two-center angle±angle cross-terms FH:C:C:H FC:C:C:H FC:C:C:C
217.7274 216.4004 221.5910
217.7274 216.4004 221.5910
C±H
Angle bend H±C±H C±C±H C±C±C Torsion H±C±C±H C±C±C±H C±C±C±C
Angle cross-terms H±C±H DH u kHH C±C±H DC u DH u kHC C±C±C DC u kCC
References [1] Tashiro K. Kobunshi Ronbunshu (Jpn J Polym Sci Technol) 1992;49:711. [2] Tashiro K. Prog Polym Sci 1993;18:377. [3] Tashiro K, Kobayashi M. Polymer 1995;37:1775. [4] Jungnitz S, Jakeways R, Ward IM. Polymer 1986;27:1651.
373
[5] Matsuo M, Sawatari C. Macromolecules 1988;21:1653. [6] Nishino T, Ohkubo H, Nakamae K. J Macromol Sci Phys B 1992;31:191. [7] Wu G, Tashiro K, Kobayashi M, Komatsu T, Nakagawa K. Macromolecules 1989;22:758. [8] Wu G, Tashiro K, Kobayashi M. Macromolecules 1989;22:188. [9] Parrinello M, Rahman A. J Chem Phys 1982;76:2662. [10] Ray JR. Comp Phys Rep 1988;8:109. [11] Tashiro K, Kobayashi M, Tadokoro H. Polym J 1992;24:899. [12] Cagin T, Kawasawa N, Dass Gupta S, Goddard III WA. In: Mark JE, editor. Computational methods in materials science, 1992. [13] Kobayashi M. J Chem Phys 1979;70:509. [14] Lacks DJ, Rutledge GC. J Phys Chem 1994;98:1222. [15] Cagin T, Karasawa N, Dasgupta S, Goddard III WA. Mater Res Soc Symp Proc 1992;278:61. [16] Palmo K, Krimm S. J Polym Sci: Polym Phys Ed 1996;34:37. [17] Strobl GR, Eckel R. J Polym Sci: Polym Phys Ed 1976;14:913. [18] Kobayashi M, Sakagami K, Tadokoro H. J Chem Phys 1983;78:6391. [19] Chang C, Krimm S. J Polym Sci: Polym Phys Ed 1979;17:2163. [20] Snyder RG, Strauss HL, Alamo R, Mandelkern L. J Chem Phys 1994;100:5422. [21] Hsu SL, Ford GW, Krimm S. J Polym Sci: Polym Phys Ed 1977;15:1769. [22] Tashiro K, Ishino K, Ohta T. Polymer 1999;40:3469. [23] Kitagawa T, Miyazawa T. Adv Polym Sci 1972;9:335. [24] Weiner JH. Statistical mechanics of elasticity. New York: Wiley, 1983. [25] Kobayashi M, Tadokoro H. J Chem Phys 1977;66:1258. [26] Kobayashi M. J Chem Phys 1979;70:4797. [27] Karasawa N, Dasgupta S, Goddard III WA. J Phys Chem 1991;95:2260. [28] Palmo K, Mirkin NG, Pietila LO, Krimm S. Macromolecules 1993;26:6831. [29] Tashiro K, Kobayashi M. Rep Prog Polym Phys Jpn 1996;39:349. [30] Sun H. J Phys Chem B 1998;102:7338. [31] Tashiro K, Sato M, Cho Y, Kobayashi M. Polym Prepr Jpn 1999;48:863. [32] Tashiro K, Sato M, Cho Y, Kobayashi M. Polym Prepr Jpn 1999;48:3574. [33] Tashiro K, Satoh M, Ike M. Polym Prepr Jpn 2000;49:2541. [34] Gasteiger J, Marsili M. Tetrahedron 1980;36:3219. [35] Pennings AJ, Zwijnenburg A. J Polym Sci: Polym Phys Ed 1979;17:1011. [36] Syi JL, Mans®eld ML. Polymer 1988;29:987. [37] Mans®eld ML, Boyd RH. J Polym Sci: Polym Phys Ed 1978;16:1227. [38] Tadokoro H. Fiber diffraction methods. In: French AD, Gardner KCH, editors. ACS Symposium Series 141Washington, DC: American Chemical Society, 1980. [39] Sumpetr BG, Noid DW, Wunderlich B. J Chem Phys 1990;93:6875. [40] Noid DW, Sumpter BG, Wunderlich B. Macromolecules 1990;23:664. [41] Liang GL, Noid DW, Sumpter BG, Wunderlich B. J Polym Sci: Polym Phys Ed 1993;31:1909. [42] Liang GL, Noid DW, Sumpter BG, Wunderlich B. Acta Polym 1993;44:219. [43] Liang GL, Noid DW, Sumpter BG, Wunderlich B. J Phys Chem 1994;98:11 739. [44] Tuzun RE, Noid DW, Sumpter BG. Makromol Chem: Theory Simul 1998;7:203. [45] Kobayashi M, Tadokoro H. Macromolecules 1975;8:897. [46] Tashiro K, Sasaki S, Kobayashi M. Macromolecules 1996;29:7460. [47] Nakamae K, Nishino T, Miyazaki H. Polym Prepr Jpn 1999;48:914. [48] Barham PJ, Keller A. J Polym Sci: Polym Lett Ed 1979;17:591. [49] Holliday L, White JW. Pure Appl Chem 1971;26:545. [50] Tashiro K, Wu G, Kobayashi M. Polymer 1988;29:1768. [51] Tashiro K, Tadokoro H. Macromolecules 1981;14:781. [52] Miyasaka K, Isomoto T, Koganeya K, Ishkawa K. J Polym Sci: Polym Phys Ed 1980;18:1047.
374
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
[53] Nakamae K, Nishino T, Hata K, Matsumoto T. Koubunshi Ronbunshu (Jpn J Polym Sci Technol) 1987;44:421. [54] Nakamae K, Nishino T, Shimizu Y, Hatano K. Polymer 1990;31:1909. [55] Ohta Y, Yasuda H. J Polym Sci B: Polym Phys 1994;32:2241. [56] Kajiyama T, Okada T, Sakoda A, Takayanagi M. J Macromol Sci Phys B 1973;7:583. [57] Kyu T, Yasuda N, Suehiro S. Polym J 1976;8:565. [58] Suehiro S, Yamada T, Inagaki H, Kyu T, Nomura S, Kawai H. J Polym Sci: Polym Phys Ed 1979;17:763.
[59] Suehiro S, Yamada T, Kyu T, Fujita K, Hashimoto T, Kawai H. Polym Engng Sci J 1979;19:929. [60] Suehiro S, Kyu T, Fujita K, Kawai H. Polym J 1979;11:331. [61] Ogita T, Yamomoto R, Suzuki N, Ozaki F, Matsuo M. Polymer 1991;32:822. [62] Usami T, Takayama S. Macromolecules 1984;17:1756. [63] Matsuo M, Sawatari C, Ohhata T. Macromolecules 1988;21:1317. [64] Mori S, Ohta Y, Aikawa Y. Polym Prepr Jpn 1994;43:3419. [65] Avitabile G, Napolitano R, Pirozzi B, Rouse KD, Thomas MW, Willis BTM. J Polym Sci B 1975;13:351.
www.elsevier.com/locate/ctps
Molecular dynamics calculation to clarify the relationship between structure and mechanical properties of polymer crystals: the case of orthorhombic polyethylene K. Tashiro* Department of Macromolecular Science, Graduate School of Science, Osaka University, 1-1 Machikaneyama-cho, Toyonaka, Osaka 560-0043, Japan Received 21 September 2000; revised 20 December 2000; accepted 26 December 2000
Abstract The molecular dynamics (MD) technique was used to calculate the temperature dependence of the structure, molecular motion, and mechanical property of the orthorhombic polyethylene (PE) crystal. The potential functional parameters reported by Karasawa et al. (J Phys Chem, 95 (1991) 2260) were re®ned further so that the vibrational frequencies of infrared and Raman bands, measured by us at ultra-low temperatures for the normal and fully deuterated PE, could be reproduced well. The ¯ip-¯op motion around the chain axis and the torsional motion of the skeletal chains were found to start above ca. 350 K and increase the amplitude of these motions progressively. Coupling these two types of chain motion resulted in a steep increase of the thermal vibration parameters or the mean-square-displacements of carbon and hydrogen atoms, corresponding well with the X-ray data. The lattice constants and the related linear thermal expansion coef®cients were also found to be in good agreement with the observed data. The calculated Young's modulus along the chain axis decreased gradually with the increasing temperature: 330 GPa at 0 K to 280 GPa at room temperature. The latter was in good agreement with the value of 280±305 GPa evaluated from the Raman measurement of the longitudinal acoustic mode. Young's modulus was found to relate intimately with the chain contraction caused by the skeletal torsional motion. Only 0.3% contraction of the chain resulted in the reduction of the modulus by ca. 35%. A similar behavior was also seen in the trigonal polyoxymethylene and nylon 6 a forms. q 2001 Elsevier Science Ltd. All rights reserved. Keywords: Molecular dynamics; Polyethylene; Crystal structure; Molecular motion; Young's modulus; Potential functions
1. Introduction In a series of papers [1±3] we have discussed the intimate relationship between the structure and mechanical properties of polymer crystals from the molecular level by way of theoretical calculations and related experiments. At the present stage we may state that, as long as the structure of a polymer substance is known reasonably well, Young's moduli in various directions of the polymer crystal can be predicted with fairly high probability. In these evaluations, however, we have not explicitly taken into consideration the effects of the thermal motion of the chains on the mechanical properties of the polymer crystals. The mechanical properties of polymer crystals may be assumed as the limiting values of the bulk samples. The mechanical properties of the bulk samples are a function of temperature. Therefore it is absolutely and essentially important to evaluate the effect of the thermal motion in the theoretical prediction of the mechanical properties of polymer crystals. Such trials were * Tel.: 181-6-6850-5455; fax: 181-6-6850-5455. E-mail address: [email protected] (K. Tashiro).
made by several researchers from both the experimental [4±8] and theoretical [9±16] points of view. The temperature dependence of the elastic constants of polyethylene (PE) crystal was calculated on the basis of the free energy treatment under the quasi-harmonic approximation [13,14], molecular dynamics (MD) method [12,15], and lattice dynamics [16]. On the experimental side, the temperature dependence of Young's modulus along the chain axis (or the crystallite modulus) was measured for PE by the X-ray diffraction method [4±6]. In this case, however, the situation became rather complicated because of the problem of the so-called stress distribution in the bulk sample, giving some uncertainty to the experimental values obtained by the X-ray diffraction method. The temperature dependence of the vibrational frequency of the longitudinal acoustic mode (LAM) may be another useful method, but requires correction to study the effect of the intermolecular interactions [17±20], the mass effect at the end groups of the chains [21], and so on. In the evaluation of the mechanical properties of polymer crystals as a function of temperature, we also need to consider the temperature dependence of the crystal structure.
1089-3156/01/$ - see front matter q 2001 Elsevier Science Ltd. All rights reserved. PII: S 1089-315 6(01)00005-8
358
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Another motivation of this paper is based on our recent study on the temperature dependence of the crystal structure of orthorhombic PE. In a previous paper [22], we carried out the X-ray diffraction measurements of highly oriented PE samples at various temperatures from 2100 to 1908C and could collect numerous X-ray re¯ections at each temperature by using the imaging plate system. By analyzing the obtained sets of re¯ection data, the temperature dependence of the crystal structure was revealed with high reliability. The lattice constants, particularly the a-axial length showed the de¯ection point around 1108C. The setting angle of the planar-zigzag chain measured from the b-axis and the thermal parameters of the constituent atoms C and H also showed similar de¯ection points in the same temperature region. The atomic thermal parameters or the mean-square amplitudes of the atomic vibrations were calculated by Kitagawa and Miyazawa [23] under the quasi-harmonic approximation, and were found to be in good agreement with the above-mentioned X-ray values in the low-temperature region. In the temperature region higher than 108C, however, a large discrepancy was detected between the observed and the calculated thermal parameters. This discrepancy led us to speculate that the large-amplitude thermal motion of the chains become signi®cant in the temperature region above 108C. However, we could not explain the type of thermal motion that occurs in the hightemperature region. This needs to be investigated from a theoretical point of view. In the theoretical calculation of the temperature dependence of the structure, thermal motion, elastic constants and so on, several methods have been presented as mentioned above. One is to treat the free energy of the crystal as a function of temperature [13,14], where the elastic constants are obtained as second derivatives of the free energy with respect to the stress component; and the other method is to introduce the anharmonic terms into the force constants: f f0 1 f 0 DX where f, f0, and f 0 are the force constant and the corresponding harmonic and anharmonic terms, respectively, and DX is the atomic displacement caused by the external perturbation. Phenomenological treatment of the anharmonic effect is to use the mode-GruÈneisen constants, which are de®ned as g 2 dn=d10 =n where n is the vibrational frequency of the mode and 1 is the strain [11,24]. The MD method was also used for the evaluation of the temperature dependence of the elastic constants [9,10]. The MD calculation on the temperature dependence of the lattice and elastic constants of orthorhombic PE crystal was made by Cagin et al. [12,15]. One serious problem in these trials is the reasonableness of the utilized potential function parameters, which could not be said to be so reasonable as to give a satisfactorily good reproduction of the observed vibrational frequencies of the infrared and Raman spectra. Kobayashi et al. re®ned the force constants of the valence force ®eld type so that the calculated vibrational frequencies could be reproduced as reasonably as possible [25,26]. Karasawa et al. presented a set of potential function
parameters (MCXX) [27] and calculated the vibrational frequencies and elastic constants of orthorhombic PE crystal. Palmo et al. re®ned the molecular mechanics potential function for n-alkane chain base on scaled ab initio force ®elds [28]. In spite of these efforts, however, the discrepancy between the calculated and observed vibrational frequencies has not yet been removed satisfactorily. As reported in the present paper, we re®ned the potential function parameters proposed by Karasawa et al. and applied them to the MD calculation of the orthorhombic PE crystal to evaluate quantitatively the thermal behaviors of lattice and elastic constants, particularly Young's modulus along the chain axis. As mentioned above, the MD calculation has already been made by Cagin et al. [12,15], but the total number of atoms included in the crystal lattice was limited to 360±504. These atom numbers may not be large enough to simulate the thermal behavior of the actual crystal, although they check the effect of the model size on the calculated results. In the present study we built the PE lattice model consisting of ca. 5200 atoms in total and calculated the trajectory (the coordinates and velocities of the constituent atoms and the lattice constants) over the long time span of 100 ps under the three-dimensional periodic boundary conditions. On the basis of the trajectory data obtained at the various temperatures, Young's modulus along the chain axis was calculated and compared with the observed data. It took about 50 days to obtain the trajectory ®le at one temperature point because of the limited power of our computer system. It took us as long as one and a half years, inclusive of the time taken for the data analysis, to obtain the results presented in this paper. A similar calculation was also made for trigonal polyoxymethylene (POM). As the details have already been already reported in Ref. [29], only the results will be compared with those of orthorhombic PE; and some common features on the relationship between the structure and mechanical property of polymer crystals will be extracted. 2. Re®nement of potential function parameters On the basis of the vibrational spectroscopic data, the force constants or the potential function parameters concerning the intra and intermolecular interactions were proposed for PE crystals [25±27]. As already mentioned above, Karasawa et al. proposed a set of potential functions (MCXX), which gave a relatively good agreement between the observed and calculated vibrational frequencies for orthorhombic PE crystal. Recently Sun reported the potential function parameters, COMPASS, which gave a good reproduction of the equations of state for many organic substances [30]. By using their proposed parameters we calculated the vibrational frequencies for orthorhombic PE crystals of both the hydrogenous and the deuterated species and compared the results with the observed values, though the comparison was not satisfactory. We then tried to
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
359
Table 1 Comparison between the calculated and observed values for the various parameters Unit cell parameters
Karasawa et al. [27]
Ê) a(A Ê) b (A Ê) c (A Vibrational spectra (cm 21) CH2
7.202 4.795 2.546 Lattice mode Lattice mode Lattice mode Rocking Bending
CD2
Lattice mode Rocking Bending
Young's modulus (GPa)
a
Karasawa et al. 84 115 132 784 798 1484 1485 ± ± ± ± ±
a-axis b-axis c-axis
This study 6.995 4.733 2.531
Observed (4 K) [65] 7.121 4.851 2.548
This study 97.1 108.2 123.3 718.7 735.2 1494.2 1496.7 89.4 513.9 525.5 1089.0 1093.0
Observed (7 K) 79.1 108.5 137.0 721.8 734.7 1460.8 1474.8 74.5 519.1 529.1 1083.8 1094.5
This study 10.9 9.9 329.5
Observed 10.2 (110) a ± See text
Ref. [47] (Young's modulus along 110 direction).
modify the potential function parameters proposed by Karasawa et al. further so that the vibrational spectroscopic data observed for various kinds of samples such as the hydrogenous PE, the deuterated PE, the random copolymers of H and D ethylene monomeric units, and the PE blends of the pure H and D species could be reproduced consistently and satisfactorily [31,32]. The re®nement of the potential function parameters was made by a trial-anderror method. The calculation of the vibrational frequencies was carried out for the energetically minimized structure models by using the commercial software Cerius 2 (Version 3.8, Molecular Simulations Incorporation) on the basis of the Hessian-matrix method. The set of potential function parameters obtained is listed in Appendix A, where the de®nitions of the potential functions and the corresponding parameters are given concretely. The comparison between the observed and calculated values is shown in Table 1, where the lattice constants, the vibrational frequencies, and the crystallite moduli are listed. We have recently collected the observed vibrational data at 7 K [33]. As seen in Table 1, the values calculated by the potential function parameters listed in Appendix A are in good agreement with the observed data. The detailed discussion on the vibrational data will be reported in a separate paper [33]. The discussion about the calculated Young's modulus will be made in a later section of the present paper.
(Version 3.8, Molecular Simulations Incorporation). In an MD cell all the 72 chains were packed, each of which was constructed by 24 CH2 monomeric units. Therefore the basic MD cell contained 5184 atoms in total. The threedimensional periodic boundary conditions were applied: the 24 CH2 chain segments were connected in®nitely along the c-axis by covalent bonding. Fig. 1(a) shows the thus-constructed starting model. By applying the abovementioned potential function parameters, the atomic charges of 20.0533 and 10.0267 electron units were assumed for all the carbon and hydrogen atoms, respectively, which were calculated by the Gasteiger method [34]. The MD calculation was made for the energetically minimized model under the constant NPT condition by the Hoover±Nose method, where the temperature (T ) was ®xed to the predetermined value (0, 100, 200, 300, 350, 400, and 500 K), and the external pressure was assumed to be 0 atm. The time interval of each step was 1 fs and the total time taken for calculation was 100 ps at each temperature. The trajectory data created in the MD calculations were saved at every 0.1 ps and analyzed by using our home-made programs to obtain the various physical quantities described in the following sections.
4. Results and discussion 4.1. Temperature dependence of lattice constants
3. Molecular dynamics calculations The MD calculation was performed by using the Cerius 2
The lattice constants a, b, and c were evaluated by averaging the trajectory data in the time region reaching the constant energy. The averaged angles a , b , and g were
360
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Fig. 3. Comparison of the temperature dependence of the b-axial length between the observed and calculated values. The horizontal bars in this ®gure indicate the temperature regions where the linear thermal expansion coef®cients were calculated.
Fig. 1. Snapshots of PE crystal model viewed along the chain axis at: (a) 0, and (b) 300 K.
Fig. 2. Comparison of the temperature dependence of the a-axial length between the observed and calculated values. The horizontal bars in this ®gure indicate the temperature regions where the linear thermal expansion coef®cients were calculated.
essentially 908. Figs. 2±4 show the temperature dependencies of the calculated a-, b-, and c-axial lengths. The absolute values of the lattice constants a, b, and c were not suf®cient to reproduce the experimental values accurately. This disagreement could probably be because of the unsatisfactory re®nement of the used potential function parameters, which is still being made to obtain a better agreement between the observed and calculated values. At the same time we cannot ignore the conditions of the experiment. It is a well-known fact that the lattice constants of polymers are dependent sensitively on the sample preparation conditions. Therefore it may be safer to compare the calculated and observed cell constants by taking into account some amount of ¯exibility. The thermal expansion itself is found to show almost the same behavior irrespective of the sample except
Fig. 4. Comparison of the temperature dependence of the c-axial length between the observed and calculated values. The horizontal bars in this ®gure indicate the temperature regions where the linear thermal expansion coef®cients were calculated.
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
361
Fig. 5. Time dependence of the setting angle of an arbitrarily chosen chain calculated at the various temperatures.
for the extremely different sample preparation conditions. Therefore, here we will focus our attention mainly on the thermal expansivity. The a-axis shows the most remarkable temperature dependence, which is not linear over the whole temperature range. The linear thermal expansion coef®cient a a was calculated from the slope of the curve in each temperature region and compared with the observed values as shown in Fig. 2. In Ref. [22] the temperature dependence of a- and b-axes was measured, but the accuracy was not very high. The measurement was made again by increasing the sample-to-camera distance for the same ultra-drawn sample reported in the previous paper. Although the aaxial length itself is a little different, the calculated a a is in good agreement with the observed value in a wide
temperature region. The temperature dependence of the baxis is one digit smaller than that of the a-axis. The linear thermal expansion coef®cient a b is compared with the observed value in Fig. 3. The b-axis is contracted in the high-temperature region, and therefore the a b is negative, which is also reproduced relatively well by the calculation. The c-axis along the chain direction contracts by increasing the temperature as shown in Fig. 4, where the linear thermal expansion coef®cient a c is negative and two digits smaller than a a. The agreement between the observed and calculated values is good, ca. 25 £ 10 26 K 21 at 250 K. When the temperature dependencies of the a-, b- and caxes are investigated carefully, the de¯ection points are found around 350 K although the degree of the de¯ection
362
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Fig. 6. Distribution of the setting angle of the chain extracted from Fig. 5.
is not very high. This might correspond to the actually observed phenomenon: the a-axial length, e.g. increased linearly with temperature but changed its slope gently above the room temperature [22]. The degree of this de¯ection is dependent on the sample. The uniaxially oriented sample with draw ratio of ca. 15 times the original length shows the clearly detectable change around 108C, while this change is small for the ultra-drawn sample with the draw ratio of ca. 150 times. As discussed in Ref. [22], one possible origin of this different behavior may be the difference in the crystallite size between these two different sample types; the ultra-drawn sample should show quite a long and wide crystallite size [5,35]. We may say that the in®nitely large crystal of PE shows a small but signi®cant change in the curve of the a-axial length vs temperature around 108C, and that the difference in the degree of the curve de¯ection might come from the effect of the crystallite surface. For the smaller crystallite, the surface of the lamella gives stress on the structure and molecular motion of the inner part of
the crystallite [36,37]. In fact, several researchers including us clari®ed the change in the unit cell parameters and atomic thermal parameters depending on the samples with a different history [22,38]. By taking it into consideration that the three-dimensional periodic boundary conditions are assumed in the present MD calculations, the small but detectable change in the curves shown in Figs. 2±4 may correspond to the behavior of a very large crystal of orthorhombic PE. However, we need to clarify the concrete origin of these de¯ections in the cell parameters, which will be discussed in the following section. 4.2. Temperature dependence of molecular motion The main purpose of the present paper is to clarify the relationship between the temperature-dependent mechanical property and the conformational disordering of chains caused by thermal motion. As for the molecular motion of PE chains, a detailed discussion was made by Noid et al. in a
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
363
Fig. 7. A series of snapshots of PE crystal taken at 400 K in the time region of 9±12 ps. The molecular chains of red and yellow are named as chains A and B in the text, respectively.
364
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
series of papers on MD calculations [39±44]. Therefore, in the present paper, the description of the thermal motion of PE chains will be made brie¯y so that the information on the essential features of this thermal motion can be helpful for the discussion on the temperature dependence of mechanical property. Fig. 1(b) shows the snapshot taken at an arbitrary time of the trajectory data obtained at 300 K. The orientation of the molecular chains viewed along the chain axis was disordered and at the same time the relative height of the neighboring chains also changed. The setting angle, which is de®ned as an angle between the zigzag plane of the chain and the b-axis, was calculated for the arbitrarily chosen chain in the cell and was plotted against time as shown in Fig. 5. The distribution of the setting angles of the chain is shown in Fig. 6. In the temperature region up to 300 K, the setting angle is almost constant, about 438 in average, but the ¯uctuation of the setting angle from the averaged value increases with the increasing temperature. Around 350 K, as seen in Fig. 5, this ¯uctuation becomes much larger and some chains were found to change their setting angle from 143 to 243 and then to 1438, i.e. the ¯ip-¯op motion was
being observed. This motion becomes more remarkable at 400 K. At 500 K, the setting angles of the chains change almost randomly without any correlation among the neighboring chains. The chain motion at 400 K is now seen in more detail. Figs. 7 and 8 illustrate a series of snapshots picked up at various times. Around 10 and 60 ps, the cooperative rotations of the neighboring chains A (red) and B (yellow) can be observed; the naming of the chains is shown in Fig. 9. For example, around 9.5 ps as shown in Fig. 7, the yellow chain B gradually changed its setting angle from 143 to 2438 (or 11378) by rotating around the chain axis. With some time delay the red chain A started to change the setting angle around 10±10.5 ps. At 11±12 ps these two chains became parallel to each other and at the same time they were parallel to the other neighboring chains. That is, the local packing structure similar to the monoclinic form was generated in the orthorhombic crystal lattice. In Fig. 8 an inverse phenomenon was observed. Around 58 ps the chain B (yellow) changed its setting angle and the red chain A changed the orientation after a delay of about 5 ps. Then the local structure returned to the original orthorhombic type.
Fig. 8. A series of snapshots of PE crystal taken at 400 K in the time region of 50±60 ps. The molecular chains of red and yellow are named as chains A and B in the text, respectively. The small ®gure in the lower right side shows the clusters of monoclinic-type structure.
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
365
Fig. 9. Naming of the PE chains discussed in the text.
How about the correlation between the rotational motion and the translational motion of the chain? Fig. 10 shows the changes in the setting angle and the relative height plotted against time for the chains A and B. For example, at 5 ps, the setting angle of chain B changed steeply and the relative height (the z coordinate) changed correspondingly. Around 55 ps, the relative height changed again after the setting angle returned to its original value. That is, the molecular motion of chain B may be described as a helical motion. The period of this ¯ip-¯op motion of the PE chain is estimated to be about 100 ps at 400 K. However, such a motion was not observed clearly in chain A; around 60 ps the relative height was kept almost unchanged even when the setting angle changed drastically. Therefore, we may say that the rotational motion around the chain axis and the translational motion along the chain axis do not have an intimate correlation with each other. At the same time we should notice that, even when the two neighboring chains A and B changed their setting angles cooperatively, the other chains adjacent to these two chains did not show such ¯ip-¯op motion over the time region of at least 100 ps, as seen in Fig. 11. Besides,
the change in the relative height or the translational motion along the chain axis also did not show any correlation as shown in Fig. 12. This way the ¯ip-¯op motion and the translational motion of the zigzag chain seem to occur with very local and small correlations among the neighboring chains. Local correlation or formation of monoclinic-type clusters can be observed at the various parts of the crystal lattice (see the snapshot of the whole lattice taken at 60 ps in Fig. 9). It should be pointed out here that this monoclinic-type or parallel packing structure could be seen even at a relatively low temperature such as 300 K when the MD calculation was made for much smaller models of ca. 10±20 chains. This is considered to come from overly strong correlation between the chains or from the large effect of the lattice boundary. Besides, the chain packing structure of the monoclinic type is energetically more stable, though only slightly, and competes with the orthorhombic type of packing structure [45]. Therefore once the chain packing changes from the orthorhombic to monoclinic structure, the latter is stabilized and is kept unchanged over a long time. Such a situation was
366
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Fig. 10. Time dependence of the setting angle and relative height of the chains A and B at 400 K.
not seen in the present large model of ca. 5200 atoms. The ¯ip-¯op-type rotation of chains occurred not only in the boundary of the repeating unit cell but also at the central parts of the cell or at an arbitrary position of the cell. At 500 K the ¯uctuation of the setting angles of the chains became drastically large and an almost free rotation was attained at various places. Another important aspect is the relation of this rotational motion with the change in the chain conformation. We calculated the time dependence of the C±C torsional angles of the skeletal chains at various temperatures. The torsional angles ¯uctuated around the center of 1808. As shown in Fig. 13, the amplitude of ¯uctuation became larger as the temperature increased, but the averaged torsional angle was almost 1808 even at 500 K. In Fig. 14 the time dependence of the setting angle, torsional angle, and the relative height of the arbitrarily chosen chain is shown at 350 and 400 K. Any time correlation could not be found for the torsional motion of the chains with the large-amplitude ¯ip-¯op motion. To some extent the ¯uctuation of the torsional angle is related with the small ¯uctuation in the rotational motion of the chain around the chain axis (or the ¯uctuation of the setting angle). The thermal parameters or the meansquare displacements from the equilibrated position were calculated for carbon and hydrogen atoms. Fig. 15 gives the calculated results. As the temperature increased, the
mean-square displacements became larger, and increased dramatically above 350 K. The X-ray structure analysis reported in Ref. [22] showed that the thermal parameters increased with the de¯ection point around 300 K. Although the temperature is a little different, a good agreement can be seen between the observed and calculated results. From Figs. 13±15, the origin of the remarkable increment of the mean-square atomic displacements is considered to be the synergistic effect of the intramolecular torsional motion and the ¯ip-¯op-type rotational motion of the chains. According to the results of the MD calculation made by Noid et al. [39±41] for the aggregation models of the chains of ®nite length without the three-dimensional periodic boundary conditions, the generation of gauche bonds could be seen at high temperatures. In our present calculation such largely twisted segments could not be detected even at 400±500 K because of the imposition of the threedimensional periodic boundary condition. That is, the chain segments were in®nitely connected by the covalent bonds along the c-axis and so the chain motion is more or less constrained. Originally the chains in the crystal should be changed to the random coils at a temperature as high as 500 K, because the equilibrium melting temperature of the PE crystal is about 410 K. However, by connecting the chain segments included in the neighboring MD cells by strong covalent bonds, the chain orientations were
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
367
Fig. 11. Time dependence of the setting angles calculated at 400 K for the various chains surrounding chain A as shown in Fig. 9.
maintained as seen in Figs. 1, 7 and 8, and the planar-zigzag chains were conformationally disordered as a result of the torsional ¯uctuation of the successive monomeric sequences. This conformational disordering became more remarkable at 400±500 K. The situation shown in Figs. 7 and 8 is apparently similar to the generation of the hexagonal phase in the temperature region immediately below the melting point. As already clari®ed experimentally, the transition from orthorhombic to hexagonal phase can be observed when the ultra-drawn PE ®ber is tightly wound around some metal holder and heated up 420 K beyond the equilibrium melting point of the orthorhombic PE crystal [35,46]. The hexagonal phase consists of the conformationally disordered chains. The tensile stress (or the hydrostatic pressure) is required to generate the hexagonal phase. This is important, because we must keep the molecular chains straight along the chain axis to allow the chains to experience the conformational disordering without any drastic contraction to the random coils in the molten state. In the temperature region of 400±500 K, the MD calculations gave structures with highly disordered chain conformations with the orientation along the c-axis as the
three-dimensional periodic boundary conditions were imposed. This situation is quite similar to the experimental observation of the generation of the hexagonal phase caused under the tensile stress. It must be noticed, however, that the structural disordering observed at such high-temperature regions does not necessarily mean that the phase transition takes place from the orthorhombic to the hexagonal phase. Here we have pointed out only the possibility of such a structural change accompanying the conformational disordering phenomenon. 4.3. Temperature dependence of Young's modulus along the chain axis The elastic constants of the crystal lattice were calculated at various temperatures by the method of Parrinello and Rahman [9] and Ray [10]. The contribution of the potential energy term (the second derivatives with respect to the strain) is overwhelmingly large compared with the other terms such as the kinetic energy term and the thermal ¯uctuation term [12,15]. Besides, the calculation of the second derivatives of the potential energy has a rather
368
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Fig. 12. Time dependence of the relative height calculated at 400 K for the various chains surrounding chain A as shown in Fig. 9.
large error in the actual calculation, about ^10% of the averaged value. Therefore we calculated only the term of the second derivatives of the potential energy, in particular Young's modulus along the chain axis. Many snapshot structures were extracted from the trajectory data and Young's moduli were calculated. The calculated Young's modulus along the chain axis is plotted against temperature as shown in Fig. 16(a). Young's modulus at 0 K is 330 GPa. This value is in good agreement with the calculated values reported in Refs. [1,2,14,16,22]. Young's modulus at a low temperature was evaluated as 235±254 GPa by the X-ray diffraction method (120 K) [47]. The bulk Young's modulus of the ultra-drawn ®ber was reported to be 288^10 GPa at liquid nitrogen temperature [48]. Young's modulus of the crystalline region at a low temperature should be higher than this. As the temperature is increased, the modulus decreases gradually. At room temperature the calculated value is about 280 GPa. The observed values at room temperature are 280±305 GPa (from the Raman scattering measurement of the LAM of the n-alkane chain [17±20]), 329 GPa (from the coherent inelastic neutron scattering measurement of the dispersion curve of the acoustic phonon mode [49]), and 210±254 GPa (from the X-ray diffraction measurement under the tensile stress [5,6]). As already pointed out, the
X-ray diffraction method contains the problem on stress distribution within the bulk sample. The evaluated modulus is different depending on the samples and the researchers. Young's modulus estimated after the correction of the heterogeneous stress distribution is about 260 GPa [50]. The neutron scattering method contains a rather large experimental error because of the low signal-to-noise ratio of the scattering intensity. Compared with these two methods, the Raman scattering method seems more reliable. The LAM frequency of the n-alkane single crystal is determined by the coupling of intra and intermolecular interactions. The calculated value of 280 GPa is close to the value observed by this LAM method, 280±305 GPa, corrected for the intermolecular interactions. Palmo and Krimm reported the modulus 303 GPa calculated at room temperature [16], where they reduced the value at 0 K (330 GPa) by changing the effective cross-sectional area and length of the chain to the values observed experimentally at room temperature. Our value of 280 GPa is a little lower than theirs. One of the reasons might be the difference in the used potential function parameters. Another reason is an error in the evaluation of the modulus. As mentioned above, the calculation error is estimated to be ca. 10%. Therefore it is a little dif®cult to say de®nitely that our 280 GPa is really different from their value of 303 GPa beyond the estimation error.
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
369
Fig. 13. Distribution of the torsional angle of a skeletal chain calculated at the various temperatures.
Anyway, Young's modulus at room temperature decreases from the value at 0 K by ca. 40^10 GPa. The factors governing the temperature dependence of Young's modulus along the chain axis are now discussed. Fig. 16 shows the temperature dependencies of: (a) Young's modulus; (b) the averaged torsional angle about the C±C bond; (c) the ¯uctuation of the skeletal torsional angle from the averaged value; and (d) the percentage of the chain contraction from the repeating period of the planar-zigzag form (0 K). The averaged torsional angle of the skeletal chain is almost 1808. Therefore the molecular chain takes essentially the planar-zigzag conformation in average at any temperature. However, the atomic displacements from the average positions become larger at higher temperatures, as seen in the increasing torsional ¯uctuation. The deviation of the torsional angle from 1808 means the contraction of the chain axis. Therefore the ¯uctuation of the skeletal torsional angle and the chain contraction have a good correlation with each other.
In Fig. 16 we also notice the intimate relationship between Young's modulus and the chain contraction. Fig. 17 shows the plot of Young's modulus against the chain contraction calculated at various temperatures. It should be emphasized here that Young's modulus decreases dramatically by quite a small degree of chain contraction from the perfectly planar-zigzag form: only 0.3% contraction gives the decrement of the modulus by 35%! This situation was already reported by Tashiro and Tadokoro for the case of nylon 6 a form [51]. Nylon 6 chain takes the planar-zigzag conformation in average, but the repeating period becomes shorter at higher temperatures because of the torsional motion around the bond between the amide group and the alkyl chain segment. Young's modulus along the chain axis decreases by 45% from 312 to 172 GPa by the chain contraction of only 0.5%! This remarkable change of the modulus was experimentally con®rmed by Miyasaka et al. [52] and Nakamae et al. [53].
370
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Fig. 14. Comparison of the time dependence of the setting angle, torsional angle and relative height calculated for chain A in Fig. 9 at 350 and 400 K.
As reported in Ref. [29], we also made the MD calculation on the temperature dependence of the chain conformation and Young's modulus along the chain axis for trigonal POM crystal. The details were described in a previous paper. Fig. 18 shows the conformational disordering of the POM chain. At higher temperatures the helical chain conformation is more disordered and the chain is
Fig. 15. Temperature dependence of the mean square displacements calculated for C and H atoms of orthorhombic PE crystal.
contracted as seen in Fig. 19, where the plot of the calculated Young's modulus and the chain contraction is made against temperature. As the temperature increases Young's modulus decreases gradually, being in good agreement with the data observed using X-ray diffraction measurement [7,54]. The contraction of the chain length by 2% results in the decrease of Young's modulus by about 50%. We now discuss a possibility of the origin of the so-called a relaxation process observed for the temperature dependence of the mechanical property of the bulk PE sample. We know the remarkable dispersion of the bulk Young's modulus in the temperature region higher than room temperature [55]. Several relaxation mechanisms were proposed for this mechanical property change [36,37,56±63]: the chain slippage between the neighboring crystallites; the anharmonic chain motion in the crystallites; and so on. The present MD calculation is considered to give a concrete and quantitative image to the anharmonic chain motion as a possible origin of the a relaxation mechanism. A similar concept was reported by Mori et al. [64], although their model was not so large as ours, having a possibility to be affected by the periodic boundary conditions as discussed in the previous
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
371
Fig. 16. Temperature dependencies of: (a) Young's modulus along the chain axis; (b) averaged torsional angle of the skeletal chain; (c) the ¯uctuation of the torsional angle from the averaged value; and (d) the chain contraction calculated for orthorhombic PE crystal.
section, and the reasonableness of the used potential functions were not checked. Of course, we do not consider that the present calculation could give the correct and ®nal answer to the a relaxation process of PE. In order to clarify the real image of this process, we are now calculating the relaxation time for the thermal motions of ¯ip-¯op, translation, and skeletal twisting of the PE chains, which might be compared with the observed relaxation time. 5. Conclusions
Fig. 17. Correlation between Young's modulus and the contraction of the molecular chain calculated for orthorhombic PE crystal at the various temperatures (see Fig. 16).
In the present paper we have mainly described the results of the MD calculations made for the orthorhombic PE crystal model consisting of 5184 atoms under the threedimensional periodic boundary condition. The cell constants were found to correspond well with the actually observed data. In particular, the linear thermal expansion coef®cients
372
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
Fig. 18. Molecular conformations in trigonal POM cystal calculated at the various temperatures.
were in good agreement with the experimental values. The molecular chains were found to show the ¯ip-¯op-type motion at high temperatures, and at the same time the torsional ¯uctuation of the skeletal chains became large. These thermal motions re¯ect on the large thermal parameters of the H and C atoms at high temperature. The thermal parameters show the de¯ection point in the vicinity of 350 K, corresponding relatively well to the X-ray analysis result. Young's modulus along the chain axis was evaluated by the second derivatives method of the potential function. The modulus at room temperature, ca. 280 GPa, is in good agreement with the experimental value 280±290 GPa obtained by the LAM measurement of the n-alkane crystals. Young's modulus was found to closely relate with the chain contraction, which was caused by the thermal ¯uctuation of the skeletal chain conformation. Young's modulus was found to decrease by 35% for quite a small contraction of the planar-zigzag chain by 0.3%. A similar situation was also found in the case of nylon 6 and trigonal POM. This way, Young's modulus along the chain axis is affected sensitively by the slight change in the chain conformation. Therefore, if we plan to produce a new polymer material with high thermal stability of the mechanical properties, we will still have to ®nd a way to ®x the chain length at higher temperatures. Acknowledgements The author wishes to thank Dr Tahir Cagin of Molecular Simulations Incorporation for his kind discussion about the numerical calculation of the elastic constants.
Fig. 19. Temperature dependence of Young's modulus along the chain axis and the c-axial length calculated for trigonal POM crystal. The open circles in the top ®gure are the observed data by the X-ray method. The bottom ®gure shows the correlation between Young's modulus and the chain contraction extracted from the above two ®gures.
functions and the re®ned parameters are listed. The total energy of the model is expressed as a summation of the following potential functions. EBond Stretch exp 2ab R 2 Rb 2 12 ; kb 2Db a2b EAngle Bend Ccos u 2 cos ua 2 =2; ku C sin2 ua ETorsion Vt 1 1 cos 3f=2
Appendix A. Potential functions used in the present MD calculation and the numerical values of the related parameters
EBond-Angle D cos u 2 cos ua R 2 Rb ; kru 2D sin ua
As stated in the text the potential functions proposed by Karasawa et al. [27] were re®ned so that the vibrational spectroscopic data could be reproduced as well as possible. In this appendix the concrete forms of these potential
EOne-Center Angle±Angle Cross-Terms
EBond±Bond krr R1 2 R1b R2 2 R2b
G cos uIJK 2 cos uaIJK cos uIJL 2 cos uaIJL
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
ETwo-Center Angle±Angle Cross-Terms F cos f cos uIJK 2 cos uaIJK cos uJKL 2 cos uaJKL Eelectrostatic Qi Qj = 10 1RIJ EvdW A exp 2BRij 2 C=R6ij DV 6 exp z 1 2 r 2 zr26 = z 2 6 van der Waals parameters for C and H atoms Karasawa (MCXX) Tashiro zC 13.0 13.0 3.8410 3.8410 RVC 0.07918 0.0792 DVC zH 11.2 11.2 3.1665 3.1665 RVH 0.0200 0.0190 DVH Valence force ®eld parameters for C and H atoms Bond stretch C±C
Rb kb Db Rb kb Db
1.4841 884.9940 85.80 1.0765 741.3720 95.10
1.5300 884.9940 85.80 1.0765 580.0000 90.00
ku ua ku ua ku ua
55.6076 119.3933 65.7301 117.7291 84.1810 121.2400
48.0000 119.3933 72.0000 117.7291 84.1810 121.2400
Vt Vt Vt
5.1686 6.1626 5.7070
5.1686 6.1626 5.7070
222.6583 3.1321 234.3195 225.9234 1.3684 254.0185 26.2187
222.6583 3.1321 234.3195 225.9234 1.3684 254.0185 26.2187
One-center angle±angle cross-terms GCC:HC GCC:HH GCH:CC GCH:CH
27.6083 25.3356 25.0824 25.3356
27.6083 25.3356 25.0824 25.3356
Two-center angle±angle cross-terms FH:C:C:H FC:C:C:H FC:C:C:C
217.7274 216.4004 221.5910
217.7274 216.4004 221.5910
C±H
Angle bend H±C±H C±C±H C±C±C Torsion H±C±C±H C±C±C±H C±C±C±C
Angle cross-terms H±C±H DH u kHH C±C±H DC u DH u kHC C±C±C DC u kCC
References [1] Tashiro K. Kobunshi Ronbunshu (Jpn J Polym Sci Technol) 1992;49:711. [2] Tashiro K. Prog Polym Sci 1993;18:377. [3] Tashiro K, Kobayashi M. Polymer 1995;37:1775. [4] Jungnitz S, Jakeways R, Ward IM. Polymer 1986;27:1651.
373
[5] Matsuo M, Sawatari C. Macromolecules 1988;21:1653. [6] Nishino T, Ohkubo H, Nakamae K. J Macromol Sci Phys B 1992;31:191. [7] Wu G, Tashiro K, Kobayashi M, Komatsu T, Nakagawa K. Macromolecules 1989;22:758. [8] Wu G, Tashiro K, Kobayashi M. Macromolecules 1989;22:188. [9] Parrinello M, Rahman A. J Chem Phys 1982;76:2662. [10] Ray JR. Comp Phys Rep 1988;8:109. [11] Tashiro K, Kobayashi M, Tadokoro H. Polym J 1992;24:899. [12] Cagin T, Kawasawa N, Dass Gupta S, Goddard III WA. In: Mark JE, editor. Computational methods in materials science, 1992. [13] Kobayashi M. J Chem Phys 1979;70:509. [14] Lacks DJ, Rutledge GC. J Phys Chem 1994;98:1222. [15] Cagin T, Karasawa N, Dasgupta S, Goddard III WA. Mater Res Soc Symp Proc 1992;278:61. [16] Palmo K, Krimm S. J Polym Sci: Polym Phys Ed 1996;34:37. [17] Strobl GR, Eckel R. J Polym Sci: Polym Phys Ed 1976;14:913. [18] Kobayashi M, Sakagami K, Tadokoro H. J Chem Phys 1983;78:6391. [19] Chang C, Krimm S. J Polym Sci: Polym Phys Ed 1979;17:2163. [20] Snyder RG, Strauss HL, Alamo R, Mandelkern L. J Chem Phys 1994;100:5422. [21] Hsu SL, Ford GW, Krimm S. J Polym Sci: Polym Phys Ed 1977;15:1769. [22] Tashiro K, Ishino K, Ohta T. Polymer 1999;40:3469. [23] Kitagawa T, Miyazawa T. Adv Polym Sci 1972;9:335. [24] Weiner JH. Statistical mechanics of elasticity. New York: Wiley, 1983. [25] Kobayashi M, Tadokoro H. J Chem Phys 1977;66:1258. [26] Kobayashi M. J Chem Phys 1979;70:4797. [27] Karasawa N, Dasgupta S, Goddard III WA. J Phys Chem 1991;95:2260. [28] Palmo K, Mirkin NG, Pietila LO, Krimm S. Macromolecules 1993;26:6831. [29] Tashiro K, Kobayashi M. Rep Prog Polym Phys Jpn 1996;39:349. [30] Sun H. J Phys Chem B 1998;102:7338. [31] Tashiro K, Sato M, Cho Y, Kobayashi M. Polym Prepr Jpn 1999;48:863. [32] Tashiro K, Sato M, Cho Y, Kobayashi M. Polym Prepr Jpn 1999;48:3574. [33] Tashiro K, Satoh M, Ike M. Polym Prepr Jpn 2000;49:2541. [34] Gasteiger J, Marsili M. Tetrahedron 1980;36:3219. [35] Pennings AJ, Zwijnenburg A. J Polym Sci: Polym Phys Ed 1979;17:1011. [36] Syi JL, Mans®eld ML. Polymer 1988;29:987. [37] Mans®eld ML, Boyd RH. J Polym Sci: Polym Phys Ed 1978;16:1227. [38] Tadokoro H. Fiber diffraction methods. In: French AD, Gardner KCH, editors. ACS Symposium Series 141Washington, DC: American Chemical Society, 1980. [39] Sumpetr BG, Noid DW, Wunderlich B. J Chem Phys 1990;93:6875. [40] Noid DW, Sumpter BG, Wunderlich B. Macromolecules 1990;23:664. [41] Liang GL, Noid DW, Sumpter BG, Wunderlich B. J Polym Sci: Polym Phys Ed 1993;31:1909. [42] Liang GL, Noid DW, Sumpter BG, Wunderlich B. Acta Polym 1993;44:219. [43] Liang GL, Noid DW, Sumpter BG, Wunderlich B. J Phys Chem 1994;98:11 739. [44] Tuzun RE, Noid DW, Sumpter BG. Makromol Chem: Theory Simul 1998;7:203. [45] Kobayashi M, Tadokoro H. Macromolecules 1975;8:897. [46] Tashiro K, Sasaki S, Kobayashi M. Macromolecules 1996;29:7460. [47] Nakamae K, Nishino T, Miyazaki H. Polym Prepr Jpn 1999;48:914. [48] Barham PJ, Keller A. J Polym Sci: Polym Lett Ed 1979;17:591. [49] Holliday L, White JW. Pure Appl Chem 1971;26:545. [50] Tashiro K, Wu G, Kobayashi M. Polymer 1988;29:1768. [51] Tashiro K, Tadokoro H. Macromolecules 1981;14:781. [52] Miyasaka K, Isomoto T, Koganeya K, Ishkawa K. J Polym Sci: Polym Phys Ed 1980;18:1047.
374
K. Tashiro / Computational and Theoretical Polymer Science 11 (2001) 357±374
[53] Nakamae K, Nishino T, Hata K, Matsumoto T. Koubunshi Ronbunshu (Jpn J Polym Sci Technol) 1987;44:421. [54] Nakamae K, Nishino T, Shimizu Y, Hatano K. Polymer 1990;31:1909. [55] Ohta Y, Yasuda H. J Polym Sci B: Polym Phys 1994;32:2241. [56] Kajiyama T, Okada T, Sakoda A, Takayanagi M. J Macromol Sci Phys B 1973;7:583. [57] Kyu T, Yasuda N, Suehiro S. Polym J 1976;8:565. [58] Suehiro S, Yamada T, Inagaki H, Kyu T, Nomura S, Kawai H. J Polym Sci: Polym Phys Ed 1979;17:763.
[59] Suehiro S, Yamada T, Kyu T, Fujita K, Hashimoto T, Kawai H. Polym Engng Sci J 1979;19:929. [60] Suehiro S, Kyu T, Fujita K, Kawai H. Polym J 1979;11:331. [61] Ogita T, Yamomoto R, Suzuki N, Ozaki F, Matsuo M. Polymer 1991;32:822. [62] Usami T, Takayama S. Macromolecules 1984;17:1756. [63] Matsuo M, Sawatari C, Ohhata T. Macromolecules 1988;21:1317. [64] Mori S, Ohta Y, Aikawa Y. Polym Prepr Jpn 1994;43:3419. [65] Avitabile G, Napolitano R, Pirozzi B, Rouse KD, Thomas MW, Willis BTM. J Polym Sci B 1975;13:351.