A constitutive model for the nonlinear viscoelastic viscoplastic behavior of glassy polymers.

INTRODUCTION

Thermoplastic materials are being used with increasing frequency in critical load-bearing structural applications in both the homogeneous state and the matrix phase of "tough" advanced composites. It is therefore of great importance to be able to understand and predict the mechanical response of these materials when they are subjected to various loading histories. The nonlinear, viscoelastic, viscoplastic response of glassy polymers is not only rate and temperature dependent, but is also sensitive to the thermomechanical prehistory of the material. The effect of loading history is a consequence of the evolution of the material microstructural state and has been monitored experimentally by changes in internal energy (e.g., 1,2), local free volume (3), and birefringence (e.g., 4). Thermal prehistory dependence, on the other hand, results from the continuous structural relaxation (also termed physical aging) of the metastable material state below the glass transition temperature (5) and can be assessed through measurements of enthalpy content (e.g., 6) and local free volume (3, 7).

A number of treatments exist for the constitutive response of glassy materials (e.g. 4, 8-12). Knauss and Emri (8) used an integral representation of non-linear viscoelasticity with a state dependent variable related to free volume content. The free volume in turn was assumed to be a function of temperature, time, and pressure histories. In the work of Nagi et al. (9), a "coupling model" of the viscoelastic behavior is proposed that also follows an integral representation where the physics of the nonlinear relaxation behavior is captured in a time and stress dependent kernel, and the evolution in structure is captured by the use of a fictive temperature. The concept of a fictive temperature as a state variable that captures the structural (volume)-related effects on mechanical behavior has been widely used to model the state-dependence of polymer deformation. In fact, the Knauss and Emri approach (8) utilizing a free-volume related parameter could be considered equivalent to the fictive temperature approach. The nonlinear behavior is typically modeled using an integral approach and either multiple relaxation times and/or stress dependent terms. However, these models do not consider viscoplastic deformation and strain softening, although the question of yield has been addressed by other workers such as Argon (10), Robertson (11), Bowden and Raha (12), and G'Sell and Jonas (13). In some of these models, the plastic strain rate has been based on an explicit mechanism for the yielding phenomenon; also, some of the authors mention or discuss the distributed nature of the plastic flow mechanism but do not formalize it within the context of a probability density function. Argon (10) uses a disclination(1) dipole to model polymer chain segment rotation. Robertson (11) considers the transfer of polymer chains from the low energy trans state to the high energy cis state occurring at a critical rate to signify yield; and Bowden and Raha (12) have based their model on the nucleation and expansion of dislocation loops. As such, there exists a need to represent the nonlinear viscoelastic and post-yield viscoplastic response of glassy polymers and the effect of thermo-mechanical prehistory on the mechanical response using a single, unified model representation. This investigation is aimed at ultimately achieving just such a model, using the experimental data of the material state (23), the results from molecular dynamics and Monte Carlo simulations on amorphous materials (e.g., 14, 15), and the mechanical behavior presented in this paper for a range of thermal prehistories as both motivation and basis.

EXPERIMENTS

Uniaxial compression, constant true strain rate experiments, and compressive creep tests under nearly constant true stress conditions were conducted. The initial states used for these tests were "annealed" (oven-cooled to room temperature from [T.sub.g] + 20 [degrees] C) and "quenched" (brine-quenched to room temperature from [T.sub.g] + 20 [degrees] C). These tests were performed on an Instron 1340 servohydraulic machine using a Macintosh for data acquisition and control via a Keithley 550 D/A and A/D interface. In order to eliminate machine compliance errors, an Instron 2620-824 extensometer was used to directly monitor the change in specimen height. All the constant strain rate stress-strain tests were done in strain control by sending a string of command strain signals to the PD controller. For the creep tests, which were done in load control mode, the desired command signals cannot be known apriori because of the continuous change in specimen dimensions. In this case, the specimen height was monitored at 1 Hz frequency (using the extensometer); then assuming constant volume, the new load required to maintain constant true stress was computed and sent to the controller. The fact that Poisson's ratio during elastic loading is about 0.33 instead of 0.5 means that the correct (absolute) value of true creep stress is consistently higher by 0.5-1.0 MPa compared to the value calculated during the test. We have included this second order correction in reporting our data. In this way, the true stress was maintained at a nearly uniform value throughout the test, even for the larger stresses. Finally, for the elevated temperature tests, a radiant furnace was used in conjunction with a West 2073 controller to maintain the specimen temperature to within [+ or -] 2 K.

Figures 1 and 2 show the results of isothermal, constant strain rate, uniaxial compression tests on annealed and quenched PMMA at two different temperatures. The essential features of the typical mechanical response observed in these experiments is the initial elastic behavior, followed by a viscoelastic response, which becomes increasingly nonlinear with increasing stress and strain. Eventually, the stress is large enough for macroscopic yield to occur, followed by isothermal, material strain softening to a steady-state flow stress. Another important feature of the experimental data shown in Figs. 1 and 2 is that global viscoplastic flow (we define the "yield" stress to be the peak stress) begins at a lower stress level in the quenched material than in the annealed polymer by virtue of its having a more "disordered" state. The effect of thermal prehistory on the mechanical response has been entirely erased after strain softening and steady-state flow occurs at the same stress level in both annealed and quenched polymer [see Hasan and Boyce (2) and Hasan et al. (3) for additional details on aging and strain softening effects on behavior and structure). The rate and temperature dependence of the mechanical response is consistent with the thermally activated nature of inelastic deformation in glassy polymers.

Figure 3 depicts the nonlinear unloading response of the material following different amounts of straining. Note that a fresh specimen is used to obtain each (load-unload) curve. Prior to "yield," full recovery (highly nonlinear) is obtained within a few minutes. Beyond the "yield" point, a substantial portion of the inelastic strain is found to be "permanent," but the unloading is still highly nonlinear and reaches a more or less stable unloading profile. The nonlinearity in unloading [also observed by Oleynik (1)] is indicative of an energy storage phenomenon that drives inelastic strain recovery upon unloading as the applied stress decreases. The presence of stored internal energy beyond that associated with linear elastic deformation has been measured by Oleynik (1) and Adams and Farris (16) during deformation using deformation calorimetry.

Figures 4 and 5 show the creep response under uniaxial compression, over a range of true stress levels at 296 and 323 K. The initial elastic response is followed by inelastic flow at a progressively decreasing rate for low values of creep stress (usually termed primary and secondary creep in the metals literature). For large stresses, however, following a certain amount of straining, the inelastic strain rate increases by orders of magnitude (a condition termed tertiary creep) and corresponds to large-scale flow in the material. Figure 6 shows a typical variation in the applied load and the true stress during the creep experiment. The rapid increase in load at larger times compensates for the changing cross-sectional area of the specimen.

CONSTITUTIVE MODEL

Model Background

Local rearrangements in metallic and glassy materials are thermally activated in nature, an observation dating back to such pioneers as Arrhenius and Boltzmann. Moreover, thermal activation is a stochastic process and is therefore formulated in terms of the activation free energy (which is based on the physical nature of the rearrangement), the thermal energy and the attempt frequency. A thorough discussion of the role and modeling of thermal activation in inelastic flow (in polycrystalline materials) has been given by Kocks, Argon, and Ashby (17). Following these considerations, the transformation rate, [Omega], typically takes the following form:

[Omega] = [[Omega].sub.0] exp {-[Delta][F.sub.f]/kT} (1)

where [[Omega].sub.o], the attempt frequency, is on the order of the Debye frequency; [Delta][F.sub.f] is the activation free energy for the transformations; k is the Boltzmann constant; and T is the absolute temperature. The activation energy is modified by any externally applied and internal stresses. Several expressions have been proposed in the literature regarding the explicit nature of this stress dependence. As discussed in Kocks, Argon, and Ashby (17), albeit with respect to polycrystalline plasticity, most of these take the form

[Delta][F.sub.f] = [Delta] [F.sub.o] [{1 - [([Tau]/[[Tau].sub.o]).sup.p]}.sup.q] (2)

where [Delta][F.sub.o] is the activation energy in the absence of any stress, [Tau] is the local stress (often taken to equal the applied stress), [[Tau].sub.o] is the limiting strength, and p and q are material dependent parameters, which are, in theory, determined by the mechanism of inelastic deformation. The effective lowering of the activation barrier in the above equation has traditionally been interpreted within the context of mechanical work done by the stress field over the transformation volume. In the model by Argon (10) for glassy polymers, p = 5/6 and q = 1; in the model proposed by Mangion et al. (18) for glasses, p = 1 and q = 3/2. Setting both p and q equal to unity results in the Eyring equation, albeit in a form slightly different from the one quoted commonly.(2) Our experimental data for glassy polymers indicates that for the range of temperatures and deformation rates mentioned above, the choice of p and q has only a marginal effect on the predictive capabilities of the thermal activation model. More specifically, we have evaluated the model to be presented below using 1) p = q = 1, 2) p = 5/6, q = 1, 3) p = 2/3, q = 1 and found little effect on the predictions of the rate and temperature dependence. For simplicity, therefore, we have based our current modeling approach on the Eyring form of the thermal activation model.

The above treatment was developed within the context of viscoplastic flow; local rearrangements in glassy polymers, on the other hand, are thermore-versible. As such, one must include a mechanism for reverse transformations (strain recovery) to occur within the material. Apart from their thermally activated nature, such reverse transformations will be opposed by the applied stress. On the other hand, they will be facilitated by the presence of internal or "transformation" strain energy stored in the region around the locally transformed sites, a feature that has often been interpreted within the context of a back stress. Inclusion of reverse transformations, or rearrangements that influence overall flow at low stresses, requires the addition of the reverse term in Eq 1, giving

[Omega] = [[Omega].sub.o] [absolute value of exp {-[Delta][F.sub.f]/kT} - exp{-[Delta][F.sub.b]/kT}] (3)

where [Delta][F.sub.b] is the activation energy for the reverse transformation.

To underscore the importance of the constitutive model presented here, we begin by briefly reviewing earlier work done by the authors in this area (3). There, the effects of physical aging and strain soften-ing on the mechanical response were modeled through a modification of the Argon thermal-activation model (10):

[Mathematical Expression Omitted]

[Delta]F = [Delta][F.sub.o] (1 - [([Tau]/[[Tau].sub.o]).sup.5/6) (5)

[Mathematical Expression Omitted]

Here [Mathematical Expression Omitted] is the viscoplastic strain rate; [Mathematical Expression Omitted] is the pre-exponential factor given by the product of the number density of shear transformation sites (D), the volume-averaged shear strain increment per transformation (2[[Gamma].sup.T][Omega]/V) and the attempt frequency ([[Omega].sub.o]); [Tau] is the applied shear stress and the other variables have been defined earlier. Using this model, the effects of aging and strain softening were captured by using D, the number density of transformable sites as a scalar state variable (3). D was phenomenologically modeled to increase with plastic straining until a "preferred" internal structure is achieved:

[Mathematical Expression Omitted]

where [Mathematical Expression Omitted] and [Mathematical Expression Omitted] indicate evolution due to aging and straining, respectively. The rate and temperature dependence of the yield stress and the strain softening after yield were well predicted using this approach (3). However, the use of a single, scalar-valued internal variable resulted in an abrupt linear elastic to viscoplastic transition at a strain of [similar] 4%, whereas yield is observed between 8 and 10% strain and is preceded by a highly nonlinear transition between linear elastic behavior and the yield point. The experimentally observed nonlinear viscoelastic to viscoplastic flow conditions suggest that a more realistic material representation must be employed. Another shortcoming of the model was that it did not account for the thermoreversible nature of the shear transformation and was, therefore, not capable of predicting nonmonotonic loading or strain recovery in the unloaded state. These are some of the issues we aim to address in the current work.

Our differential scanning calorimetric (DSC) (2) and positron annihilation lifetime spectroscopic (PALS) (3) measurements on glassy polymers indicate that the microstructural state is highly distributed in nature. This state may be characterized in a number of (mechanistically and operationally) equivalent ways, such as a distribution of local free volume, a distribution of internal stress due to nonuniform polymer chain packing, or a distribution of activation energies for the thermalized, local inelastic shear transformations that occur in the glassy state. The equivalence of these different treatments is a consequence of the fact that one should, at least in principle, be able to relate one method of characterizing material state to the others. In fact, limited success has been obtained in molecular dynamics simulations of simplified glassy systems, where a direct correlation has been observed between local free volume, transformation frequency and local stress (14). In this work, we focus on using a distribution of activation energies to characterize the energy barrier to localized shear transformations. These transformations correspond to the rearrangements that occur in the material during inelastic deformation and structural relaxation; however, in our modeling approach, we will not be considering any specific mechanism of these shear events. The typical bond energy changes associated with polymer chain rotation, in conjunction with the specific enthalpy changes associated with inelastic deformation, suggest that the shear transformations are localized in nature. Atomistic simulations on glassy polymers by Mott (19) suggest that the transforming volume is on the order of a few ten thousand cubic angstroms. In Mott's simulations, the number of degrees of freedom in the simulation cell were restricted for reasons of computational expediency; therefore, we expect the actual transformation volume to be perhaps somewhat smaller than that obtained using his simulations.

Proposed Model

Prior to outlining our proposed model, we consider the pictorial representation of the evolution in material "state" during nonlinear viscoelastic and vis-coplastic deformation [ILLUSTRATION FOR FIGURE 7 OMITTED]. Here, a typical global stress-strain curve up to strains of about 35% and the corresponding distribution of material state has been shown at different stages of deformation. We take the distribution in material state to be directly related to the distribution in, for example, local free volume, local energy barrier to deformation, or associated shear transformation sites where we consider high local free volume to correspond to low activation barrier and therefore greater probability of a local shear transformation event. In Figure 7, the size of the circles is proportional to the probability of occurrence of a local, strain-producing transformation event and the dots represent the locally stored transformation strain energy. Figure 7-i shows the initial configuration with no net strain and little or no transformation strain energy. For this state, there are a number of sites with meaningful transformation probability within the time frame of the experiment. During the initial stages of deformation (at stresses that are low compared with the flow stress), only the regions of high local free volume (low activation energy) can transform at meaningful rates [ILLUSTRATION FOR FIGURE 7-ii OMITTED] accompanied by a departure from linearity in the stress-strain response. The corresponding transformation strain energy is stored elastically in the relatively rigid surrounding material, where the local free volume is much lower. The surrounding material, in turn, exerts a back stress on the transformed regions with two consequences. First, the back stress prevents multiple shear transformations in the same material volume by increasing the effective activation energy required for subsequent transformations. Second, upon unloading, it enables rapid recovery of the strain associated with these transformations, thus resulting in the nonlinear unloading behavior. As the stress is increased, regions of progressively higher activation energies (and lower local free volume) become accessible for shear transformations [ILLUSTRATION FOR FIGURE 7-iii OMITTED] and the mechanical behavior becomes increasingly nonlinear. Moreover, the "matrix" material available for the elastic storage of the transformation strain energy is depleted to the point that additional energy can no longer be stored in the material by this mechanism. When this happens [ILLUSTRATION FOR FIGURE 7-iv OMITTED], further energy is now stored in the material by the creation of new defects (high local free volume sites). The creation of new "soft" sites results in a decrease in the macroscopc flow stress (material strain softening) and is reflected in the shift of the transformation probability spectrum to lower activation energies. During strain softening, the material state evolves to a steady-state flow condition for which the regions of liquid-like mobility are sufficiently numerous to enable indefinite, "plastic" flow [ILLUSTRATION FOR FIGURE 7-v OMITTED].(3) This scenario is supported by the molecular dynamics simulations on atomic level glasses of Deng et al. (14), the more recent Monte Carlo simulations of Bulatov and Argon (15), and our own experimental results.

We now translate the pictorial representation of nonlinear deformation into a mathematical model. We begin by proposing generic local thermoreversible rearrangements (time, temperature, strain and stress dependent) with [Psi] denoting the fraction of material available for such rearrangements. Then, assuming thermal activation and first order kinetics, we can denote the net rate of change of available processes as

[Mathematical Expression Omitted]

where [[Omega].sub.o] is the attempt frequency and [Delta][F.sub.f] and [Delta][F.sub.b] are the activation energies for the forward and reverse transformations, respectively. However, the above equation is valid only if all possible molecular rearrangements have the same activation energy. In order to represent the distributed nature of local rearrangements, [Psi] must be distributed continuously over activation energy space. Hence, the above rate equation must be modified to:

[Mathematical Expression Omitted]

where [Psi]([Delta][F.sub.f], [Delta][F.sub.b]) is now a probability density function in the mathematical sense.

[Mathematical Expression Omitted]

Setting [Psi] = 1 means that the full material volume is available for shear transformations, at all times. This does not result in any loss of generality because the probability of transformations of some of the material volume (corresponding to high thermal activation energies) will always be extremely low.

If we associate an inelastic strain increment of +[Delta][Gamma] with the forward transformation and that of -[Delta][Gamma] with the reverse transformation, the macroscopic inelastic strain rate is given by

[Mathematical Expression Omitted]

where

[Mathematical Expression Omitted]

Here we make the assumption that the strain increment per transformation is independent of the activation energy of the transformation. The total number of available processes, D, is computed from the reciprocal of the activation volume [Omega]; [Delta][Gamma] is the transformation strain increment measured over the material volume; [[Gamma].sup.T) is the local transformation strain on the barrier configuration and V is the total volume. Additionally, if we assume an Eyring form for the stress dependence of the activation free energy,(4) the rate equation modifies to

[Mathematical Expression Omitted]

where [Tau] is the equivalent shear stress and [Mathematical Expression Omitted] is the shear activation volume.(5) which we have observed to scale with the shear modulus, [Mu], for a fairly wide range of temperatures below [T.sub.g], i.e.

[Mathematical Expression Omitted]

This temperature scaling is also in accord with the scaling of [[Tau].sub.o] in the Argon model for yield (Eq 5), where [[Tau].sub.o] is linearly proportional to [Mu]. We also note that this shear activation volume is more accurately interpreted as the product of the transformation strain and the volume of the transforming site [see Mott (19)]. The second term in Eq 13 characterizes the probability of the reverse transformation, which is opposed by the externally applied stress but facilitated by the local transformation strain energy, S, stored around the transformed sites. Note that S is taken only to affect the probability of reverse transformations and not that of forward transformations for two reasons: 1) the presence of S hinders forward flow, and therefore sites without S are the primary contributors to forward flow, and thus S is not included in the first term; 2) the presence of S aids reverse flow, and therefore sites with S are the primary contributors to reverse flow, and thus S is included in this term. Making an explicit choice regarding the stress dependence of the activation energy and the role of the transformation strain energy in modifying the rate equation enables [Psi] to be expressed as a function of a single random variable, [Delta]F, instead of two random variables, [Delta][F.sub.f] and [Delta][F.sub.b].

To evaluate the performance of the rate equation presented above, a choice must be made regarding the explicit form of the probability density function.(6) For the moment, we choose a "pseudo-Gaussian" distribution of the type shown in Fig. 8. Unlike the Gaussian, this distribution has a finite "bandwidth" and thus does not give a non-zero probability for the occurrence of processes with negative activation energies. The probability density function is constructed from the monotonic parts of [e.sup.x]sin(x) and [e.sup.-x]sin(x) type terms in the range 0 to [Pi] and is completely specified by the two parameters a and [[Alpha].sup.-1], which denote the "position" and "standard deviation," respectively. Thus, we have

[Psi] ([Delta]F) = 1/2[A.sub.1] exp {[Alpha]([Delta]F-a)} sin [[Alpha]([Delta]F - a)];

a [less than or equal to] [Delta]F [less than or equal to] a + 3[Pi]/4[Alpha]

= 1/2[A.sub.2] exp{-[Alpha]([Delta]F - d)} sin [[Alpha]([Delta]F - d)];

a + 3[Pi]/4[Alpha] [less than or equal to] [Delta]F [less than or equal to] a + 3[Pi]/2[Alpha]

= 0 otherwise

(15)

where [A.sub.1] and [A.sub.2] normalize the area under the probability density function;

[A.sub.1] = 1/2[Alpha] [1 + [square root of 2] [e.sup.3[Pi]/4]] = [e.sup.[Pi]][A.sub.2] (16)

and

a[prime] = a + [Pi]/2[Alpha] (17)

Substituting the form of [Psi]([Delta]F) from Eq 15 into the expression for inelastic strain rate (Eq 13) and integrating over the rage of activation energies gives the expression for overalll inelastic strain rate as a function of temperature, T, equivalent shear stress, [Tau], and the internal variables, a, [[Alpha].sup.-1], S

[Mathematical Expression Omitted]

[Mathematical Expression Omitted]

where d is defined in Eq 17 and [Eta] is the reciprocal of the thermal energy, kT. The overall inelastic rate is then given by the sum of the two rate equations

[Mathematical Expression Omitted]

In our model, we have sought to capture the evolution of microstructural disorder using the three internal state variables, two of which (a and [[Alpha].sup.-1]) characterize the mean and standard deviation of the distribution of disorder, respectively, and a third, S, which is related to the mean local transformation strain energy in the material. Referring to our pictorial representation of Fig. 7, we note that increasing the effective activation energy using the evolution of [[Alpha].sup.-1] models the nonlinear viscoelastic response to be a consequence of regions of relatively high local free volume being rapidly depleted with increasing inelastic strain. On the other hand, by decreasing the effective activation energy, the evolution of a represents the creation of additional regions of high local free volume due to free volume redistribution in the vicinity of previously transformed sites. This can occur at an appreciable rate only when the fraction of transformed material has reached a sufficiently large value (3). The combined evolution of a and [[Alpha].sup.-1] corresponds to a broadening of the activation energy distribution. We believe that the distributed nature of microstructural disorder is a key feature of the material state. As such, one must explicitly address the evolution of this distribution with inelastic deformation, something which to our knowledge has not been done so far within the framework of a general-purpose constitutive model. In fact, simply including a distribution but not accounting for its evolution results in identical behavior as that found using a single-valued effective activation energy.

In considering the evolution of the mean stored local strain energy, S, we note that its equilibrium value under zero applied stress must be zero to prevent spontaneous shear deformation in the unloaded state; however, when the equivalent shear inelastic strain is non-zero, it will govern strain recovery. During the initial stages of deformation (of previously undeformed material), S will increase rapidly as energy is stored locally by inelastic shear transformations of spatially well separated sites. With increasing inelastic strain, the proximity of the sites will result in a partial dissipation of the stored energy.(7) Now the locally stored transformation strain energy will actually decrease somewhat, before achieving a constant value during steady-state flow conditions.

Based on the above understanding of local material state, we now consider the evolution laws for the internal variables; if we use s to denote a generic internal variable (e.g. a, [Alpha], S), then the overall rate of evolution of s is given by

[Mathematical Expression Omitted]

where the first term represents the change resulting from inelastic deformation and the second term, that due to physical aging. Assuming first order kinetics, we can denote the contribution of deformation and aging to the evolution of the internal variables as

[Mathematical Expression Omitted]

and

[Delta]s/[Delta]t = -[s - [s.sub.eq](T)] [[Omega].sub.aging](T, s) (23)

where [Omega] is the effective frequency of inelastic shear transformations and [[Omega].sub.aging](T, s) is the rearrangement rate during physical aging. For the mechanical tests presented here, the contribution of physical aging to the overall evolution rate, over the time scale of the experiments, is negligible. Hence, we simplify the evolution of the three internal variables as

[Mathematical Expression Omitted]

[Mathematical Expression Omitted]

[Mathematical Expression Omitted]

where [Omega], the effective frequency for inelastic shear transformations, follows directly from Eqs 3 and 13

[Mathematical Expression Omitted]

where

[Phi] ([Delta]F, T) = [Psi]([Delta]F) exp{-[Delta]F/kT} (28)

denotes the probability distribution of available shear transformation events and f([[Gamma].sup.p]) and [S.sub.eq] are given by

f([[Gamma].sup.p]) = exp(- [Zeta] exp(- [Zeta][[Gamma].sup.p])) (29)

[Beta] = [[Beta].sub.1][1 + [[Beta].sub.2] exp(-[[Beta].sub.3] [[Gamma].sup.p])] (30)

Here, the double exponential function in Eq 29 reflects the fact that creation of high local free volume sites occurs at a negligible rate until the transformed sites are sufficiently numerous (or the inelastic strain is sufficiently large). Thus, f([[Gamma].sup.p]) is almost zero for small strains and increases rapidly to unity in the vicinity of the yield phenomenon. Next, the evolution of the locally stored transformation energy is the result of two competing effects: an increase in S due to (partial) storage of the inelastic work of deformation and a decrease in S due to recovery phenomena. While the first effect dominates during the initial stages of deformation, once the value of the locally stored transformation strain energy is sufficiently large, secondary relaxation events start occurring at an appreciable rate and the value of S falls some, what; during steady-state flow conditions, the two effects virtually cancel each other to give no net change in the value of S. Finally, the variation of the parameter fi in Eq 30 is a reflection of the experimental observation that the rate of storage of inelastic work falls off with increasing inelastic strain.

The overall strain rate is given as the sum of the elastic and inelastic contributions

[Mathematical Expression Omitted]

where [Mu] is the shear modulus. Although the above formalism is in terms of shear stress and strain, we can easily predict uniaxial compression data by equating the work done under uniaxial compression to the work done under simple shear loading. Using a Mises equivalent representation, the shear stress is scaled by [square root of 3]/2 and the shear strain by 2/[square root of 3]. Also, pressure effects can be included in the model in the manner described earlier.

This completes the formal outline of the model (where the essential equations are Eqs 15 through 20 and 24 through 31). We now evaluate the model for a practical case.

RESULTS AND DISCUSSION

In order to obtain the needed material properties for the model and also to assess its predictive capabilities, both stress-strain and creep tests were conducted on PMMA, poly(methyl methacrylate): These have been presented earlier in this paper. We now use the uniaxial compression data to determine the model parameters and then evaluate the model predictions for both the nonmonotonic constant strain rate tests and the creep tests.

Material Property Determination

We shall use the data in Figs. 1 and 2 to determine the model parameters, which are [a.sub.eq], [Zeta], [[Alpha].sup.-1](0), [Mathematical Expression Omitted], [[Beta].sub.i], [Mathematical Expression Omitted] and [Mathematical Expression Omitted]. We also note that the fundamental attempt frequency, [[Omega].sub.o], is on the order of the Debye frequency ([10.sup.11] to [10.sup.13] [s.sup.-1]) and we shall use a value of [10.sup.11] [s.sup.-1] without any adjustment. The shear modulus [Mu] has been determined from dynamic mechanical testing to be about 1.25 GPa at room temperature. The pre-exponential factor [Mathematical Expression Omitted] is determined from the intercepts of a 1n [Mathematical Expression Omitted] vs. [Sigma]/(kT) plots of the steady-state flow stress data in Fig. 1, assuming a Dirac delta function in Eq 11. This yields a value of 1.51 x [10.sup.10]/s for [Mathematical Expression Omitted] and 1324 cubic angstroms for the shear activation volume [Mathematical Expression Omitted] at 296 K. f([[Gamma].sub.p]) can be made to approximate the unit step function at the typical yield strain of glassy polymers using a value of 30 for [Zeta]. Finally, in determining the parameters specifying the distribution, we make use of the experimental measurements of local free volume distribution (22). These results indicate that it is primarily the mean local free volume, and not the standard deviation of the local free volume, that is sensitive to thermal prehistory. Thus, we take the effect of thermal prehis-tory to only affect the initial value of a. The values of the distribution parameters, determined by a least-squares fit of the experimental data, are summarized in Table 1. Future work will be aimed at quantifying this aging dependence of the internal variables and providing the necessary evolutionary laws explicitly.

[TABULAR DATA FOR TABLE 1 OMITTED]

Stress-Strain-Test Simulations

Using the above values of the material properties, we simulate the experimental data of Figs. 1 and 2. The results of the stress-strain simulations are shownin Figs. 9 and 10. The corresponding evolution of the distribution of available shear transformations and the internal variables for the annealed material (at room temperature) are shown in Figs. 11 and 12.

The model simulations shown in Figs. 9 and 10 indicate that the model is capable of describing the nonlinear transition from viscoelastic deformation toviscoplastic flow. It also captures the temperature, rate, and thermal prehistory dependence of the material response, as well as the material strain softening very well. However, it is somewhat deficient in the precise quantitative prediction of the onset of nonlinear viscoelastic deformation prior to yield. This is a consequence of the simple probability density function we have used in this work and may be rectified by using a distribution with more than two parameters. Figure 11 depicts the evolution of the probability of transformation of processes with a given activation energy, as given by [Phi]([Delta]F). Following the definition given in Eq 28, this will depend not only on temperature and activation energy, but also on the number density of available processes with the specified activation energy. During the initial stages of deformation, the evolution of [Phi]([Delta] F, T = 296 K) corresponds to the rapid depletion of regions of (relatively) high local free volume available for shear transformations, coinciding with the observed and predicted non-linear viscoelastic behavior prior to yield; consequently, only regions of lower local free volume (higher activation free energy) are available for subsequent transformations and are accessed as the stress level increases. However, once the fraction of transformed sites is sufficiently large for substantial interaction and free volume redistribution to occur, macroscopic yield occurs (the strain is now between 5% and 10%). Thereafter, there is an increase in the density of high local free volume sites available for shear transformations (3), and this is accompanied by an increase in the fraction of material available for transformations at low activation energies [ILLUSTRATION FOR FIGURE 11 OMITTED], thus producing the observed and predicted strain softening behavior.

The features of the distribution of processes available for shear transformations are also depicted clearly by the changes in the internal variables a and [[Alpha].sup.-1], representing "position" and "standard deviation," respectively (the evolution of S is considered later, in conjunction with the nonmonotonic loading simulations), with increasing true strain [ILLUSTRATION FOR FIGURE 12 OMITTED]: During the first 1% of strain, the stress is too low for thermally activated evolution of the internal variables to occur over the time scale of this test. Thereafter, a very rapid increase in [[Alpha].sup.-1] is observed. This corresponds to the depletion in the number of low activation energy sites available for transformations, as also shown in Fig. 11. Once the number of transformed sites is sufficiently high (at 6-8% strain), their mutual interaction dissipates the back stress (as shown later) and permits the evolution of regions of high local free volume sites, corresponding to the dramatic reduction in a. This results in macroscopic yield and strain softening. Eventually, the material attains a steady-state flow condition, with virtually uniform values of the internal state variables. Figure 12 also compares the evolution of the internal variables for annealed and quenched initial states. Our assumption that physical aging primarily affects the density of high local free volume sites and not the standard deviation of local free volume size distribution means that the evolution of [[Alpha].sup.-1] is virtually identical for the annealed and quenched polymer. However, the initial value of a is lower for the quenched material than for the annealed material for the reason just mentioned. For both annealed and quenched polymer, the value of the internal variables is the same at the end of strain softening, which is consistent with experimental characterization of the material state (2, 3). Finally, the evolution of the internal variables during the 50 [degrees] C simulations is virtually identical to the room temperature simulation results. This is consistent with our current assumption that the initial and equilibrium values of the internal variables are the same at both temperatures.

Next we model the nonmonotonic mechanical response shown earlier in Fig. 3. The simulations are shown in Fig. 13 and the corresponding evolution of the internal variables is shown in Fig. 14a-c. The results in Fig. 13 indicate that the model captures the nonlinear unloading behavior rather well. The nonlinear unloading is driven by the stored energy, S, which activates reverse flow as the applied stress decreases. The strain recovery, in turn, is responsible for further evolution in the state variables. Figure 14a shows the decrease in a during monotonic inelastic straining, and the partial recovery during unloading, and a more gradual recovery in the unloaded state. Similar changes are observed for the other two internal variables in Figs. 14b-c. The evolution of the internal variables suggests that inelastic strain recovery acts to reverse the local changes induced in the material during monotonically increasing strain. Moreover, the finite value of stored energy S, at the end of unloading, is responsible for additional stress-free strain recovery, as observed in our experiments as well as those of Oleynik (1). These results indicate the ability of the proposed modeling framework to predict nonmonotonic deformation and strain recovery in glassy polymers. As a further test of the model, the parameters determined from the uniaxial compression tests were used to predict experimental results for uniaxial creep tests.

Creep Test Simulations

Using the model properties determined from the stress-strain data in Fig. 1, we predict the deformation response during constant (compressive) stress creep tests on annealed PMMA. The experimental data (solid lines) and the model predictions (dotted lines) are shown in Figs. 15a-b over a range of stress levels at 296 K and 323 K. The corresponding evolution of the internal variables a and [[Alpha].sup.-1] is shown in Figs. 16a-b (the evolution of S is similar to that shown in Fig. 14). It can be seen that the predictive capability of the model is quite reasonable. The model predicts the initial elastic response followed by inelastic flow. If the creep stress is low ([less than] 52 MPa for 323 K, [less than] 74 MPa for 296 K), the rate of inelastic flow eventually stabilizes to a value consistent with the rejuvenation rate of sites that have just been transformed, as shown by the evolution of a and [[Alpha].sup.-1] in Figs. 16a-b. The stabilization strain for these stress levels is not precisely predicted as a result of our not quite precise prediction of the nonlinear transition in the stress-strain behavior discussed earlier. As mentioned there, a distribution more sophisticated than the two-parameter distribution used here would rectify these discrepancies. We have not done that here because the primary aim of the current work is not to achieve a remarkable fit to the experimental data using ever increasing number of parameters, but to emphasize the importance of addressing the distributed nature of microstructural response in developing constitutive models for glassy materials. For high stresses, the model correctly predicts the gradual acceleration of the inelastic flow rate to a substantially higher value, as well as the temperature dependence of this transition. This behavior is quite similar (at the local material scale) to the yield phenomenon observed during uniaxial compression testing. In both cases, a sufficiently large number of shear transformations have occurred by this strain level to result in a transition from an essentially frozen (glassy) state with isolated regions of liquid-like mobility to one in which contiguous regions of liquid-like mobility permeate the microstructure.

The evolution of the internal variables for the creep tests [ILLUSTRATION FOR FIGURE 16A-B OMITTED] is consistent with the trends observed for the case of constant strain rate tests. For short times, the progressive decrease in strain rate with time signifies the depletion of regions of high local free volume (modeled by an increase in [[Alpha].sup.-1]). If the stress is not high enough, the strain virtually levels off as all the sites transformable within the time frame of the experiment are depleted. In this case, the value of a remains virtually unchanged. If, on the other hand, the creep stress is close to the typical yield stress at that temperature, a sufficiently large number of sites are transformed to enable "yield" to occur; this is characterized by a dramatic increase in the strain rate and is captured in the model by the reduction in the value of a. The rapid evolution in the value of S during the initial stages of deformation plays a significant role in the strain recovery during unloading. Note again that the creep response was predicted using the property values determined from the uniaxial compression data, without any adjustment.

CONCLUSIONS AND FUTURE WORK

A simple model of the glassy state based on the distributed nature of microstructural disorder has been developed. The state is characterized by the distribution of activation energies for processes available for inelastic shear transformations and by the local strain energy stored in the materiall following such transformations. By using three internal state variables and the corresponding evolution equations, the temperature, rate, and thermal prehistory dependence of the nonmonotonic deformation response during uniaxial compression and creep tests have been satisfactorily modeled. Equally important is the fact that this has been done using a very simple, physically motivated constitutive framework using two internal state variables to characterize the distributed nature of microstructural disorder and one additional internal variable to represent the stored local energy. The material properties needed in the model are obtained from a small set of stress-strain curves, and the model was found to be predictive of the nonlinear viscoelastic-viscoplastic behavior under nonmonotonic loading conditions over a range of strain rates and temperature as well as predictive of the creep behavior over a wide range in stress level at different temperatures. The model results and predictions are very satisfactory, although not perfect, and indicate the promise of the proposed framework of a thermal activation model considering an evolutionary distribution in activation energies to provide a single unified model to predict nonlinear viscoelastic and viscoplastic deformation of amorphous polymers. Better descriptions of the actual distribution in the energy barrier should provide a more precise correlation between the model and the experimental results. Experiments probing this distributed nature of the material structure [e.g., Jean (21) and Li (22)] indicate the promise of obtaining such descriptions in the near future.

ACKNOWLEDGMENTS

This work was supported by the NSF through a PYI award to MCB with matching funds from the DuPont Company through a DuPont Faculty Award to MCB, and, in part, by the NSF through grant number MSS-9215805.

1 Disclinations are rotational defects as opposed to dislocations, which result from translational mismatch.

2 In the Eyring formulation, the stress-reduced activation energy is typically expressed as [Delta][F.sub.f] = [Delta][F.sub.o] - [Tau][Delta]v. This provides an equivalent representation to that of Eq 2, where p = q = 1, when the activation volume can simply be represented as [Delta]v([equivalent] [Delta][F.sub.o]/[[Tau].sub.o]).

3 One may also interpret this evolution of material state as the "percolation" [see e.g., Stauffer and Aharony (13) of regions of liquid-like mobility through the glassy matrix, to eventually form a continuous network.

4 In this current form, we neglect the pressure dependence of yield. However, it can easily be included either by introducing a pressure activation volume giving [Mathematical Expression Omitted] as the stress modified activation energy or by using a pressure dependent shear activation volume. In the experiments to be presented and discussed, this term represents a second order effect and is not included in the current formulation.

5 As discussed earlier, the modification of the activation energy by an applied stress can equivalently be viewed either within the context of a limiting strength or as the work done by the applied stress.

6 Our analysis indicates that in principle, [Psi]([Delta]F) can be determined experimentally; unfortunately, the technique requries determining the initial slope (initial inelastic strain rate) of constant true stress, isothermal creep tests. This is not meaningful given the fact that the initial slope of creep data varies by orders of magnitude over a small time interval.

7 This dissipated energy is used to create additional regions of easily transformable material (experimentally verified by an evolution (increase) of mean local free volume with deformation (3).

REFERENCES

1. E. Oleynik, High Performance Polymers, p. 79, E. Baer and S. Moet, eds., Hauser, Munich (1990).

2. O. A. Hasan and M. B. Boyce, accepted for Polymer (1993).

3. O. A. Hasan, M. B. Boyce, X. S. Li, and S. Berko, J. Polym. Sci., Polym. Phys. Ed., 31, 185 (1993).

4. E. M. Arruda and M. B. Boyce, to appear in Intn. J. Plas.

5. L. E. Struik, Physical Aging in Amorphous Polymers and Other Materials, Elsevier, Amsterdam (1973).

6. I. M. Hodge, Macromolecules, 16, 898 (1983).

7. X. S. Li, PhD thesis, Brandeis University (1992).

8. W. Knauss and R. Emri, ASME Winter Meeting (1989).

9. R. W. Rendell, K. L. Ngai, G. R. Fong, A. F. Yee, and R. J. Bankert, Polym Eng. Sci., 27, 2 (1987).

10. A. S. Argon, Phil. Mag., 28, 839 (1973).

11. R. E. Robertson, J. Chem. Phys., 44, 3950 (1966).

12. P. B. Bowden and S. Raha, Phil. Mag., 22, 463 (1974).

13. C. G'Sell and J. J. Jonas, J. Mater. Sci., 16, 1956 (1981).

14. D. Deng, A. S. Argon, and S. Yip, Phil. Trans. Roy. Soc., A329, 613 (1989).

15. V. Bulatov and A. S. Argon, "A Stochastic Model of Elasto-Plastic Continuum," accepted for Model. Simul. Matr. Sci. Eng.

16. G. W. Adams and R. J. Farris, J. Polym Sci., 26, 433 (1988).

17. U. F. Kocks, A. S. Argon, and M. F. Ashby, Thermodynamics and Kinetics of Slip, Prog. Mater. Sci., Vol. 19, Pergamon Press, Oxford, England (1975).

18. M. B. M. Mangion, J. Y. Cavaille, and J. Perez, Phil. Mag. A, 66, 773 (1992).

19. P. H. Mott, PhD thesis, MIT (1992).

20. D. Stauffer and A. Aharony, An Introduction to Percolation Theory, Taylor & Francis, London (1992).

21. Y. Jean, Microchem J., 42, 72 (1990).

22. X. S. Li, private communication on distribution of local free volume measurements, 1993.