A Monte Carlo study of the glass transition in three-dimensional polymer melts

JOURNAL OF

ELSEVIER

Journal of Non-Crystalline Solids 201 (1996) 199-210

A Monte Carlo study of the glass transition in three-dimensional polymer melts M. Wittkop *, Th. H~51zl, S. Kreitmeier, D. G~Sritz Unit,ersittt Regensburg, Institut ftir Experimentelle und Angewandte Phvsik D-93040 Regensburg, German3' Received 19 October 1995; revised 13 December 1995

Abstract

The glass transition of densely packed linear polymer chains has been investigated by means of the three-dimensional bond-fluctuation model including van der Waals interaction and bond-length potential. The computer simulations presented here are based on different chain-lengths, cooling rates, and a wide range of temperatures. All monitored static and dynamic properties of the polymer melt change distinctly at one characteristic temperature, identified as Tg. Variations of the chain-length from 50 up to 200 monomers per chain yield nearly identical behavior, while an increase of the cooling rate from FQ = 2.7 × 10 -8 MCS- t (Monte Carlo steps) to infinity causes pronounced differences below Tg both in static and dynamic properties.

1. Introduction

In the process of rapid cooling of liquids or melts, nearly every material vitrifies continuously into an amorphous solid, conventionally called 'glass' [1]. The process of vitrification occurs in a narrow temperature interval near Tg, the glass transition temperature, depending on the nature of the liquid and the cooling rate. Polymer melts especially have a marked tendency to vitreous solidification, even at moderate cooling rates. Due to the enormous scope of application for polymer glasses, manifold efforts of experimental and theoretical research have been made over the last

* Corresponding author. Tel.: +49-941 943 3194; fax: +49-941 943 3196; e-mail: [email protected].

decades in order to study the nature of the glass transition. The results of this work have been documented in numerous publications [2-13]. More recently, the dynamics of structural rearrangement were studied by quasielastic neutron scattering (QUENS) [14]. An overview can be found in the article of Gentile and Suter [15]. Computer simulations have impressively demonstrated their possibilities in investigating liquids, the glass transition and the structure of glasses. The static and dynamic properties of liquids were studied by various authors both with Monte Carlo techniques [16-31] and molecular dynamics [32-34]. Molecular mechanics, sometimes combined with molecular dynamics, was often used to investigate atomistic models and detailed molecular structures of special polymer glasses: for polypropylene see Theodorou and Suter [35] as well as Mott, Argon and Suter [36], for

0022-3093/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved. PII S 0 0 2 2 - 3 0 9 3 ( 9 6 ) 0 0 14 1 -X

200

M. Wittkop et al./ Journal of Non-Crystalline Solids 201 (1996) 199-210

polymethylene refer to Boyd and Pant [37] and for aromatic polysulfone to Fan and Hsu [38]. Polycarbonate was investigated by Fan et al. [39] as well as Hutnik and co-workers [40]. The statics and dynamics of glasses and structural rearrangement during the glass transition were the scope of the molecular dynamics simulations by Roe and co-workers [4143]. Monte Carlo simulations were also used to investigate the glass transition [44-52]. Since glass transition is recognized as a universal phenomenon [1], chemical details of specific polymers are negligible in obtaining the common features of the transition. In the sense of such coarse-grained models, the glass transition was simulated via a geometrical frustration process by Wittmann, Baschnagel, Ray, Lobe and Binder [45-52]. These simulations were successful in describing cooling rate and chain-length dependencies of the glass transition of short chains. They were also promisingly interpreted within the framework of the idealized mode coupling theory (MCT) [53]. Despite all these efforts, the problem of a strictly theoretical description is not solved. As a consequence, additional computer simulations must be regarded as an appropriate means to investigate these complicated interrelations. Since molecular dynamics suffers from its limitation to short time regime, and, since we are not interested in chemical details of a particular polymer, we took advantage of the efficiency of a coarse-grained lattice model which enabled us to obtain larger time scales. Using the bond-fluctuation model we intend to simulate the glass transition of systems consisting of long chains, including van der Waals interaction and bond-length potential. This is in contrast to the work of Binder and co-workers [45-52] who studied glasses with short chains ( N < 15) without intermolecular interaction. Due to the short chain length, chain-end effects are pronounced. The glassy state in their work is obtained by an artificial geometric frustration process introduced via the bond-length potential. In our simulations we will take a different approach. Our chains are long (49 < N < 199) and chain end effects should be negligible. In case of the longest chain length the initial liquids will access the 'reptation' regime. Also the glass transition should be initiated by the van der Waals interaction and not by a geometric frustration. Moreover the usage of a van

der Waals potential seems to be reasonable since there exist hints that non-bonded interactions are not completely screened even in liquids [54]. As our next goal is to study uniaxial deformation we want to generate samples which - - after cooling to below the glass transition - - should keep its form even when enlarging the simulation box. This process will be managed by a proper choice of the van der Waals interaction. To prevent the system from entirely adopting the long range order of the underlying lattice, we introduced a bond-length potential. In the next section the simulation method will be described. The presentation of the results and a discussion will follow.

2. Simulation method In order to consider intra- as well as intermolecular interactions, we combined the well-known threedimensional bond-fluctuation model [45,55-57] with a van der Waals potential and a suitable bond-length potential. In this model the polymer chains are created by mutually- and self-avoiding walks on a cubic lattice, each monomer occupying eight comers of the unit cell (for details see Refs. [55-57]). A major advantage of the bond-fluctuation model in contrast to simple cubic lattice models is its flexibility which leads to a relatively high acceptance rate even at low temperatures and high densities. In addition to the excluded-volume interaction we take into account a truncated Lennard-Jones potential, V(r), with a range of three grid units between all points of the monomers ( r is the distance between the monomers) and a bond-length potential, V'(I), describing the increasing backbone stiffness with decreasing temperature:

= E{{r)-12 2( /,.) -6} V(r)

v'(t)=

~

ro

c0+c a-+c2

-~o

' +c3

(1) ,

(2)

where ~ is the unit of energy, r 0 = 1.1 a and a is the lattice spacing. The choice of the cut-off parameter ensures both the attractive and repulsive forces. The form of the van der Waals equation assures that the minimum is at the point r 0. The parameters, c o =

M. Wittkop et al. / Journal of Non-Crystalline Solids 201 (1996) 199-210

-207.12, c t = 342.88, c 2 = - 163.52, c 3 = 24.32, were estimated from a static energy-minimization simulation, including 10 unified PE-units interacting by bond stretching-, torsional- and van der Waals forces. The identification of a bond with a Kuhn segment of 10 PE-segments is somewhat arbitrary, but is in approximate agreement with values used in literature [46,58,59]. Note that the exact form of V'(l) has no influence on the qualitative results. Hence V'(l) is equivalent to the potential used for modeling chain stiffness in Refs. [45-52]. With these potentials, the jumps of monomers in our model were accepted using the Metropolis criterion [60]. As a consequence of the model hamiltonian, a competition arises between the attractive V(r) term and the bond energy, V'(l). The dominating attractive van der Waals interaction, in conjunction with non-periodic boundary conditions, results in a compact sample below Tg. Below Tg the sample coheres to a self-supported cluster which does not fill the complete box anymore. The non-periodic boundary conditions were modeled by hard walls. The unusual choice of boundary conditions arises from the intention to create samples suitable for uniaxial deformation. But the non-periodic boundary conditions yield another advantage. Variations of volume with changing temperature are permitted. The initial configuration was randomly created at infinite temperature. After sufficient equilibration the continuous cooling-process starts with cooling rate, F o. F o denotes the reciprocal temperature step per time unit of the Monte Carlo simulation. The inverse temperature, 13 = 1/kaT, at time, t, follows from 13e = FQt. The simulated systems contain P chains with N segments in a cubic box of linear dimension L = 90. All parameters are listed in Table 1. Therefore, every simulation box contains 40000 monomers, corresponding to an initial density of

0.44 at infinite temperature. It is known that, at this density, the bond-fluctuation model should represent the properties of a liquid [29,46]. During cooling, configurations were extracted at intervals of A13e = 0.05 and afterwards subjected to isothermal annealing for a duration of 1.1 × 1 0 7 timesteps. Time will be measured in units of MCS (Monte Carlo steps), denoting one attempted move per monomer. The last 5 )< l 0 6 MCS of the annealing process were used for the determination of time-dependent and time-averaged quantities.

3. R e s u l t s

N + 1 P F . ( M C S -~ )

50 800 _ 2.7)<10 -7

100 400 2.7)<10 -7

200 200 2.7×10 -8 2.7>(10 -7

and discussion

Before starting the cooling process we tested whether our athermal systems resemble melts. For this purpose we analyzed the athermal segment diffusion, more precisely the mean square displacement, g~(t), of the monomers N/2 + 2

g,(t) = ½

E

((Ti(t) - r~(O))2),

(3)

i= N/2 - 2

~(t) being the position of monomer i at time t. The brackets ( ) denote an average over all P chains in the sample. Only the five innermost monomers of a chain are summed to reduce possible chain-end effects [22,32,33]. As long as the mean square displacement is much smaller than the box size the problems due to the non-periodic boundary conditions can be neglected. Furthermore we made simulations with periodic boundary conditions and found the same behavior. For short chains in a melt g~(t) should obey the following power laws according to the Rouse model [61,62]

g,(t) ~ Table 1 Parameters o f the different samples

201

t I,

t<

t I/2 ,

"Co < t < ~'N

t 1,

t>rN.

TO

(4)

For times shorter than %, the monomer diffuses freely and -rN ~ N 2 is the largest relaxation time of the chain (the Rouse relaxation time). However for large chains the reptation concept [63,64] modifies this behavior. For chains which significantly exceed

M. Wittkop et al. / Journal of Non-Co'stalline Solids 201 (1996) 199-210

202

s e g m e n t diffusion: 13e=0 I

I

I

I

I

[

103

I

,"

...'"" , f * .,.'" J ,,: f . t

..*"

ii

5 O

"~

I

• N---49 P=800 • N--99 P--400 • N=199 P=200

2

e ~ 102

-'"

I

2

104

5

t

I

I

2

5

I

105

I

I

I

I

2

5

106

2

5

t (units o f M C S ) Fig. I. Segment diffusion, gl(t), for chain lengths N = 49, 99 and 199.

a critical entanglement length Are, gl(t) should be described by the following power-law sequence: Ct 1 '

gl(t)'~

t < %,

tl/2 '

70_
2,

t I/4,

Te
~,

/1/2,

Tu < t < rd "~ N 3 / N e ,

t I,

~'d <__t.

(5)

In the new time-range, Te < t < "rd, the chain should exhibit a Rouse diffusion along its coarse-grained random walk structure. The chain is confined in a so-called tube built by the surrounding chains. In the bond-fluctuation model Ne is 45 for a density of about 0.44 [29]. Plots of gt(t) for N = 49, 99 and 199 are shown in Fig. 1. Since the initial ~ t behavior is hard to observe, this time-range is omitted in the figures. In the case N = 49 the ~ tl/2-regime can be observed for a long time-range without any slowing down. Hence the diffusion follows the Rouse prediction. This was expected since these chains do not significantly exceed Ne. For larger times a broad crossover to ~ t takes place. In the case of N = 99 a slowing down of the diffusion for intermediate times can be

seen. This regime is more pronounced in the case of N = 199. All in all the diffusion results confirm that our initial systems are liquids and show some 'reptation' behavior for the longest chains. This behavior can be seen in Fig. 1 where the diffusion slows down to a power law with exponent ~ 1 / 3 in an intermediate time regime ~.

3.1. Static properties After we have established that the initial state behaves like a melt, we intend to investigate the static properties obtained by cooling these melts.

3.1.1. Free volume One of the theories concerning glass transition tries to derive the molecule's diffusion coefficient from the available free volume. Following this approach, we calculate the average free volume of a monomer, < v). According to Wittmann et al. [45], ( v > was evaluated by checking the possibilities of

i For those who are interested in the discussion about reptation, we do not find a time regime with an exponent of 1/4. But this is a common finding in many simulations, see, for example, [29] and references therein.

M. Wittkop et al. / Journal of Non-Crystalline Solids 201 (1996) 199-210

203

free volume vs. [3

~

.

'

'

'

"":~'L:

'

I

.

.

.

.

I

.

.

.

.

I

0.05

li:,

' -',::.

0.15

'

'

; ....

'

'

~ ....

I

, ....

~;;;'~

"":z

....

.

-

. ...

0.035 ~;; il,

. . . . . A. . . . . . .

0,03 . . . . 0.3

~ .... 0.35

\k

x N---49, FQ=2.7 10 "7 o N - - ~ , FQ=2.7 10 .7

' .... 0.4

• "7

v N = 1 9 9 , FQ=2.7 10 .7

0

i

i

I

0.1

i

0.5

-..

. . . . . . . . . . . . i

i

i

I

0.2

i

t

i

i

. - C

~k.z-.-.';.-.~a . . . . . . . g . . . . . . . o . . . . .

a N = 1 9 9 , FQ=2.7 10 "s i

' .... 0.45

"-'%,;

o N = 1 9 9 , FQ= oo

i

~ .... ...o'"

"O

0.04

0.1

'

...-

0.045 - ' : -

A V

0.05

'

~:,., '-

0

'

'

I

i

i

0.3

i

i

I

A . . . . . . . . i

i

0.4

i

A

i

0.5

inverse temperature Be Fig. 2. Free v o l u m e , ( v ) , v e r s u s g e for all chain lengths and c o o l i n g rates.

the monomers to jump in an arbitrary spatial direction disregarding the Metropolis criterion. ( v ) is normalized to 0 < ( v ) < 1 and was averaged over all monomers in a sample. In Fig. 2 ( v ) is plotted versus the inverse temperature. The included error bars represent the time averaged standard deviation. In most cases the error is about the symbol size, only for temperatures slightly above Tg larger fluctuations can be observed. ( v ) shows a distinct change around 13~ = 0.25. Thus 13e = 0.25 can be used for an estimate of the glass transition temperature Tg. A precise determination of Tg would require the finding of an intersection point of two straight lines, which are extrapolated from the glassy, respectively, liquid region [1,4]. But in the case of 'S-shaped-curves' like Fig. 2 or the following figures of the static and dynamic properties, this can not be performed unambiguously [46,51]. However a precise determination of Tg is not the aim of this paper. Varying the chain length from 49 to 199 segments has no significant influence on the behavior 2. Since

2 ( c ) is slightly larger for shorter chains. This is d u e to the fact that short c h a i n s are m o r e m o b i l e .

our chains are relatively long, this result is in agreement with experiments [66] and simulations [50,51] which state, that for sufficiently long chains Tg is nearly independent of the chain length. For example in the simulation of Lobe et al. [50] Tg becomes nearly constant for N>_ 12 (Fig. 11 in Ref. [50]). Changing the cooling rate from Fo = 2.7 x 10 -8 MCS -~ to 2.7 x 10 - 7 M C S - j and to the limit FQ = ~ results in a distinct increase of differences for inverse temperatures above 13e = 0.25. The more slowly cooled sample reveals the lowest free volume. Cooling below Tg, the system leaves equilibrium and gets locked in a structure which depends on the cooling rate. These structural differences can not vanish completely during the isothermal annealing process. 3.1.2. Energy

While cooling the sample, the Metropolis procedure [60] decreases the energy of the system. For the systems under consideration, the temperature dependence of the mean potential energy per monomer ( E ) exhibits (like all other quantities) the most pronounced change at 13~ = 0.25 (see Fig. 3). As expected from the free volume, the slowest cooled sample is energetically best saturated. The error bars

204

M. Wittkop et aL / Journal of Non-Cry'stalline Solids 201 (1996) 199-210

potential energy vs. I

.

.

.

.

I

.

.

.

.

I

'

-10.5

"11

. . . . I ;"..

-5

I

. . . .

I

0

r

I

. . . . .°., .''

. O "

I~:..".

-11

. . . .

...0

"-',

.-°

"o- ....

-11.5

"'"-.A..

\

-7.5

-12 . . . . 0.3

A

J .... 0.35

, , . . , i , , , , 0.4 0.45 0.5

V x

-10

N--49, FQ=2.7 10 -7

o N--99, FQ=2.7 10 -7 o

"":;" •.* . . . . ~-.'-,?,; . . . . . . . . . . -o . . . . . .

N=199, FQ= ,,~

v N=199, FQ=2.7 10 "7 a N=199, FQ=2.7 10 "8 e

n

e

0

e

I

......

e

e

I

0.1

I

e

I

J

J

0.2

I

n

e

e

-o . . . . A

A ........ e

I

0.3

I

....c .

........

e

a

A " e

0.4

0.5

inverse temperature ~e Fig. 3. Potential e n e r g y , ( E ) / e ,

v e r s u s 13e for all chain lengths and c o o l i n g rates.

are determined in the same way as pointed out in the Section 3.1.1. Analyzing the contributions of the van der Waals

potential, V(r), and the bond-length potential, V'(I), we find that the pronounced change in Fig. 3 results from the van der Waals interaction. The mean bond2

radius of gyration /N vs. .

.

.

.

I

.

...

.

-""

.

.

I

"~ ........

.

.

.

.

I

.

.

.

.

I

.

.

.

.

"X

1.6 .o|..-°" o.O.--

[] .......

u.

....... ~)....... o

~

',

1.5 • o~

O

x.J

. . oOV.

.

1.4 .....

0

~,~...::::::~-2 . . . . . x A "

V

1.3

',

'

,'"

. .El- .......

[] ........

ra

x N=49, FQ=2.7 10 "7 o N=99, FQ=2.7 10 -7

" " ",.[]"

o N=199, I"Q= ,~

\,

i',

,"

,

,V"

....re.

.......

,q- .......

.

v N=199, FQ=2.7 10 "7

1.2

a N=199, FQ=2.7 10 "s .

0

.

.

.

I

0.1

i

0"" . . . . . ~" . . . . . . . /k- . . . . . . . A. . . . . . . . A. . . . . . . . ,fi i

i

I

0.2

i

e

i

e

I

e

0.3

e

a

e

I

,

,

,

0.4

inverse temperature ~e Fig. 4. R a d i u s o f gyration, ( R~ } / N , v e r s u s 13e for all chain lengths and c o o l i n g rates.

,

I

0.5

M. Wittkopet al./Journal of Non-CrystallineSolids 201 (1996) 199-210 length energy gradually decreases with decreasing temperature by a total amount of about 2.5e in comparison to total loss of about 7.5e. Hence not only the shape of the energy behavior in Fig. 3 but also the major contribution to the total energy change results from the van der Waals interaction.

3.1.3. Radius of gyration In terms of structural properties, one can also see the characteristic glass transition. This can be obtained by plotting the normalized mean square radius of gyration, (R~)/N, averaged over all chains (see Fig. 4). ( R ~ ) versus 13e shows distinctly smaller values in the glass regime. As a hint for the possible uncertainties we included a representative error bar for each chain length. These error bars were determined from the calculation of the standard deviation of R2g/N for the final configuration, averaged over all chains in the system. As the monomer number stays constant the number of chains increases with decreasing chain length. This leads to a better averaging for the shorter chains, being expressed by the smaller error bars. When determining the error from the time evolution the error bars are more or less the same in the liquid regime, but are diminishing considerably for temperatures below Tg. This change is not surprising since the dynamics of the glassy system is frozen. In general the following two points should be considered. First, the independence of the configurations for the averaging process must be determined. Unfortunately this is a rather vague guess for configurations around and below T~. The occurrence of large correlations in time is one of the most important properties of a glass transition. A second point is that for a more quantitative study more independent systems should be simulated. But considering that the presented data are based on about 14000 CRAY-equivalent hours of CPU-time this is not feasible up to now. Variation of the cooling rate for the system with chain length 199 reveals an interesting aspect: the higher the cooling rate, the more pronounced the hindrance of global relaxation, especially for temperatures far below Tg. Only in the vicinity of the glass transition temperature (and naturally above) structural differences - - as far as they are interrelated with (R2~) - - are completely reduced by isothermal

205

annealing. Hence a minimum of (R2g) arises at Tg which becomes more pronounced with higher cooling rates. Below Tg the sample is locked in a nonequilibrium structure. It should be mentioned, however, that the sample with chain-length 199 is affected by a finite-size effect, which can be recognized in the non-compliance of the scaling law, ( R ~ ) ~ N, for 13e = 0.0. Due to the screening of excluded-volume in dense systems the chains should behave like non-reversal random walks (( R2g ~ N) in contrast to single chains which behave as self-avoiding walks ((R~ ~ N~18). In Fig. 4 ( R ~ ) / N nearly coincides for N = 49 and N = 99 (for 13e = 0.0) but the values for N = 199 are smaller, As our samples are clearly relaxed in the athermal case, the relaxation period of 6000000 MCS is much larger than the typical relaxation time of 10000 to 100000 MCS 3 this difference indicates, that the size of the simulation box should have been larger. But considering the correctness of the diffusion behavior, which will be discussed next, this point should not be overestimated. Finally, note that the same characteristics were found for the mean square end-to-end distance.

3.2. Dynamic properties Let us now turn to the dynamic behavior and determine whether there are similar changes as for the static properties.

3.2.1. Diffusion properties In order to study the temperature dependence of g i(t) we plotted the mean square displacement, g i(t) in Fig. 5 for various temperatures for the case of N = 199 and FQ = 2 . 7 × 10 -7 MCS -1. The errors are in the range of the symbol size. The curves for the other N and FQ values obey the same qualitative behavior, The mean square displacement, g~(t), shows a significant reduction of segment diffusion for lower temperatures and a maximum change at 13e = 0.25. This is obviously an onset to freezing.

See for example Fig. 6. where all bond correlations are lost within this time frame.

206

M. Wittkop et al. / Journal of Non-Crystalline Solids 201 (1996) 199-210

segment diffusion: N=199, P=200, FQ=2.7 10-7 I

I

I

I

I

+ ~=o.oo ® 1~=o.o5

102

~*~ 101

-~'--.-'~

B ~e---0.15

~

o 13~--o.2o * 1~=0.25 lie--0.30 13~=0.35 o 13~=0.40

/ J ~ . / -

=

/

/

/ -

15~=0.45

3 =7

lO0

10 -1 101

102

103

104 t

105

106

(units of MCS)

Fig. 5. Temperature dependence of the segment diffusion gl(t) for N = 199, cooling rate FQ = 2.7 × 10 -7 M C S - t and inverse temperature varying from 13~ = 0 to 13e = 0.50.

bond-autocorrelation: N=199, P=200, .1"°=217" -10-7 5

1002

+®~--0.0013E=0.0'5

DI~--~'2:

~~ ~ :~ ~/:, "~"~ ~' :..~ ~ ' g ~ ~":' ..~ . " . ~ . . . . . .

/ / /

/

/

- . . . . . . . .~ -

/

-

2 10-1

5

- - - -.

13e.---0230

//fjf

v t

100

-

101

I

I

102

10 3 t

I

104

I

10 5

- - - KWW-fit I

106

(units of MCS)

Fig. 6. Bond-autocorrelation function, - I n qb(t), versus t for N = 199, FQ = 2.7 × 10 -7 M C S - i and inverse temperature varying from 13e = 0 to [3~ = 0.50. Simulation data (solid lines) are compared with the KWW-fit formula ~ ( t ) = exp - (t/T) y (dashed lines).

M. Wittkop et al. / Journal of Non-Crystalline Solids 201 (1996) 199-210

207

exponent y vs. 13 0.4

t

I{ i

. . . .

.......

0.35

0.2 I I .

.

.

,

'

'

,

'

'

. . . .

,

,

,

0.175

.

"""~,

0.3

"

o 0.25 O

0.15

-.~

...

-.:-xt. " - - . ~ .

..m.:...

-..

0.125 ."

,

0.I

"

,

0.25

,

l

. . . .

0.3

I

,

,

i

0.35

•

I

. . . .

0.4

I

i

i

I

0.45

0.5

0.2 x N=49, FQ=2.7 10"7

i~'.

"":::~::

u N=99, FQ=2.7 10.7 0.15

0.1

o N=199, FQ. . . v N=199, FQ=2.7 10.7 /x N=199, rQ=2.7 10 .8 . . . . . . . . 0 0.I

.

.

.

..... i';::--~.'.:::l'_"--..{~..~::

.

.... I . : : - ~ ' I ....... .......I

I . . . . 0.2

xX. . . . . . . . t . . . . 0.4 0.5

I . . . . 0.3

inverse temperature 13e Fig. 7. Plot of the exponent, y, versus ~ ( as derived from the KWW-fit formula for all chain lengths and cooling rates.

~'VS. .

1013

.

.

.

I

'

'

1011

x N---49, FQ=2.7 10.7 [] N=99, FQ=2.7 10 "7 o N=199, FQ= oo v N=199, FQ=2.7 10 "7

1010

A N=199, 1"Q=2.7 104

1012

'

'

I

.

.

.

.

I

T

,

,

,

I

.

.

.

.

I

.{. ,{-""

..~

.,"" ,.~'"

..>.~'"

.... ~-

t . . . . 0.4

0.5

~.) 109 10 g

. ::;;:F: ....

O 10 7 "~10 6

•

,::.4};:::

o•

105

104

,~""

103 102 101 100

. . . . ,i~' --'• ....... ~ ....... ~ . . . . 0

,

0.1

4]""

. . . .

,

0.2

. . . .

,

. . . .

0.3

,

inverse temperature 13e Fig. 8. Plot of the "r versus 13¢ as derived from the KWW-fit formula for all chain lengths and cooling rates.

M. Wittkop et al. / Journal of Non-Crystalline Solids 201 (1996) 199-210

208

3.2.2. Relaxation properties Generalizing the idea of the Edwards-Andersonorder parameter for dense polymer liquids according to Wittmann et al. [45], the suitable standardized bond-autocorrelation function, qb(t), may be employed for a convenient description of relaxation 4 processes :

qb(t) = - ~

~_~ 6[bi(O), bi(t)] - - ~

,

i=1

with

6[x, y]={O' 1,

x~y, x=y,

(6)

with M = PN and bi(t) denoting the ith bond at timestep t. qb(t) is normalized between 0 and 1. In an intermediate time regime, where the relaxation of the monomer as a part of the polymer chain can be monitored, it is possible to fit the decay of correlations according to the K o h l r a u s c h Williams-Watts [65] formula by a stretched exponential function (see Fig. 6): ~(t) = exp[--(t/T)Y].

(7)

Since in the studied system several relaxation processes are superimposed the Kohlrausch-WilliamsWatts formula fits the data only for a certain time window. The curves for the other N and FQ values obey the same qualitative behavior as shown in Fig. 6. The determined relaxation time "r as well as the exponent y show a distinct change around [3¢ = 0.25 (see Figs. 7 and 8) with significant cooling rate dependence for larger inverse temperatures 5. This change indicates a glass transition near [3¢ = 0.25 in agreement with the results of the static properties. Fig. 8 shows that -r increases by 11 decades when the temperature is decreased from 13¢ < 0.25 to 13¢ > 0.25. This increase demonstrates the freezing of the melt at Tg and the sharp rise of relaxation times when cooling below Tg. The distinct changes due to the different cooling rate below the glass transition temperature are in accordance with the idea that highly quenched sam-

4 The n u m b e r 108 arises f r o m the set of permitted bonds. 5 The included error bars were determined b y a variation o f the fit parameters within reasonable limits.

pies are far from equilibrium and have a certain tendency towards structural relaxation, more than moderately cooled configurations.

4. Conclusions In this paper we presented a three-dimensional model of a polymer melt undergoing the glass transition. To our knowledge it is the first dynamical Monte Carlo simulation with long chains interacting by a van der Waals potential, reaching temperatures far below Tg. After an athermal generation of initial configurations consisting of 40 000 monomers we carried out a continuous cooling procedure with subsequent isothermal annealing. Below Tg the combination of inter- and intramolecular interactions results in compact samples keeping their form even when enlarging the simulation box. At the end of the annealing process we determined several static and dynamical quantities, all of them indicating a pronounced glass transition. We found that, within our simulations, the glass transition is initiated by the energetic component, mainly the van der Waals interaction. This process is in contrast to the geometrically induced glass transition, investigated extensively by Binder and coworkers [45-52]. Variation of chain-length from 49 up to 199 reveals no significant effects. However, this asymptotic independence of Tg from the chain-length is in agreement with the results of the short chain simulation of Lobe et al. [50] which indicates that Tg becomes nearly constant for N > 12. Using Tg(N)

= Tg(oo) -

A/N

with a = 0.23 +_ 0.01

(8) from Ref. [50] we expect a Tg-variation of less than 13e = 0.25 _ 0.02. This variation cannot be resolved by the simulation presented here. A comparison with the molecular weight dependence of polystyrene determined by differential scanning analysis (DSC) and dilatometry [66,67] confirms the above statement. By a rough estimation, interpreting each segment of our model as a group of about seven chemical bonds [68] along the backbone, we obtain a M w range of about

M. Wittkop et aL / Journal of Non-Crystalline Solids 201 (1996) 199-210

36000 to 140000. In this range the changes in the glass temperature are small [66,67]. On the other hand, we were able to reveal the dominant mechanism of vitrification by varying the cooling rate: the faster the samples are cooled, the more pronounced is the tendency that monomers are arrested locally, keeping the system far from a global energy minimum. By sufficient isothermal annealing near Tg, this marked effect of the cooling rate can be removed completely, but below Tg this is not possible.

Acknowledgements Parts of the simulations were performed at the Leibniz Rechenzentrum der Bayerischen Akademie der Wissenschaften Miinchen, H/Schstleistungsrechenzentrum f'tir Wissenschaft und Forschung KFA Jiilich and on the facilities of the Rechenzentrum der Universit~it Regensburg. The authors are grateful for a generous grant of computing time. They also want to thank F. Gotsis for his technical support of this work. Moreover, M.W. (Go 287/181) and T.H. (Go 287/'20-1) would like to thank the Deutsche Forschungsgemeinschaft (DFG) for financial support.

References [1] R. Zallen, The Physics of Amorphous Solids (Wiley, New York, 1983). [2] K. Kawasaki, Ann. Phys. 61 (1970) 1. [3] R.G. Palmer, Adv. Phys. 31 (1982) 669. [4] J. J~ickle, Rep. Prog. Phys. 49 (1986) 171. [5] M.H. Cohen and D. Turnbull, J. Chem. Phys. 31 (1959) 1164; D. Turnbull and M.H. Cohen, J. Chem. Phys. 34 (1961) 120; J. Chem. Phys. 52 (1970) 3038. [6] M.H. Cohen and G.S. Grest, Phys. Rev. B20 (1979) 1077. [7] J.H. Gibbs, J. Chem. Phys. 25 (1956) 185. [8] J.H. Gibbs and E.A. DiMarzio, J. Chem. Phys. 28 (1958) 373; E.A. DiMarzio and J.H. Gibbs, J. Chem. Phys. 28 (1958) 807. [9] J.H. Wendorff, Polymer 23 (1982) 543. [10] D.J. Plazek. J. Non-Cryst. Solids 131-133 (1991) 836. [11] D.J. Plazek, E. Schlosser, A. SchSnhals and K.L. Ngai, J. Chem. Phys. 98 (1993) 6488. [12] R.D. de la Batie, J.-L. Viovy and L. Monnerie, J. Chem. Phys. 81 (1984) 567.

209

[13] B. Frick, D. Richter and C.L. Ritter, Europhys. Lett. 9 (1989) 557. [14] M. B~e, Quasielastic Neutron Scattering (Hilger, London, 1988). [15] F.T. Gentile and U.W. Suter, in: Materials Science and Technology, A Comprehensive Treatment, ed. R.W. Cahn, P. Haasen and E.J. Kramer, Vol. 12, Structure and Properties of Polymers, Volume ed. E.L. Thomas (VCH, Weinheim, 1993) p. 33. [16] J.G. Curro, J. Chem. Phys. 61 (1974) 1203. [17] F.T. Wall, J.C. Chin and F. Mandel, J. Chem. Phys. 66 (1977) 3243; F.T. Wall and W.A. Seitz, J. Chem. Phys. 67 (1977) 3722. [18] H.L. Frisch, M. Bishop, D. Cepedey and M.H. Kalos, J. Macromol. Sci. Phys. BI8 (1980) 453. [19] M. Bishop, D. Ceperley, H.L. Frisch and M.H. Kalos, J. Chem. Phys. 72 (1980) 3228. [20] A. BaumgSstner and K. Binder, J. Chem. Phys. 75 (1981). [21] M. Bishop, D. Ceperley, H.L. Frisch and M.H. Kalos, J. Chem. Phys. 76 (1982) 1557. [22] J.M. Deutsch, Phys. Rev. Lett. 49 (1982) 926. [23] M. Bishop, M.H. Kalos and H.L. Frisch, J. Chem. Phys. 79 (1983) 3500. [24] K. Kremer, Macromolecules 16 (1983) 1632. [25] A. BaumgSa-tner, Ann. Rev. Phys. Chem. 35 (1984) 419. [26] C.C. Crabb and J. Kovac, Macromolecules 18 (1985) 1430. [27] C.C. Crabb, D.F. Hoffmann, M. Dial and J. Kovac, Macromolecules 21 (1988) 2230. [28] A. Kolinski, J. Skolnick and R. Yaris, J. Chem. Phys. 84 (1986) 1922; 86 (1987) 1567, 7164, 7174; J. Skolnick, R. Yaris and A. Kolinski, J. Chem. Phys. 88 (1988) 1407. [29] W. Paul, K. Binder, D.W. Heermann and K. Kremer, J. Phys. II 1 (1991) 37; J. Chem. Phys. 95 (1991) 7726. [30] J. Wittmer, W. Paul and K. Binder, Macromolecules 25 (1992) 7211. [31] J.J. de Pablo, M. Laso and U.W. Suter, J. Chem. Phys. 96 (1992) 2395. [32] K. Kremer, G.S. Grest and I. Carmesin, Phys. Rev. Lett. 61 (1988) 566. [33] K. Kremer and G.S. Grest, J. Chem. Phys. 92 (1990) 5057. [34] H. Takeuchi and R.-J. Roe, J. Chem. Phys. 94 (1991) 7446. [35] D.N. Theodorou and U.W. Suter, Macromolecules 18 (1985) 1467; 19 (1986) 379. [36] A.S. Argon, P.H. Mott and U.W. Suter, Phys. Status Solidi (b)172 (1992) 193; P.H. Mott, A.S. Argon and U.W. Suter, Philos. Mag. 67 (1993) 931; 68 (1993) 537. [37] R.H. Boyd and P.V.K. Pant, Macromolecules 24 (1991) 4073, 4078. [38] C.F. Fan and S.L. Hsu, Macromolecules 24 (1991) 6244; 25 (1992) 266. [39] C.F. Fan, T. (~a~in, Z.M. Chen and K.A. Smith, Macromolecules 27 (1994) 2383. [40] M. Hutnik, F.T. Gentile, P.J. Ludovice, U.W. Suter and A.S. Argon, Macromolecules 24 (1991) 5962; M. Hutnik, A.S. Argon and U.W. Suter, Macromolecules 24 (1991) 5956, 5970; 26 (1993) 1097. [41] D. Rigby and R.-J. Roe, J. Chem. Phys. 87 (1987) 7285; J.

210

M. Wittkop et al./Journal of Non-Crystalline Solids 201 (1996) 199-210

Chem. Phys. 89 (1988) 5280; Macromolecules 22 (1989) 2259. [42] H. Takeuchi and R.-J. Roe, J. Chem. Phys. 94 (1991) 7458. [43] R.-J. Roe, J. Chem. Phys. 100 (1994) 1610. [44] S. Yashonath, R.J. Rao and C.N.R. Rao, Phys. Rev. B31 (1985) 3196. [45] H.-P. Wittmann, K. Kremer and K. Binder, J. Chem. Phys. 96 (1992) 6291. [46] J. Baschnagel, K. Binder and H.-P. Wittmann, J. Phys.: Condens. Matter 5 (1993) 1597. [47] P. Ray, J. Baschnagel and K. Binder, J. Phys.: Condens. Matter 5 (1993) 5731. [48] J. Baschnagel and K. Binder, Physica A204 (1994) 47. [49] J. Baschnagel, Phys. Rev. B49 (1994) 135. [50] B. Lobe, J. Baschnagel and K. Binder, Macromolecules 27 (1994) 3658. [51] J. Baschnagel and B. Lobe, Macromol. Symp. 81 (1994) 63. [52] J. Baschnagel and R. Dickman, J. Chem. Phys. 101 (1994) 3326. [53] W. G&ze J. Phys.: Condens. Matter 2 (1990) 8485; W. G&ze and L. Sj~Sgren, Rep. Prog. Phys. 55 (1992) 241; M. Fuchs, W. G&ze, S. Hildebrand and A. Latz, J. Phys.: Condens. Matter 4 (1992) 7709. [54] A. Milchev, W. Paul and K. Binder, J. Chem. Phys. 99 (1993) 4786.

[55] Carmesin and K. Kremer, Macromolecules 21 (1988) 2819. [56] A.L. Rodriguez, H.-P. Wittmann and K. Binder, Macromolecules 23 (1990) 4327. [57] H.P. Deutsch and K. Binder, J. Chem. Phys. 94 (1991) 2294. [58] J. Baschnagel, K. Binder, W. Paul, M. Laso, U.W. Suter, I. Batoulis, W. Jilge and T. BUrger, J. Chem. Phys. 95 (1991) 6014. [59] J. Baschnagel, K. Qin, W. Paul and K. Binder, Macromolecules 25 (1992) 3117. [60] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys. 21 (1953) 1087. [61] P.E. Rouse, J. Chem. Phys. 21 (1953) 1272. [62] M. Doi and S.F. Edwards, The Theory of Polymer Dynamics (Clarendon, Oxford, 1988). [63] S.F. Edwards, Proc. Phys. Soc. 92 (1967) 9. [64] P.G. de Gennes, J. Chem. Phys. 55 (1971) 572; Scaling Concepts in Polymer Physics (Cornell University, Ithaca, NY, 1979). [65] G. Williams and D.C. Watts, Trans. Faraday Soc. 66 (1970) 80. [66] A. Rudin, D. Burgin, Polymer 16 (1975) 291. [67] M.J. Richardson and N.G. Savill, Polymer 18 (1977) 3. [68] M. Wittkop, S. Kreitmeier and D. G~iritz, J. Chem. Phys. 104 (1996) 3373.