Viscoelasticity and viscoplasticity of semicrystalline polymers: Structure–property relations for high-density polyethylene

Computational Materials Science 39 (2007) 729–751 www.elsevier.com/locate/commatsci

Viscoelasticity and viscoplasticity of semicrystalline polymers: Structure–property relations for high-density polyethylene A.D. Drozdov a

a,*

, J. deC. Christiansen

b

Department of Chemical Engineering, Ben-Gurion University of the Negev, P.O. Box 653, Beer-Sheva 84105, Israel b Department of Production, Aalborg University, Fibigerstraede 16, DK-9220 Aalborg, Denmark Received 7 August 2006; accepted 18 September 2006

Abstract Observations are reported on two commercial grades of high-density polyethylene (HDPE) in uniaxial tensile tests, relaxation tests, creep tests and cyclic tests with a strain-controlled deformation program. Constitutive equations are derived for the viscoelastic and viscoplastic responses of semicrystalline polymers at three-dimensional deformation with small strains. A polymer is modeled as a two-phase continuum consisting of a crystalline skeleton and an amorphous phase treated as a transient network of chains. Its viscoelastic response is associated with thermally activated rearrangement of strands in the temporary network. The viscoplastic behavior reflects fine and coarse slip of lamellar stacks and sliding of junctions between chains in the network. Adjustable parameters in the stress–strain relations are found by fitting the experimental data. The study focuses on the effect of molecular weight of HDPE on its mechanical properties.  2006 Elsevier B.V. All rights reserved. Keywords: High-density polyethylene; Viscoelasticity; Viscoplasticity; Cyclic deformation

1. Introduction This paper is concerned with the experimental investigation and numerical simulation of the viscoelastic and viscoplastic responses of semicrystalline polymers. Correlations between crystalline morphology and mechanical properties of polymers have attracted substantial attention in the past decade, see reviews [1–4] and the references therein. Previous studies of structure–property relations focused mainly on structure, whereas mechanical properties remained of secondary importance (the authors confined themselves to the evaluation of the Young’s modulus, yield stress and elongation to break). Not arguing the significance of these works for characterization of semicrystalline polymers, we would like to emphasize their shortcoming: the lack of a reliable model that allows changes in the time- and rate-

*

Corresponding author. Tel.: +972 8 647 2146. E-mail address: [email protected] (A.D. Drozdov).

0927-0256/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2006.09.001

dependent behavior to be predicted quantitatively when internal structure of a polymer is altered. Development of structure–property relations for semicrystalline polymers is a complicated task due to insufficient knowledge about • the influence of internal structure of chains (concentration of entanglements, content and length of side chains, etc.) on crystalline morphology and mechanical properties of crystallites at the stage of cooling of a melt, • the effect of crystalline morphology and mechanical properties of crystallites on the viscoelastic and viscoplastic responses of a polymer at ambient temperature. To avoid (at least, partially) these complications, two commercial grades of high-density polyethylene (HDPE) with different molecular weights are chosen for the experimental investigation. As HDPE contains linear chains with a small amount of side-chains, we suppose that mass

730

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

average molecular weight may be considered as the only internal parameter responsible for the mechanical properties of amorphous and crystalline phases. This is a rather crude approximation, which implies that the crystalline morphology is not affected by distribution of chains with various lengths. This approach is adopted with the only aim to simplify interpretation of results. Samples for mechanical tests were injection-molded under relatively close thermal conditions, which implies that the effect of processing on mechanical properties may be neglected. For a discussion of the influence of processing conditions on the crystalline morphology of injection-molded HDPE, the reader is referred to [5,6]. We do not dwell on characterization of the crystalline structure of HDPE. This may be explained by the following reasons: 1. The effect of morphology of a polymer network on the crystalline structure of polyethylene was investigated in a number of publications [7–10]. 2. Although some correlations were established between morphology of a semicrystalline polymer and its elastic and plastic deformations, no quantitative relations are available. Even prediction of an elastic modulus based on the theory of two-phase composite media remains a complicated and non-reliable procedure [11]. Unlike conventional studies on the mechanical response of HDPE that focused mainly on uniaxial tension with constant cross-head speeds [12–16], our experimental program includes (i) uniaxial tension with various cross-head speeds, (ii) standard relaxation tests at various strains, (iii) standard creep tests at various stresses, and (iv) cyclic tensile tests with a strain-controlled program and various maximum strains. This program allows the viscoelastic and viscoplastic behavior to be analyzed at active loading (uniaxial tension), neutral loading (creep), unloading and reloading (cyclic tests). The viscoelastic response of high-density polyethylene was studied in [17–21] (tensile creep), [22] (tensile and compressive creep), [23–26] (tensile relaxation), [20,27] (dynamic tests with small strains) and [13,18,28] (strain recovery). The effect of strain rate on the stress–strain curves of HDPE was investigated in [14–16,18,19,21,25,26,28] (uniaxial tension), [29,30] (uniaxial compression), [31] (simple shear) and [9,10,32] (plain strain compression). Cyclic deformation of high-density polyethylene was analyzed in [13,21,28]. Stress–strain relations for the viscoelastic and viscoplastic responses of HDPE were developed in [13–15,19,21, 24–26,29,33–35]. Among constitutive equations for semicrystalline polymers that may be applied to modeling the time- and rate-dependent behavior of high-density polyethylene, it is worth mentioning [36–48]. Two shortcomings of the above relations (or, at least, most of them) are to be emphasized:

1. When material constants in the governing equations are found by fitting observations in standard relaxation tests, they fail to predict adequately the experimental data in creep tests (in particular, when the creep tests are performed at stresses relatively close to the apparent yield stress). 2. The constitutive equations do not reproduce correctly observations in cyclic tests with more than one cycle of loading–retraction. The objective of the present study is four-fold: 1. To report experimental data in uniaxial tensile tests, relaxation tests, creep tests and cyclic tests on two commercial grades of HDPE. 2. To derive a constitutive model for the viscoelastic and viscoplastic responses of a semicrystalline polymer at arbitrary three-dimensional deformations with small strains and to find adjustable parameters by matching the observations. 3. To discuss the effect of molecular weight of HDPE on the material constants in the stress–strain relations. 4. To compare results of numerical simulation in tests with sophisticated time-dependent deformation programs with observations on semicrystalline polymers. Following [7,14,21,24], a semicrystalline polymer is thought of as a two-phase continuum formed by a crystalline skeleton surrounded by meso-regions of chains in the rubbery state. The crystalline phase is modeled as an isotropic viscoplastic medium. The plastic flow of crystallites reflects fine (homogeneous shear of crystal blocks) and coarse (heterogeneous inter-lamellar sliding) slip of lamellar blocks. With reference to [24,29], the amorphous phase is treated as a temporary network of strands bridged by junctions (entanglements and physical cross-links on the lamellar surfaces). The network is assumed to be strongly heterogeneous [38,39,41]. Its inhomogeneity is attributed to the heterogeneity of interactions between chains in the amorphous phase and lamellar blocks with various lengths and thicknesses (chains in the rubbery state are located between spherulites and inside spherulites between lamellae stacks). Following [49], the viscoelastic response is associated with thermally activated rearrangement of strains in the network (separation of active strands from their junctions and attachment of dangling strands to the network). The rate of detachment is governed by the Eyring equation [50], where different meso-regions are characterized by different activation energies. Viscoplasticity of the amorphous phase is thought of as flow of junctions in the network with respect to their reference positions driven by plastic deformation of surrounding crystallites. To simplify derivation of the stress–strain relations, the crystalline and amorphous phases are treated as incompressible media. The incompressibility hypothesis reflects the fact that high-density polyethylene is a weakly com-

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

pressible polymer (its Poisson’s ratio is conventionally estimated as m = 0.41 [28] or m = 0.42 [51]). Two starting points in the development of constitutive equations are to be mentioned. The first is the concept of viscoplasticity introduced in [41], according to which the strain-rate tensor for plastic deformations is proportional to the strain-rate tensor for macro-deformation (not to the stress tensor as conventional theories of plasticity presume). The other is the concept of pseudo-elasticity [52], which postulates that some parameters in the stress– strain relations may be considered as piece-wise constant functions of strain: these parameters remain constant along each loading and retraction path and change these values when the sign of strain rate is altered. We treat the coefficient of proportionality between the rate of strain tensors for plastic deformation and macro-deformation as the only parameter with such a property. This implies that not only the strain tensor for plastic deformation is a continuous function of time (as conventional approaches require), but its derivative with respect to time as well. The exposition is organized as follows. Observations on two grades of HDPE with different molecular weights are reported in Section 2. Constitutive equations for the viscoelastic and viscoplastic responses of a semicrystalline polymer are developed in Section 3. Adjustable parameters in the stress–strain relations are found in Section 4 by fitting the experimental data. The effect of molecular weight on the mechanical response is briefly discussed in Section 5. Some results of numerical simulation are presented in Section 6. Concluding remarks are formulated in Section 7.

engineering stress r was determined as the ratio of the axial force to the cross-sectional area of specimens in the stressfree state. The first series of experiments involved two tensile tests with constant strain rates. The specimens were deformed with constant cross-head speeds of 10 and 50 mm/min (which corresponded to the strain rates _ ¼ 0:002 and _ ¼ 0:01 s1 ) up to the maximum strain max = 0.1. The chosen strain rates ensure nearly isothermal experimental conditions. No necking of samples was observed. Given a strain rate _ , at least three tests were performed on different specimens. The experimental data demonstrate good reproducibility of measurements. The stress–strain diagrams (the engineering stress r versus engineering strain ) are plotted in Fig. 1 (HDPE-A) and Fig. 9 (HDPE-B). These stress–strain curves are strongly nonlinear. The stress r monotonically increases with  (the yield strains of both grades of polyethylene exceed 0.1). Given a strain , the engineering stress noticeably grows with strain rate _ . For a fixed , the longitudinal stress of HDPE-B substantially exceeds that of HDPE-A. The other series of experiments consisted of four standard relaxation tests with the strains  = 0.025, 0.05, 0.075 and 0.10 for HDPE-A and three relaxation tests with the strains  = 0.025, 0.05, 0.075 for HDPE-B. In each relaxation test, a specimen was stretched with a constant strain rate _ ¼ 0:01 s1 up to a given longitudinal strain  that was preserved constant at relaxation. In accord with the ASME protocol for short-term relaxation tests, the relaxation time tr = 20 min was chosen. The dimensionless tensile stress

2. Experimental procedure ¼ r High-density polyethylene Eraclene MM 95 (density 0.953 g/cm3, melt flow index 4 g/10 min [190 C, 2.16 kg], melting temperature Tm = 134 C) was supplied by Polimeri Europa SpA (Italy). Its mass average molecular weight is roughly estimated as 100 kg/mol based on the results for a similar polymer (Eraclene MP 90) reported in [53]. We treat this polymer as HDPE with a low molecular weight and designate it as HDPE-A. High-density polyethylene Lupolen 5261 Z (density 0.954 g/cm3, melt flow rate 0.1 g/10 min [190 C, 21.6 kg], melting temperature Tm = 134 C) was purchased from Basell Polyolefins (Zaventem, Belgium). Its mass average molecular weight is estimated as 420 [54] to 450 kg/mol [55]. This polymer is treated as HDPE with a high molecular weight, and it is designated as HDPE-B. Dumbbell specimens for uniaxial tensile tests (ASTM standard D638) with length 140 mm, width 9.85 mm and thickness 3.74 mm were prepared by using injection-molding machine Ferromatic K110/S60-2 K. Mechanical tests were conducted at room temperature with the help of a universal testing machine Instron-5568 equipped with electro-mechanical sensors for the control of longitudinal strains in an active zone of samples. The tensile force was measured by a standard load cell. The

731

r ; r0

ð1Þ

Fig. 1. The engineering stress r versus engineering strain  in tensile tests with various cross-head speeds. Symbols: experimental data on HDPE-A. Solid lines: results of numerical simulation.

732

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

where r0 is the engineering stress at the beginning of the relaxation process, is plotted versus the logarithm (log = log10) of time t (the initial instant t = 0 corresponds to the beginning of relaxation) in Figs. 2 (HDPE-A) and 10 (HDPE-B). These figures show that the stress r monotonically decreases with time, and the rate of stress relaxation is practically independent of strains  at which the tests were conducted. The third series of experiments involved four creep tests with the engineering stresses r = 10, 15, 17.5 and 20 MPa on HDPE-A and the stresses r = 20, 30, 35 and 40 MPa on HDPE-B. These stresses (approximate values are provided) were chosen to characterize the time-dependent response of polyethylenes in the interval of strains where the mechanical response is strongly nonlinear. In each creep test, a specimen was stretched with a constant strain rate _ ¼ 0:01 s1 up to a given engineering stress r that was preserved constant during the interval of creep. With reference to the ASME protocol for short-term creep tests, the creep time tc = 20 min was set. At relatively high engineering stresses, the creep tests were stopped when the strain  exceeded 0.1, to avoid necking of specimens. The strain  is plotted versus the logarithm of time t (the initial instant t = 0 corresponds to the beginning of creep) in Figs. 3 (HDPE-A) and 11 (HDPE-B). These figures show that given r, the strain  monotonically increases with time. For both grades of HDPE, the rate of increase in strain is strongly affected by r and it grows with engineering stress. Our results are in agreement with the observations in tensile creep tests on HDPE reported in Fig. 6 of [18] and Fig. 1 of [19]. They differ, however, from the creep curves depicted in Fig. 9 of [21] (which show that the strain  linearly increases with log t).

 versus time t in tensile relaxation tests at Fig. 2. The dimensionless stress r various strains . Symbols: experimental data on HDPE-A. Solid line: results of numerical simulation.

Fig. 3. The engineering strain  versus time t in tensile creep tests with various engineering stresses r (MPa). Symbols: experimental data on HDPE-A. Solid lines: results of numerical simulation.

The last series of tests consisted of four cyclic tests (a strain-controlled program) with the maximum strains max = 0.025, 0.05, 0.075 and 0.1 (HDPE-A) and max = 0.03, 0.05, 0.07 and 0.09 (HDPE-B). In each cyclic test, a specimen was loaded up to a maximum strain max with the strain rate _ ¼ 0:002 s1 , unloaded down to the zero stress with the strain rate _, reloaded up to the maximum strain with the strain rate _ , etc. Each test consisted of N = 10 cycles of deformation. To avoid buckling of samples, the minimum stresses at retraction rmin = 0.4 MPa and rmin = 0.04 MPa were chosen for HDPE-A and HDPE-B, respectively. The stress–strain diagrams (the engineering stress r versus engineering strain ) for the first loading, first retraction and first reloading are depicted in Figs. 4 (HDPE-A) and 12 (HDPE-B). Fig. 4 demonstrates that the hysteresis loop of HDPE-A is strongly affected by the maximum strain max. With an increase in max, the slope of the hysteresis curve decreases, whereas its area noticeably grows. On the contrary, Fig. 12 shows that the slope of the hysteresis curve of HDPE-B remains practically constant, whereas its area increases with max rather weakly. Similar observations (but in stress-controlled cyclic tensile tests with finite strains) were reported in Fig. 9 of [9]. The difference between our results and those of Ref. [9] is that the strong effect of maximum strain on the area of hysteresis curves is found for HDPE-A with a relatively small molecular weight (and, as a consequence, with a relatively low concentration of entanglements), whereas the same conclusion was reached in [9] for an annealed and quenched polyethylene, which was presumed to have a higher number of entanglements than another polyethylene crystallized under high pressure. This disagreement may be

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

733

in stress-controlled cyclic compressive tests reported in [34]. Fig. 7 of [34] shows that in increase in degree of crystallinity of HDPE prevents changes in the hysteresis curves measured at various maximum strains. Detailed dependencies of the engineering stress r on engineering strain  in cyclic tests with various maximum strains max are plotted in Figs. 5–8 (HDPE-A) and 13–16 (HDPE-B). These figures show that

Fig. 4. The engineering stress r versus engineering strain  in cyclic tensile tests with various maximum strains max and the minimum stress rmin = 0.4 MPa. Symbols: experimental data on HDPE-A. Solid lines: results of numerical simulation.

resolved, however, if we suppose that concentration of entanglements in a melt affects the viscoplastic response of a semicrystalline polymer not only directly, but through changes in the structure of crystallites as well, in other words, if we presume mechanical properties of spherulites in quenched and slowly crystallized under high pressure polyethylenes (used for testing in Ref. [9]) to distinguish substantially. The latter is in agreement with observations

Fig. 5. The engineering stress r versus engineering strain  in a cyclic tensile test with the maximum strain max = 0.025 and the minimum stress rmin = 0.4 MPa. Circles: experimental data on HDPE-A. Solid line: results of numerical simulation.

Fig. 6. The engineering stress r versus engineering strain  in a cyclic tensile test with the maximum strain max = 0.05 and the minimum stress rmin = 0.4 MPa. Circles: experimental data on HDPE-A. Solid line: results of numerical simulation.

Fig. 7. The engineering stress r versus engineering strain  in a cyclic tensile test with the maximum strain max = 0.075 and the minimum stress rmin = 0.4 MPa. Circles: experimental data on HDPE-A. Solid line: results of numerical simulation.

734

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

the plastic strain tensor ^p characterizes (i) rotation and twist of individual lamellae and lamellar stacks in spherulites, (ii) fine (homogeneous shear of crystal blocks) and coarse (heterogeneous inter-lamellar sliding) slip of lamellar blocks, and (iii) micro-necking of lamellae. The rate-of-strain tensor for plastic deformation is proportional to the rate-of-strain tensor for macrodeformation, d^p d^ ¼/ ; dt dt

^p ð0Þ ¼ ^0;

ð3Þ

where t stands for time, and / is a scalar function to be described in what follows. The initial condition in Eq. (3) (where ^0 denotes the zero tensor) expresses the fact that plastic strain vanishes in a stress-free medium. The strain energy of the crystalline phase (per unit volume of a semicrystalline polymer) is given by 1 W c ¼ lc^e : ^e ; 2 Fig. 8. The engineering stress r versus engineering strain  in a cyclic tensile test with the maximum strain max = 0.1 and the minimum stress rmin = 0.4 MPa. Circles: experimental data on HDPE-A. Solid line: results of numerical simulation.

1. The experimental stress–strain diagrams at subsequent reloadings and retractions are strongly nonlinear. 2. The apparent residual strain (the strain measured at the instant when the stress r vanishes at retraction) increases with number of cycles. The most pronounced growth of the apparent residual strain occurs during the first three cycles of deformation. 3. The maximum stress per cycle (the stress measured at the instant when  = max) and the area of hysteresis curves monotonically decrease with number of cycles. 3. Constitutive model For the sake of generality, stress–strain relations for the viscoelastic and viscoplastic responses of a semicrystalline polymer are developed for an arbitrary three-dimensional deformation with small strains. The polymer is thought of as a two-phase continuum consisting of crystalline and amorphous phases. The crystalline phase is treated as a skeleton of coupled crystal blocks, whereas the amorphous phase is modeled as a network of entangled chains [24]. 3.1. Crystalline phase The crystalline skeleton is treated as an incompressible viscoplastic medium whose deformation coincides with macro-deformation of the polymer. At small strains, the strain tensor ^ is split into the sum of strain tensors for elastic, ^e , and plastic, ^p , deformations, ^ ¼ ^e þ ^p :

ð2Þ To ensure the incompressibility condition, the tensors ^e and ^p are assumed to be traceless. At the micro-level,

ð4Þ

where lc stands for an elastic modulus of the crystalline phase, and colon denotes convolution of tensors. Differentiating Eq. (4) with respect to time and using Eqs. (2) and (3), we find that dW c d^ ðtÞ ¼ lc ð1  /ðtÞÞð^ðtÞ  ^p ðtÞÞ : ðtÞ: dt dt

ð5Þ

3.2. Amorphous phase Following [29], the amorphous phase is treated as a temporary network of flexible chains bridged by junctions (entanglements between chains in the rubbery state and physical cross-links on the lamellar surfaces). With reference to the concept of temporary networks [49], we suppose that active strands (whose ends are connected to contiguous junctions) separate from these junctions at random times when the strands are thermally activated. An active strand whose end detaches from a junction is transformed into a dangling strand. A dangling strand returns into the active state when its free end captures a nearby junction at a random instant. At the micro-level, rearrangement of strands reflects (i) disentanglement of chains in the amorphous phase and (ii) detachment of chain folds and loops from surfaces of lamellar blocks. The network of chains is modeled as an ensemble of meso-regions with various activation energies v for separation of strands from their junctions. According to the theory of thermally activated processes [50], the rate of detachment in a meso-region with an activation energy v reads   v C ¼ c exp  ; ð6Þ kBT where kB is Boltzmann’s constant, T is the absolute temperature, and the pre-factor c is independent of energy v and temperature T. The activation energies of meso-regions v

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

and the attempt rate c for transition of active strands into the dangling state are independent of mechanical factors. Confining ourselves to isothermal deformations at some temperature T0 and introducing the dimensionless activation energy v ¼ v=ðk B T 0 Þ, we present Eq. (6) in the form CðvÞ ¼ c expðvÞ:

ð7Þ

Denote by N the number of active strands per unit volume of a semicrystalline polymer and by n0(v) the number of active strands (per unit volume) in meso-regions with dimensionless energy v. The distribution function p(v) for active meso-regions with various activation energies v is given by pðvÞ ¼

n0 ðvÞ : N

ð8Þ

This function is independent of mechanical factors and satisfies the normalization condition Z 1 pðvÞ dv ¼ 1: ð9Þ 0

To approximate experimental data, we adopt the random energy model [56] with the quasi-Gaussian distribution of meso-regions " # 2 ðv  V Þ pðvÞ ¼ p0 exp  ðv P 0Þ; pðvÞ ¼ 0 ðv < 0Þ; 2R2 ð10Þ where V is an apparent average activation energy for separation of strands from temporary junctions, R is an apparent standard deviation of activation energies, and the pre-factor p0 is determined by Eq. (9). Rearrangement of a temporary network is described by the function n(t, s, v) that equals the number of active strands at time t (per unit volume) in meso-regions with activation energy v that have last rearranged before instant s 2 [0, t]. In particular, n(t, t, v) is the number (per unit volume) of active strands in meso-regions with activation energy v at time t, nðt; t; vÞ ¼ n0 ðvÞ; and the amount u(s, v) ds, where   on uðs; vÞ ¼ ðt; s; vÞ os t¼s

ð11Þ

ð12Þ

equals the number (per unit volume) of dangling strands in meso-regions with activation energy v that merge with the network within the interval [s, s + ds]. The rate of detachment C equals the ratio of the number of active strands that separate from the network per unit time to the current number of active strands. Applying this definition to active strands that separate (for the first time) from temporary junctions within the interval [t, t + dt] and to those that merged with the network during the interval [s, s + ds] and separate from their junctions within the interval [t, t + dt], we find that

on ðt; 0; vÞ ¼ CðvÞnðt; 0; vÞ; ot

735

o2 n on ðt; s; vÞ ¼ CðvÞ ðt; s; vÞ: ot os os ð13Þ

The solutions of Eqs. (13) with initial conditions (11) and (12) read nðt; 0; vÞ ¼ n0 ðvÞ exp½CðvÞt; on ðt; s; vÞ ¼ uðs; vÞ exp½CðvÞðt  sÞ: os

ð14Þ

The function u(t, v) is found from Eqs. (14) and the identity Z t on ðt; s; vÞ ds: ð15Þ nðt; t; vÞ ¼ nðt; 0; vÞ þ 0 os Simple algebra implies that uðt; vÞ ¼ CðvÞn0 ðvÞ:

ð16Þ

Insertion of expressions (8) and (16) into Eqs. (14) results in nðt; 0; vÞ ¼ NpðvÞ exp½CðvÞt; on ðt; s; vÞ ¼ NpðvÞCðvÞ exp½CðvÞðt  sÞ: os

ð17Þ

Separation of active strands from their junctions and attachment of dangling strands to the network describe the viscoelastic behavior of the amorphous phase. Its viscoplastic response is associated with sliding of junctions with respect to their reference positions. At the micro-level, sliding of junctions reflects (i) chain slip through the crystals and (ii) sliding of tie chains along the surfaces of lamellar blocks. At small strains, we write by analogy with Eq. (2), ^ ¼ ^ee þ ^ep ;

ð18Þ

where the strain tensor ^ee describes elastic deformation of the amorphous phase, and the strain tensor ^ep characterizes viscoplastic flow of junctions with respect to their reference positions (i.e., their positions in a stress-free network). The difference between Eqs. (2) and (18) is that for a homogeneous macro-deformation, all tensors in Eq. (2) depend on time t only, whereas the tensors on the righthand side of Eq. (18) are functions of two times: the current instant t and the instant s 2 [0, t] when appropriate active strands merged with the network. In the rigorous form, Eq. (18) reads ^ðtÞ ¼ ^e0e ðtÞ þ ^e0p ðtÞ;

^ðtÞ ¼ ^ee ðt; sÞ þ ^ep ðt; sÞ;

ð19Þ

where ^e0e ðtÞ and ^e0p ðtÞ are strain tensors for elastic and viscoplastic deformations of strands that have not separated from the network within the interval [0, t], while ^ee ðt; sÞ and ^ep ðt; sÞ are strain tensors for elastic and viscoplastic deformations of strands that have attached to the network at instant s 2 [0, t]. We suppose that ^e0p ðtÞ and ^ep ðt; sÞ are independent of v, which means that the viscoplastic flow of junctions occurs in the same way in meso-regions with different activation energies for separation of strands. The strain energy stored in an active strand is determined by analogy with Eq. (4),

736

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

1 w0 ðtÞ ¼ ls^e0e ðtÞ : ^e0e ðtÞ; 2

1 wðt; sÞ ¼ ls^ee ðt; sÞ : ^ee ðt; sÞ; 2

where ls stands for an elastic modulus per strand. The strain energy of the amorphous phase Wa (per unit volume of a semicrystalline polymer) equals the sum of strain energies of individual strands (following the conventional approach, the energy of interactions between strands is accounted for by means of the incompressibility condition only). Multiplying the strain energy per strand by the number of active strands (per unit volume) in meso-regions with various activation energies v, we arrive at the formula for the strain energy density of the amorphous phase Z 1 1 W a ðtÞ ¼ ls nðt; 0; vÞ^e0e ðtÞ : ^e0e ðtÞ 2 0  Z t on ðt; s; vÞ^ee ðt; sÞ : ^ee ðt; sÞ ds dv: þ ð20Þ 0 os ^e0p ðtÞ

To introduce kinetic equations for the strain tensors and ^ep ðt; sÞ, it is postulated that flow of junctions in the amorphous phase is entirely determined by viscoplastic deformation of the crystalline skeleton. From the physical standpoint, this assumption means that plastic deformation of amorphous meso-regions located between lamellar blocks is entirely determined by fine and coarse slip of these blocks. It follows from this hypothesis that d^e0p d^ ¼/ ; dt dt

o^p d^ ¼/ ; dt ot

ð21Þ

where / is the same function as in Eq. (3). The initial conditions for Eqs. (21) read  ^e0p ð0Þ ¼ ^ 0; ^ep ðt; sÞt¼s ¼ ^ðsÞ: ð22Þ The first equality in Eqs. (22) means that no slip of junctions with respect to their reference positions occurs in a stress-free network. The other equality follows from Eq. (19), where we set t = s and keep in mind that stresses totally relax in dangling chains before their attachment to the network [49] (the latter implies that ^ee ðs; sÞ ¼ 0). Integrating Eqs. (21) with initial conditions (22) and using Eq. (3), we arrive at the formulas

Z 1 dW a ðtÞ ¼ Y ðtÞ þ ls ð1  /ðtÞÞ nðt;0;vÞð^ðtÞ ^p ðtÞÞ dt 0  Z t on d^ ðt; s;vÞ½ð^ðtÞ ^p ðtÞÞ  ð^ðsÞ ^p ðsÞÞds dv : ðtÞ; þ dt 0 os

ð25Þ where ( Z 1 1 CðvÞ nðt; 0; vÞð^ðtÞ  ^p ðtÞÞ : ð^ðtÞ  ^p ðtÞÞ Y ðtÞ ¼ ls 2 0 Z t on ðt; s; vÞ½ð^ðtÞ  ^p ðtÞÞ  ð^ðsÞ  ^p ðsÞÞ þ 0 os ) : ½ð^ðtÞ  ^p ðtÞÞ  ð^ðsÞ  ^p ðsÞÞds dv:

ð26Þ

3.3. Stress–strain relations The strain energy density per unit volume of a semicrystalline polymer W equals the sum of the strain energy densities of the amorphous and crystalline phases, W ðtÞ ¼ W a ðtÞ þ W c ðtÞ:

ð27Þ

At isothermal deformation of an incompressible medium with small strains, the Clausius–Duhem inequality reads QðtÞ ¼ 

dW d^ ^0 ðtÞ : ðtÞ P 0; ðtÞ þ r dt dt

ð28Þ

where Q stands for energy dissipation per unit time and ^ is the stress tensor, and prime denotes the unit volume, r deviatoric component of a tensor. Substitution of Eqs. (5), (25) and (27) into Eq. (28) results in  ^0 ðtÞ  ð1  /ðtÞÞ ðlc þ ls N Þð^ðtÞ  ^p ðtÞÞ QðtÞ ¼ Y ðtÞ þ r  Z 1 Z t on d^ ðt; s; vÞð^ðsÞ  ^p ðsÞÞds ls dv : ðtÞ; ot dt 0 0 ð29Þ

Insertion of expressions (19) and (23) into Eq. (20) results in Z 1 1 W a ðtÞ ¼ ls nðt; 0; vÞð^ðtÞ  ^p ðtÞÞ : ð^ðtÞ  ^p ðtÞÞ 2 0 Z t on ðt; s; vÞ½ð^ðtÞ  ^p ðtÞÞ  ð^ðsÞ  ^p ðsÞÞ þ 0 os  : ½ð^ðtÞ  ^p ðtÞÞ  ð^ðsÞ  ^p ðsÞÞds dv: ð24Þ

where we used Eqs. (8), (9), (11) and (15). It follows from Eq. (26) that the function Y(t) is non-negative: the rate of detachment of active strands from their junctions C(v) is positive in accord with Eq. (7), the functions n(t, 0, v) and on/os(t, s, v) are positive according to Eqs. (17), whereas convolution of any tensor with itself is non-negative. Eq. (29) implies that the Clausius–Duhem inequality is satis^ is given by fied, provided that the stress tensor r  ^ðtÞ ¼ P ðtÞ^I þ ð1  /ðtÞÞ ðlc þ ls N Þð^ðtÞ  ^p ðtÞÞ r  Z 1 Z t on ðt; s; vÞð^ðsÞ  ^p ðsÞÞds ; ls dv ð30Þ 0 0 ot

Differentiating Eq. (24) with respect to time and using Eqs. (3) and (13), we find that

where P is an unknown pressure, and bI stands for the unit tensor. Substituting expression (17) into Eq. (30) and introducing the notation

^e0p ðtÞ ¼ ^p ðtÞ;

^ep ðt; sÞ ¼ ^p ðtÞ  ^p ðsÞ þ ^ðsÞ:

ð23Þ

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

l ¼ lc þ ls N ;

j¼

ls N ; lc þ ls N

ð31Þ

we arrive at the constitutive equation for the viscoelastic and viscoplastic responses of a semicrystalline polymer " b ^ðtÞ ¼ P ðtÞ I þ lð1  /ðtÞÞ ð^ðtÞ ^p ðtÞÞ r j

Z

Z

1

pðvÞCðvÞdv 0

#

t

expðCðvÞðt  sÞÞð^ðsÞ ^p ðsÞÞds :

0

ð32Þ Eqs. (3) and (32) are satisfied for an arbitrary isothermal deformation program with small strains. 3.4. Evolution of the function / Viscoplastic flow is characterized by the function /(t) in the stress–strain equations. To describe this function, three regimes of deformation are distinguished: 1. Active loading of a stress-free medium, when intensity of plastic strains grows with time. The stress intensity may increase (in uniaxial tensile tests below the yield point) or decrease (in tensile tests above the yield point or in relaxation tests which are treated as active loading with the zero strain rate). 2. Neutral loading. In experiments, this regime corresponds to creep tests. It is characterized by the fact that the stress tensor remains constant, and the stress intensity equals the maximum stress intensity reached at active loading. The intensity of plastic strains increases with time due to viscoplastic flow. 3. Unloading and reloading. In experiments, this regime corresponds to subsequent reloadings and retractions in cyclic tensile tests. It is characterized by the fact that the stress intensity and the intensity of plastic strain change with time.

737

stands for intensity of elastic strains, and a and b are positive parameters. Our choice of the stretched exponential function (33) is explained by the following reasons: 1. Formula (33) involves only two material constants a and b. 2. This relation is rather flexible to ensure good agreement with observations for a number of semicrystalline polymers. 3. According to Eqs. (3) and (33), the rate of plastic strain equals zero at the initial instant, it monotonically increases with elastic deformation, and reaches the strain rate for macro-deformation when the elastic strain becomes relatively large. All these conclusions appear to be physically plausible. 4. A similar function (although in another context) was previously employed to describe the mechanical response of solid polymers [57] (the so-called Matsuoka equation). The exponent b 2 (0, 1] in Eq. (33) serves as a measure of nonlinearity of the mechanical behavior, whereas the prefactor a characterizes an elastic strain, at which transition occurs to a developed viscoplastic flow. Evolution of the plastic strain tensor ^p at active loading with a constant strain rate (Figs. 1 and 9) is described by Eqs. (3) and (33). These relations are noticeably simplified for relaxation tests. According to Eq. (3), in the latter case, the strain tensor for plastic deformations remains timeindependent, d^p ¼ 0: dt

ð34Þ

This classification of deformation processes slightly differs from the standard classification of active loading, neutral loading and unloading of an elastoplastic medium with a yield surface. A reason for our modification is to account for (i) the viscoelastic response (rearrangement of strands in a transient network) and (ii) evolution of hysteresis curves at cyclic loading. 3.4.1. Active loading At active loading of a virgin polymer, the function / monotonically increases with time following the phenomenological relation / ¼ 1  expðaJ be Þ; where rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ^e : ^e Je ¼ 3

ð33Þ Fig. 9. The engineering stress r versus engineering strain  at tension with various cross-head speeds. Symbols: experimental data on HDPE-B. Solid lines: results of numerical simulation.

738

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

3.4.2. Neutral loading At neutral loading, the function / reads / ¼ /0 þ D/;

Eq. (42) is chosen for the following reasons: ð35Þ

where /0 equals the value of / reached along the loading path, whereas the increment D/ is determined by the formula similar to Eq. (33), b

D/ ¼ 1  exp½AðDJ e Þ :

ð36Þ

Here A is a positive adjustable parameter, b is the same exponent as in Eq. (33), DJ e ¼ J e  J e0

ð37Þ

is the increment of strain intensity for elastic deformation, and Je0 is the intensity of elastic strain reached at active loading. The following explanation may be provided for our choice of Eqs. (35) and (36). Differentiation of Eq. (33) with respect to Je implies that at active loading, the function / obeys the nonlinear differential equation d/ ¼ abð1  /ÞJ eb1 ; dJ e

/ð0Þ ¼ 0:

ð38Þ

Eqs. (35) and (36) mean that after transition to neutral loading, the function / proceeds to grow with DJe, and evolution of this function is described by the kinetic equation similar to Eq. (38), d/ b1 ¼ Abð1  /ÞðDJ e Þ : dDJ e

ð39Þ

The initial condition for Eq. (39) reflects the continuity condition for the function / at transition from active to neutral loading. The coefficient A in Eq. (36) depends on the entire history of deformation up to an instant t0, when active loading is replaced with neutral loading, and remains constant afterward. The effect of deformation history on A is described by the phenomenological relation log A ¼ A0  A1 log W p ; where A0 and A1 are material parameters, and Z t0 d^p ^0 ðtÞ : ðtÞ dt r Wp ¼ dt 0

ð40Þ

ð41Þ

stands for the work of external forces on plastic deformations at active loading.

1. It contains only two adjustable parameters a and b. 2. Eq. (42) involves the derivative (with respect to time) of the strain energy Wc of the crystalline phase as an input, which seems physically plausible: evolution of the rate of plastic strain is determined by the rate of changes in mechanical energy stored in the crystalline skeleton under deformation. 3. Formula (42) ensures that the current rate of plastic strain is affected by the entire history of deformation. 4. Eq. (42) implies that the rate of changes in / is higher at retraction than at active loading. Formula (42) differs from the kinetic equation introduced in our previous study on the viscoplastic response of semicrystalline polymers [58] by the presence of modulus on the right-hand side and alternation of the sign of the second term in the left-hand side at loading and retraction. Insertion of expression (5) into Eq. (42) results in    d/ d^   ^ ð43Þ  b/ ¼ að1  /ðtÞÞð^ðtÞ  p ðtÞÞ : ðtÞ: dt dt Following the concept of pseudo-elasticity [52], we treat a and b in Eq. (43) as constants along each path of the stress–strain diagram at cyclic loading. This means that a and b remain constant during each subsequent retraction and reloading, but they alter their values at the instants when the strain rate changes its sign. To complete description of the constitutive model, the following parameters are to be determined: 1. a ¼ a0þ and b ¼ b0þ for the first retraction path, that is, the quantities a and b at the first instant when the strain tensor reaches its ‘‘maximum value’’ and the strain rate changes its sign, 2. a ¼ a0 and b ¼ b0 for the first reloading path, i.e., the quantities a and b at the first instant when the stress ten^ reaches its ‘‘minimum value’’ and the strain rate sor r changes its sign, 3. a rule that allows the parameters a0þ , b0þ and a0 , b0 to be transformed into the coefficients a and b in Eq. (43) for each subsequent cycle of deformation. By analogy with Eq. (40), we set

3.4.3. Cyclic loading At subsequent reloading and retraction paths at cyclic deformation, evolution of the function /(t) is governed by the first-order kinetic equation   d/ a dW c  ; ð42Þ  b/ ¼  dt lc dt  where a and b are positive quantities, and the second term in the left-hand side is positive for active loading ð^ r0 : d^=dt > 0Þ and negative at retraction ð^ r0 : d^=dt < 0Þ.

log a0þ ¼ a00þ  a01þ log W p ;

log b0þ ¼ b00þ  b01þ log W p ;

log a0 ¼ a00  a01 log W p ;

log b0 ¼ b00  b01 log W p ; ð44Þ

where Wp is the work of external forces on plastic deformations until the instant when the first retraction (first reloading) starts, and a0i , b0i (i = 0, 1) are material constants. To describe evolution of the coefficients a and b in Eq. (43) with number of cycles, the following phenomenological relations are proposed for retraction and reloading:

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

a ¼ k þ a0þ ;

b ¼ k þ b0þ ;

k  a0 ;

k  b0 :

a¼

b¼

ð45Þ

Eqs. (45) mean that the coefficients a and b are affected by the deformation history in the same fashion. Insertion of Eqs. (45) into Eq. (43) implies that    d/ d^ 0 0  ð46Þ  b / ¼ a ð1  /Þð^  ^p Þ : ; dt dt where a0 and b0 are the coefficients in the kinetic equation for the first retraction (reloading) path, and the internal time t± is governed by the differential equation dt ¼ k: dt

ð47Þ

Eqs. (3), (46) and (47) are similar to appropriate relations in the endochronic theory of plasticity [36,45]. The difference between the present approach and the conventional concept is that the functions k± are piece-wise constant: the values of k± are altered only when the strain rate changes its sign. To account for the effect of cyclic deformation on the coefficients k±, we assume these quantities to decrease linearly with the increment of intensity of elastic strain DJe, k  ¼ 1  K  DJ e ;

ð48Þ

where DJe is given by Eq. (37), and Je0 is the intensity of elastic strain when the first retraction (reloading) starts. Eq. (48) implies that the parameters k± obey the zero-order kinetic equations dk  ¼ K  ; dDJ e

log K  ¼ K 0  K 1 log W p ;

4. Material parameters To describe an algorithm for the determination of material constants, we begin with simplification of the stress–strain relation for uniaxial tension of a specimen. At uniaxial deformation of an incompressible medium, the strain tensor reads   1 ^ ¼ ðtÞ e1  e1  ðe2  e2 þ e3  e3 Þ ; ð51Þ 2 where (t) stands for longitudinal engineering strain, ek (k = 1, 2, 3) are unit vectors of a Cartesian coordinate frame, whose vector e1 coincides with the direction of deformation, and  denotes tensor product. The plastic strain tensor ^p is presented in the form (51),   1 ^p ¼ p ðtÞ e1  e1  ðe2  e2 þ e3  e3 Þ ; ð52Þ 2 where p(t) is a function to be found. Insertion of Eqs. (51) and (52) into Eq. (3) results in dp d ¼/ ; dt dt

p ð0Þ ¼ 0:

ð53Þ

Substituting expressions (51) and (52) into Eq. (32) and excluding the unknown pressure P from the boundary condition on the lateral surface of the specimen, we find the engineering tensile stress, rðtÞ ¼ Eð1  /ðtÞÞ ððtÞ  p ðtÞÞ  Z t Z 1 pðvÞCðvÞdv expðCðvÞðt  sÞÞððsÞ  p ðsÞÞds ; j 0

0

ð54Þ

ð49Þ

where the increment of intensity of elastic strain DJe plays the role of internal time. The rates of evolution of the coefficients k± are assumed to obey the phenomenological relations similar to Eqs. (40) and (44), ð50Þ

where K0± and K1± are material constants. Constitutive Eqs. (3), (32), (33), (36) and (43) together with phenomenological relaxation (7), (10), (40), (44), (45), (48) and (50) involve 21 adjustable parameters l; V ; R; c; j; a; b; A0 ; A1 ; a00þ ; a01þ ; b00þ ; b01þ ; a00 ; a01 ; b00 ; b01 ; K 00þ ; K 01þ ; K 00 ; K 01 : Although this number is not small, it is comparable with the number of material constants in other models for the viscoelastic and viscoplastic responses of semicrystalline polymers. An advantage of the constitutive equations is that the adjustable parameters can be found one after another by fitting the experimental data reported in Section 2 in such a way that no more than three parameters are determined in the approximation of each interval of a stress–strain curve.

739

where E ¼ modulus.

3 l 2

stands for an analog of the Young’s

4.1. Relaxation tests First, we consider the standard relaxation test with an engineering strain 0. Insertion of the expression  0; t 6 0; ðtÞ ¼ ð55Þ 0 ; t > 0 into Eq. (53) yields ( 0; t 6 0; p ðtÞ ¼ 0p ; t > 0;

ð56Þ

where 0p stands for a plastic strain reached at the instant when the relaxation process starts. Substitution of expressions (55) and (56) into Eq. (54) implies that  Z rðtÞ ¼ Eð1  /0 Þ 1  j

Z

1

t

pðvÞCðvÞdv

0

0

 expðCðvÞðt  sÞÞds ð0  0p Þ;

0

where / is the value of / at the beginning of the relaxation process. Introducing the notation r0 ¼ Eð1  /0 Þð0  0p Þ;

740

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

using Eq. (1) and calculating the integral over s, we arrive at the formula Z 1 ðtÞ ¼ 1  j pðvÞ½1  expðCðvÞtÞ dv: ð57Þ r 0

Eq. (57) implies that the dimensionless relaxation curves  plotted versus time (the dimensionless engineering stress r t) are independent of strains 0 at which relaxation tests are performed. This conclusion is confirmed by the experimental data depicted in Figs. 2 and 10 for HDPE-A and HDPE-B, respectively. To determine the material constants V, R, c and j in Eqs. (10) and (57), the following algorithm is applied. We begin with the analysis of observations in tensile relaxation tests on HDPE-A and fix some intervals [0, V0], [0, R0] and [0, c0], where the ‘‘best-fit’’ parameters V, R and c are assumed to be located. Each of these intervals in divided into J subintervals by the points V(i) = iDV, R(j) = jDR and c(l) = lDc with DV = V0/J, DR = R0/J and Dc = c0/J (i, j, l = 0, . . . , J  1). For each pair {V(i), R(j)}, the pre-factor p0 in Eq. (10) is determined from Eq. (9). The integral in Eq. (9) is evaluated numerically by Simpson’s method with the step Dv = 0.1 and M = 200 steps. For each triplet {V(i), R(j), c(l)}, the coefficient j is found by the least-squares method from the condition of minimum of the function F¼

X

2

num ðtm ÞÞ ; ðð rexp ðtm Þ  r

m

where the sum is calculated over all instants tm at which the exp is the dimensionless experimental data are reported, r num  stress measured in a test, and r is given by Eq. (57). e and ~c are determined from the The best-fit values Ve , R

condition of minimum of the function F on the set {V(i), R(j), c(l)}. When the best-fit parameters are found, the initial intervals [0, V0], [0, R0] and [0, c0] are replaced with the new e  DR; R e þ DR and ½~c intervals ½ Ve  DV ; Ve þ DV , ½ R Dc; ~c þ Dc, and the calculations are repeated. To ensure good agreement between the observations and the results of numerical simulation, this procedure is repeated three times with J = 10. After matching the observations on HDPE-A, we proceed with fitting the experimental data in relaxation tests on HDPE-B. A variant of the above algorithm with only three adjustable parameters is employed (to reduce the number of material constants to be compared for polyethylenes with different molecular weights, the same value of c is used for HDPE-A and HDPE-B). The best-fit material parameters are collected in Table 1. 4.2. Tension with a constant strain rate We consider uniaxial tension of a virgin specimen with a constant strain rate _ . It follows from Eqs. (33) and (51) that the function / reads b

/ ¼ 1  exp½að  p Þ :

ð58Þ

To find the material constants E, a and b, we match the experimental stress–strain curves obtained at loading with a cross-head speed of 10 mm/min ð_ ¼ 0:002 s1 Þ and reported in Figs. 1 and 9. These observations are approximated with the help of the following algorithm. First, some intervals [0, a0] and [0, b0] are fixed, where the bestfit parameters a and b are assumed to be located. These intervals are divided by the points a(i) = iDa and b(j) = jDb Table 1 Adjustable parameters in the stress–strain relations Parameter

 versus time t in tensile relaxation tests Fig. 10. The dimensionless stress r at various strains . Symbols: experimental data on HDPE-B. Solid line: results of numerical simulation.

HDPE-A

HDPE-B

E (GPa) a b A0 A1 V R j c (s1)

1.38 9.37 0.79 0.43 0.17 2.36 4.28 0.67 1.50

2.86 3.98 0.50 0.32 0.34 4.69 2.60 0.62 1.50

a00þ

2.50

3.76

a01þ

0.50

0.41

b00þ

2.17

1.04

b01þ

0.81

1.00

a00

2.87

3.73

a01

0.21

0.64

b00

1.34

0.65

b01

0.21

0.54

K0+ K1+ K0 K1

0.13 1.15 0.82 0.22

0.98 0.80 1.04 0.45

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

with Da = a0/J and Db = b0/J (i, j = 1, . . . , J  1). For each pair {a(i), b(j)}, stress–strain relations (53) and (54) are integrated numerically together with phenomenological equations (7), (10) and (58) (by the Runge–Kutta method with the time step Dt = 0.001 s) from  = 0 to  = 0.1. The integrals over v are evaluated by using the Simpson method with the step Dv = 0.1, M = 200 steps, and with the material constants V, R and j listed in Table 1. The pre-factor p0 is found from the normalization condition (9). The elastic modulus E is determined by the least-squares method from the condition of minimum of the function X 2 ðrexp ðm Þ  rnum ðm ÞÞ ; ð59Þ F¼ m

where the sum is calculated over all points m at which observations are reported, rexp is the stress measured in a test, and rnum is given by Eq. (54). The best-fit parameters ~ are chosen from the condition of minimum of the ~ a and b function F on the set of pairs {a(i), b(j)}. After finding these quantities, the initial intervals [0, a0] and [0, b0] are replaced with the new intervals ~  Db; b ~ þ Db, and the above calcu½~ a  Da; ~ a þ Da and ½b lations are repeated. To ensure an acceptable quality of fitting observations, this procedure is repeated three times with J = 10. The best-fit material constants E, a and b are found separately for HDPE-A and HDPE-B. These quantities are listed in Table 1. To show that the constitutive equations can predict the mechanical response of HDPE in a tensile test with one strain rate when the material parameters are found by fitting data in a tensile test with another strain rate, we perform numerical simulation of uniaxial tension with a cross-head speed of 50 mm/min ð_ ¼ 0:01 s1 Þ. The results of numerical analysis are depicted in Figs. 1 and 9 together with appropriate experimental data. These figures demonstrate good agreement between the observations and predictions of the model.

741 b

/ðtÞ ¼ /0 þ ½1  expðAðe ðtÞ  0e Þ Þ

ðt P t0 Þ:

ð61Þ

To perform numerical integration of the stress–strain relations with r = r0, Eq. (54) is transformed. Setting Z t expðCðvÞðt  sÞÞe ðsÞ ds ð62Þ Zðt; vÞ ¼ CðvÞ 0

with e ðtÞ ¼ ðtÞ  p ðtÞ

ð63Þ

and resolving Eq. (54) with respect to e(t), we find that Z 1 r0 þj e ðtÞ ¼ pðvÞZðt; vÞ dv: ð64Þ Eð1  /ðtÞÞ 0 Differentiation of Eq. (63) with respect to time implies that oZ ðt; vÞ ¼ CðvÞ½e ðtÞ  Zðt; vÞ: ot

ð65Þ

The initial condition for Eq. (65) reads Zðt0 ; vÞ ¼ Z 0 ðvÞ;

ð66Þ

0

where Z (v) is determined by integration of the constitutive equations along the loading path of the stress–strain diagram. To establish a connection between the rates of elastic and plastic deformations, we differentiate Eq. (63) with respect to time, use Eq. (53) and obtain dp / de ¼ : 1  / dt dt

ð67Þ

To determine the best-fit value of A in Eq. (61), we introduce some interval [0, A0], where the adjustable parameter A is located, and divide this interval by the points A(i) = iDA with DA = A0/J (i = 1, . . . , J  1). For each A(i), the governing equations are integrated numerically

4.3. Creep tests To find the adjustable parameter A, the experimental data are approximated in tensile creep tests with various stresses r0 (Figs. 3 and 11). The deformation program for a standard creep test reads ðtÞ ¼ _ t

ð0 6 t < t0 Þ; 0

rðtÞ ¼ r0

ðt P t0 Þ;

ð60Þ 0

where the instant t is found from the condition r(t ) = r0. Each creep curve is fitted separately with the help of the following algorithm. First, we integrate governing Eqs. (53), (54) and (58) along the loading path of the stress– strain curve (from  = 0 to r = r0) with the material parameters reported in Table 1 and find the elastic strain 0e , the plastic strain 0p and the value /0 at the beginning of the creep process. Afterwards, numerical integration is performed of the governing Eqs. (53) and (54) together with phenomenological relations (35) and (36), which imply that in a tensile creep test,

Fig. 11. The engineering strain  versus time t in tensile creep tests with various engineering stresses r (MPa). Symbols: experimental data on HDPE-B. Solid lines: results of numerical simulation.

742

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

(by the Runge–Kutta method with the time step Dt = 0.001 s). The integrals over v are calculated by the Simpson method with the step Dv = 0.1 and M = 200 steps. Integration of kinetic equation (65) with initial condition (66) is performed for any v = mDv (m = 0, . . . , M  1). Given an instant t, the elastic strain e(t + Dt) is determined by Eq. (63). The strain (t + Dt) at macro-deformation is found from Eq. (53), ðt þ DtÞ  ðtÞ ¼

1 ½e ðt þ DtÞ  e ðtÞ; 1  /ðtÞ

ð68Þ

where the function /(t) is given by Eq. (61). The best-fit e is chosen from the condition of minimum parameter A (on the set {A(i)}) of the function X 2 F¼ ðexp ðtm Þ  num ðtm ÞÞ ; m

where the sum is calculated over all instants tm at which the experimental data are reported, exp is the strain measured in the test, and num is given by Eq. (68). e is found, the initial interval [0, A0] is replaced When A e  DA; A e þ DA, and the calculawith the new interval ½ A tions are repeated. This procedure is repeated three times with J = 10. Given a stress r0, we calculate the work Wp of tensile force on plastic deformations by means of Eq. (41), Z t0 dp ðtÞ dt; ð69Þ rðtÞ Wp ¼ dt 0 where the instant t0 corresponds to the beginning of an appropriate creep test, and the derivative dp/dt is given by Eq. (67). The adjustable parameter A is plotted versus Wp in Fig. 17. The experimental data are approximated by Eq. (40), where the coefficients A0 and A1 are found by the least-squares technique. These quantities are listed in Table 1.

/ðt0 Þ ¼ /0 : The function Z(t, v) is governed by kinetic equations (63) and (65) with initial condition (66). The quantities /0 and Z0(v) are determined by integration of the constitutive equations along the loading path of the stress–strain curve (from  = 0 to  = max). To find the coefficients a0þ and b0þ in Eq. (72), we match the first retraction path of each stress–strain curve depicted in Figs. 4 and 12 separately. For this purpose, we fix some intervals [0, a0] and [0, b0], where the best-fit values of a0þ and b0þ are assumed to be located, and divide these intervals by the points a(i) = iDa and b(j) = jDb with Da = a0/J, Db = b0/J (i, j = 1, . . . , J  1). For each pair {a(i), b(j)}, Eqs. (53), (71) and (72) are integrated numerically (by the Runge–Kutta method with the step Dt = 0.001 s) from  = max to r = 0. The best-fit parameters ~a and ~b are determined from the condition of minimum of function (59) on the set of pairs {a(i), b(j)}. When these values are found, the initial intervals [0, a0] and [0, b0] are replaced with the new intervals ½~a  Da; ~a þ Da and ½~b  Db; ~b þ Db, and the calculations are repeated. This procedure is repeated three times with J = 10. When the best-fit parameters a0þ and b0þ are found for each max, we calculate the work Wp of tensile force at plastic deformations along the loading path of each stress– strain curve by means of Eq. (69) and plot the quantities a0þ and b0þ versus Wp in Fig. 18. The experimental data are approximated by phenomenological relations (44), where the coefficients a0iþ and b0iþ (i = 0, 1) are determined by the least-squares method. These parameters are collected in Table 1.

4.4. First retraction For the first cycle of loading–retraction with a strain rate _ and a maximum strain max, the deformation program reads ðtÞ ¼ _ t

ð0 6 t < t0 Þ;

ðtÞ ¼ max  _ t

ðt P t0 Þ;

ð70Þ

0

where the instant t is found from the condition (t0) = max. It follows from Eqs. (54) and (62) that the engineering stress is given by   Z 1 rðtÞ ¼ Eð1  /ðtÞÞ ððtÞ  p ðtÞÞ  j pðvÞZðt; vÞ dv : 0

ð71Þ At the first retraction, the function /(t) obeys evolution equation (43) with a ¼ a0þ and b ¼ b0þ , d/ ¼ a0þ ð1  /ðtÞÞððtÞ  p ðtÞÞ_ þ b0þ / dt and the initial condition

ð72Þ

Fig. 12. The engineering stress r versus engineering strain  in cyclic tensile tests with various maximum strains max and the minimum stress rmin = 0 MPa. Symbols: experimental data on HDPE-B. Solid lines: results of numerical simulation.

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

4.5. Second cycle of deformation To determine the quantities a ¼ a0 and b ¼ b0 , we approximate the experimental data at the second cycle of deformations. For the first reloading path of each stress–strain diagram, Eq. (43) reads d/ ¼ a0 ð1  /Þð  p Þ_  b0 /: dt

743

algorithm with two adjustable parameters k and k+, where k and k+ are coefficients in Eqs. (43)–(45). According to these equations, the evolution equation for the function / reads d/ ¼ k  ½a0 ð1  /Þð  p Þ_  b0 / dt

ð75Þ

ð73Þ

For the second retraction path, the evolution equation is presented in the form (43) and (45), d/ ¼ k þ ½a0þ ð1  /Þð  p Þ_ þ b0þ /; dt

ð74Þ

where a0þ and b0þ are found by matching the first retraction path. Each stress–strain curve is fitted separately with the help of the following algorithm. We fix some intervals [0, a0] and [0, b0], where the best-fit parameters a0 and b0 in Eq. (73) are located, and an appropriate interval [0, k0] for the parameter k+ in Eq. (74). These intervals are divided by the points a(i) = iDa, b(j) = jDb and k(l) = lDk with Da = a0/J, Db = b0/J and Dk = k0/J (i, j, l = 1, . . . , J  1). For any triplet {a(i), b(j), k(l)}, the governing equations are integrated numerically (integration over time is performed by the Runge–Kutta method with the step Dt = 0.001 s). First, Eqs. (53), (65), (71) and (73) are integrated from r = rmin to  = max with the coefficients a(i) and b(j). Afterwards, Eqs. (53), (65), (71) and (74) are integrated from  = max to r = rmin with the coefficient k(l). The best-fit parameters ~ a, ~ b and ~k are determined from the condition of minimum of functional (59) on the set {a(i), b(j), k(l)}. Then, the initial intervals [0, a0], [0, b0] and [0, k0] are replaced with the new intervals ½~ a  Da; ~a þ Da, ½~ b  Db; ~ b þ Db and ½~k  Dk; ~k þ Dk, and the calculations are repeated. To guarantee an acceptable accuracy of fitting, this procedure is repeated three times with J = 10. When the adjustable parameters a0 and b0 are found by fitting the second cycle of the stress–strain diagrams with each maximum strain max, the work Wp of tensile force on plastic deformations is calculated from Eq. (69), where t0 is the instant when the first reloading starts. The quantities a0 and b0 are plotted versus Wp in Fig. 19. The experimental data are approximated by phenomenological relations (44), where the coefficients a0i and b0i (i = 0, 1) are found by the least-squares technique. The material constants a0i and b0i are collected in Table 1.

Fig. 13. The engineering stress r versus engineering strain  in a cyclic tensile test with the maximum strain max = 0.03 and the minimum stress rmin = 0 MPa. Circles: experimental data on HDPE-B. Solid line: results of numerical simulation.

4.6. Other cycles of deformation We proceed with fitting the experimental stress–strain diagrams depicted in Figs. 5–8 for HDPE-A and Figs. 13– 16 for HDPE-B. Each stress–strain curve is approximated separately. A loading–retraction path of a stress–strain diagram corresponding to the nth cycle of deformation (n = 3, . . . , 10) is matched by using a version of the above

Fig. 14. The engineering stress r versus engineering strain  in a cyclic tensile test with the maximum strain max = 0.05 and the minimum stress rmin = 0 MPa. Circles: experimental data on HDPE-B. Solid line: results of numerical simulation.

744

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

be located, and divide these intervals by the points k(i) = iDk and k(j) = jDk with Dk = k0/J (i, j = 1, . . . , J  1). For each pair {k(i), k(j)}, Eqs. (53), (65), (71) and (75) are integrated from r = rmin to  = max with the coefficient k(i), and, afterwards, Eqs. (53), (65), (71) and (76) are integrated from  = max to r = rmin with k(j). Numerical simulation is performed by the Runge–Kutta method with the time step Dt = 0.001 s. The best-fit parameters ~k  and ~k þ are determined from the condition of minimum of functional (59) on the set of pairs {k(i), k(j)}. When these quantities are found, the initial intervals are replaced with ½~k   Dk; ~k  þ Dk and ½~k þ  Dk; ~k þ þ Dk, respectively, and the calculations are repeated. This procedure is repeated three times with J = 10. The best-fit values of k+ and k are plotted versus increment of elastic strain De [the difference between the current value of e and its value 0e for the first retraction (reloading)] in Fig. 20 for HDPE-A and in Fig. 21 for HDPE-B. The experimental data are approximated by Eq. (48), Fig. 15. The engineering stress r versus engineering strain  in a cyclic tensile test with the maximum strain max = 0.07 and the minimum stress rmin = 0 MPa. Circles: experimental data on HDPE-B. Solid line: results of numerical simulation.

k þ ¼ 1  K þ De ;

k  ¼ 1  K  De ;

ð77Þ

where the coefficients K± are determined by the lastsquares algorithm. Given a maximum strain in a cyclic test max, we calculate the work Wp on tensile force on plastic deformations along the first loading path (for K+) and the first loading–retraction path of the stress–strain diagram (for K) and plot the coefficients K± versus Wp in Fig. 22. The experimental data are approximated by Eq. (50), where Ki± (i = 0, 1) are found by the least-squares method. These coefficients are collected in Table 1. 4.7. Accuracy of fitting

Fig. 16. The engineering stress r versus engineering strain  in a cyclic tensile test with the maximum strain max = 0.09 and the minimum stress rmin = 0 MPa. Circles: experimental data on HDPE-B. Solid line: results of numerical simulation.

for reloading and d/ ¼ k þ ½a0þ ð1  /Þð  p Þ_ þ b0þ / dt

ð76Þ

for retraction. To find the best-fit values of k and k+, we fix some intervals [0, k0], where these parameters are assumed to

Our aim now is to assess the accuracy of approximation of the experimental data by the constitutive equations. Figs. 2 and 10 demonstrate that Eq. (57) with three adjustable parameters V, R and j correctly describes stress relaxation in HDPE-A and HDPE-B. Small deviations of the relaxation curves in tests with various strains 0 from each other and from the results of numerical simulation may be ascribed to the effect of loading paths on the decay in stresses. Figs. 1 and 9 show that the constitutive equations with three adjustable parameters E, a and b adequately describe the stress–strain curves at tension with a cross-head speed of 10 mm/min. The model with the same material parameters correctly predicts the stress–strain curves in tensile tests with a cross-head speed of 50 mm/min, which implies that no additional dependencies are needed to account for the effect of strain rate of these parameters. Figs. 3 and 11 demonstrate that the stress–strain relations with the only adjustable parameter A correctly describe the response of HDPE in tensile creep tests with various stresses r0. According to Fig. 17, the pre-factor A monotonically changes with Wp and its dependence on Wp is adequately approximated by Eq. (40).

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

745

Fig. 17. The adjustable parameter A versus work Wp of tensile force on plastic deformations. Symbols: treatment of observations in creep tests on HDPE-A and HDPE-B. Solid lines: approximation of the experimental data by Eq. (40).

Fig. 19. The adjustable parameters a0 and b0 versus work Wp of tensile force on plastic deformations. Symbols: treatment of observations in cyclic tensile tests on HDPE-A and HDPE-B. Solid lines: approximation of the experimental data by Eq. (44).

Figs. 4 and 12 reveal that the constitutive equations with four adjustable parameters (two of them are found by fitting the retraction path and two others are determined by matching the reloading path of each stress–strain curve) correctly describe observations on HDPE-A and HDPEB in cyclic tensile tests with all maximum strains max under

consideration. Figs. 18 and 19 show that the dependencies of the quantities a0 and b0 on Wp are adequately described by phenomenological equation (44). Figs. 5–8 and 13–16 demonstrate excellent agreement between the experimental data in cyclic tensile tests on HDPE-A and HDPE-B and the results of numerical simulation. The stress–strain diagram for each cycle of deformation is determined by only two adjustable parameters, k+ and k. According to Figs. 20 and 21, the dependencies of these parameters on De are correctly approximated by phenomenological relations (77). Weak deviations between the experimental data and the results of numerical analysis are observed at the maximum strain max = 0.1 for HDPE-A (they may be ascribed by the fact that this maximum strain is located in the close vicinity of the yield strain for this grade of HDPE) and at the maximum strain max = 0.03 for HDPE-B (these inconsistencies may be explained by the fact that the rate of strain did not remain constant along reloading and retraction paths of a stress–strain diagram with such a small maximum strain). Fig. 22 shows that phenomenological Eq. (50) adequately describes the effect of work Wp of external forces on plastic deformations on the adjustable parameters K+ and K. This is an important advantage of the model, as the coefficients K± are obtained by means of a rather complicated procedure of fitting: (i) first, we approximate the experimental stress–strain curves to calculate k+ and k, (ii) afterwards, the dependencies of these quantities on De are approximated by relations (77), and (iii), finally, the coefficients K± are fitted with the help of Eq. (50).

Fig. 18. The adjustable parameters a0þ and b0þ versus work Wp of tensile force on plastic deformations. Symbols: treatment of observations in cyclic tensile tests on HDPE-A and HDPE-B. Solid lines: approximation of the experimental data by Eq. (44).

746

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

Fig. 20. The adjustable parameters k+ and k versus increment of elastic strain De. Symbols: treatment of observations on HDPE-A in cyclic tensile tests with various maximum strains max. Solid lines: approximation of the experimental data by Eq. (77).

Fig. 22. The adjustable parameters K+ and K versus work Wp of tensile force on plastic deformations. Symbols: treatment of observations in cyclic tensile tests on HDPE-A and HDPE-B. Solid lines: approximation of the experimental data by Eq. (50).

The following conclusions are drawn from Table 1:

Fig. 21. The adjustable parameters k+ and k versus increment of elastic strain De. Symbols: treatment of observations on HDPE-B in cyclic tensile tests with various maximum strains max. Solid lines: approximation of the experimental data by Eq. (77).

5. Discussion Our aim now is (i) to discuss the influence of molecular weight of HDPE on the material constants listed in Table 1, and (ii) to speculate regarding the physics at the microlevel that may explain these effects.

1. The elastic modulus E of HDPE-B noticeably (more than by twice) exceeds that of HDPE-A. This may be explained by two reasons: (i) the growth of molecular weight results in an increase in the content of entanglements in the amorphous phase, and as a consequence, an increase in the concentration of strands per unit volume, and (ii) the growth of molecular weight leads to an increase in the number of segments between entanglements and, as a consequence, to a more regular organization of segments in crystalline lamellae, which, in turn, induces growth of their rigidity. 2. It appears that both explanations are important, because the values of j are rather close for HDPEA and HDPE-B. In accord with Eq. (31), the latter means that lc grows with molecular weight being (roughly) proportional to N. 3. The parameter a of HDPE-A strongly exceeds that of HDPE-B. According to Eqs. (53) and (58), this implies that plastic flow in HDPE-A starts at substantially lower longitudinal strains than in HDPEB. Bearing in mind that the beginning of plastic flow is conventionally associated with fine slip (homogeneous shearing) of lamellar stacks, the fact that this flow begins earlier in HDPE-A with a low molecular weight means that less entangled chains in amorphous regions between lamellae are weaker than appropriate strongly entangled chains in HDPE-B with a high molecular weight. 4. The exponent b in Eq. (58) is higher for HDPE-A than for HDPE-B. As b serves as a measure of non-

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

5.

6.

7.

8.

linearity of the mechanical response, this means that the nonlinearity of the stress–strain diagram of HDPE-B at tension is more pronounced than that of HDPE-A (in agreement with the observations depicted in Figs. 1 and 9). The smaller value of b for HDPE-B may be ascribed to the fact that b is responsible for the effect of elastic strain on the rate of developed viscoplastic flow. This developed flow (associated with coarse slip of lamellar blocks and necking of individual lamellae) occurs more pronouncedly in HDPE-B with a high molecular weight and a better organization of crystallites (higher rigidity lc). The average activation energy V for detachment of strands from the network in HDPE-B exceeds (about by twice) that in HDPE-A. As separation of strands from their junctions is associated with disentanglements of chains in amorphous regions and detachment of tie chains and chain loops from the lamellar surfaces, this conclusion seems to be physically plausible: (i) an increase in molecular weight makes more difficult disentanglement of longer chains in a network, whereas an improvement of quality of crystallites (their more regular organization) prevents detachment of tie chains. The standard deviation of activation energies R is larger for HDPE-A than for HDPE-B. Keeping in mind that R characterizes a degree of inhomogeneity of a transient network, this conclusion appears to be in accord with our previous findings: more regular organization of crystalline lamellae in HDPE-B with a higher molecular weight results in a less pronounced inhomogeneity of the amorphous phase. Fig. 17 demonstrates that given a work Wp of external forces on plastic deformations, the parameter A in Eq. (36) is higher for HDPE-A than for HDPE-B. This is in agreement with the above conclusion for the pre-factor a in Eq. (33), and may be associated with a more pronounced fine slip of lamellar blocks in polyethylene with a lower molecular weight. The only difference between these results is that A reflects homogeneous shearing of lamellar stacks observed in creep tests, whereas a is responsible for the same process at uniaxial tension. Fig. 18 shows practically the same effect of the work Wp of external forces at plastic deformations on the parameters a0þ and b0þ of HDPE-A and HDPE-B. These quantities monotonically decrease with Wp and the slopes of the curves a0þ ðW p Þ and b0þ ðW p Þ (plotted in the double-logarithmic coordinates) are similar for polyethylenes with high and low molecular weights. At the same time, given a Wp, the quantities a0þ and b0þ for HDPE-B exceed those for HDPE-A by an order of magnitude. It follows from Eq. (43) that an increase in a0þ and b0þ is tantamount to the growth of rate of viscoplastic flow at the first retraction. In other words, HDPE with a more regular organization

747

of crystallites returns to its initial state at unloading easier than polyethylene with a non-regular structure of crystallites, which appears to be physically plausible. 9. According to Fig. 19, the adjustable parameters a0 and b0 decay with work Wp of external forces at plastic deformations. In other words, the more severely the crystalline structure was altered at active loading, the less substantial changes occur in the rate of viscoplastic deformations at the first reloading. The rates of decrease in a0 and b0 with Wp practically coincide for each grade of polyethylene, but these rates are higher for HDPE-B than for HDPE-A. This means that HDPE-B with a higher molecular weight (more regular structure of crystallites) suffers stronger changes in its crystalline morphology that HDPE-A with a lower molecular weight (less regular structure of crystallites). Given a work Wp, the quantities a0 and b0 found by matching observations on HDPE-B are higher than those on HDPE-A, that is they follow the same trend as the parameters a0þ and b0þ in Fig. 18. 10. Fig. 22 demonstrates that K+ and K decrease with work Wp of external forces at plastic deformations. The decay in K+ is more pronounced for HDPE-B than for HDPE-A, whereas an inverse dependence is observed for K. Given a work of tensile force Wp, the parameters K+ and K for HDPE-B exceed those for HDPE-A. As the coefficients K+ and K are responsible for slowing down of the rates of viscoplastic deformation with number of cycles, we conclude that ‘‘tuning’’ of the crystalline morphology (driven by cyclic deformations) in HDPE-B with a higher molecular weight is stronger than that in HDPE-A with a lower molecular weight. 6. Numerical simulation To reveal some characteristic features of the time- and rate-dependent behavior of HDPE, (i) numerical simulation is performed of the constitutive equations under more sophisticated deformation programs than those reported in Section 2, and (ii) the stress–strain diagrams are compared with available experimental data on semicrystalline polymers. All calculations are performed with the material constants for HDPE-A (Table 1). We begin with the analysis of the mechanical response in uniaxial tensile tests. Following [30], two deformation programs are considered: (i) tension with a constant strain rate _ ¼ 0:01 s1 , and (ii) tension with a piece-wise strain rate 8 0:01 s1 ; 0 6  < 0:03; > > > > > > < 0:001 s1 ; 0:03 6  < 0:06; ð78Þ _ ¼ > 1 > ; 0:06 6  < 0:09; 0:1 s > > > > : 0:01 s1 ; 0:09 6  < 0:1:

748

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

The results of numerical simulation are plotted in Fig. 23. The stress–strain curves are in qualitative agreement with the observations depicted in Fig. 3 of [16] and Fig. 3 of [18] at uniaxial tension and in Fig. 11 of [30] at uniaxial compression of HDPE with large strains. Quantitative comparison of our results with experimental data is rather difficult, as conventional equipment does not allow changes in the strain rate (by an order of magnitude) to be performed at uniaxial tension within small intervals of strains with the length of 0.03. We proceed with the study of relaxation curves in uniaxial tensile tests with the strain 0 = 0.1, when the stress relaxation followed active loading immediately, as well as when relaxation is performed after several cycles of loading–retraction down to the zero stress. Although observations in such tests on HDPE are absent in the literature, experimental data on isotactic polypropylene were reported in [59]. The results of numerical simulation are presented in  is plotted versus Fig. 24, where the dimensionless stress r the logarithm of time (the instant t = 0 corresponds to the beginning of the relaxation process). This figure demonstrates that given a time t, the amount of relaxing stress ðtÞ decreases with number of cycles. With a D rðtÞ ¼ 1  r high level of accuracy, the relaxation curves may be shifted along the time axis to construct a master-curve. Both conclusions are in accord with the observations presented in [59]. Our next purpose is to assess the effect of strain rate at active loading on creep curves in conventional tensile creep tests. The results of numerical simulation are depicted in Fig. 25, where the strain  is plotted versus the logarithm of time t (the instant t = 0 corresponds to the beginning

Fig. 23. The engineering stress r versus engineering strain  in tensile tests. Symbols: results of numerical simulation. Circles: deformation with the constant strain rate _ ¼ 0:01 s1 . Solid line: deformation with the piecewise constant strain rate (78).

 versus time t in a standard relaxation Fig. 24. The dimensionless stress r test at  = 0.1 (curve 1), in a relaxation test after one cycle of tension– retraction with the maximum strain max = 0.1 and the minimum stress rmax = 0 MPa (curve 2) and in a relaxation test after two cycles of tension– retraction with the maximum strain max = 0.1 and the minimum stress rmax = 0 MPa (curve 3). Solid lines: results of numerical simulation. The strain rate at tension and retraction is _ ¼ 0:002 s1 .

of creep). The creep curves start at the same tensile strain  = 0.05, which is reached with various strain rates _ . Fig. 25 shows that given an instant t, the creep compliance D = (t)  (0) is noticeably reduced with an increase in

Fig. 25. The engineering strain  versus time t in tensile creep tests that start when the strain  reaches 0.05. Solid lines: results of numerical simulation. Along the loading part of the stress–strain curves, deformation occurs with the strain rates _ ¼ 0:1 (curve 1), 0.01 (curve 2), 0.001 (curve 3) and 0.0001 s1 (curve 4).

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

strain rate. This conclusion is in qualitative agreement with the experimental data presented in Fig. 10 of [21]. It should be noted, however, that our creep curves differ substantially from those reported in [21] on poly(ethylene-co-vinyl acetate). The experimental dependencies (log t) depicted in Figs. 3 and 11 are strongly nonlinear, whereas those plotted in [21] may be approximated fairly well by straight lines. We proceed with the study of the time-dependent response of HDPE in tensile creep tests with  = 0.05 followed the active loading immediately, as well as when creep starts after a few cycles of loading–retraction with the maximum strain max = 0.05 and the minimum stress rmin = 0 MPa. The results of numerical analysis are depicted in Fig. 26. This figure demonstrates that given a creep time t, the creep compliance decreases with number of cycles. The decay in engineering strain (which may be associated with an apparent hardening of the polymer) is relatively weak when loading is performed with a high strain rate, and it becomes substantial at slow cyclic deformation. The creep curves (log t) can be superposed (with a high level of accuracy) by shifting along the time axis to construct a master-curve, but the master-curves obtained at different strain rates distinguish noticeably. Although appropriate tests on HDPE were not conducted, an apparent hardening of semicrystalline polymers in creep tests that followed several cycles of cyclic deformation was

Fig. 26. The engineering strain  versus time t in tensile creep tests that start when the strain  reaches 0.05. Symbols: results of numerical simulation in a standard creep test (curve 1), in a creep test that follows a cycle of loading–retraction with the maximum strain max = 0.05 and the minimum strain rmin = 0 MPa (curve 2), and in a creep test starting after two cycles of loading–retraction with the maximum strain max = 0.05 and the minimum strain rmin = 0 MPa (curve 3). The strain rates at tension and retraction are _ ¼ 0:0005 (unfilled circles) and _ ¼ 0:0015 s1 (filled circles).

749

observed in our experiments on isotactic polypropylene (to be published). To evaluate the mechanical response of HDPE in uniaxial tensile tests interrupted by intervals of stress relaxation (with a fixed duration of 3 min), simulation of the constitutive equations is performed for (i) tension with constant strain rates _ and (ii) tension with piece-wise strain rates of the form 8 1 > < 0:001 s ; 0 6  < 0:04; _ ¼ 0:01 s1 ; 0:04 6  < 0:08; ð79Þ > : 1 0:1 s ;  P 0:08 and

8 1 > < 0:0001 s ; 0 6  < 0:04; _ ¼ 0:001 s1 ; 0:04 6  < 0:08; > : 0:01 s1 ;  P 0:08:

ð80Þ

The results of numerical analysis are presented in Figs. 27 and 28. According to Fig. 27, an increase in strain rate results in more pronounced peaks on the stress–strain diagram. The amount of relaxing stress reaches its maximum at the first interval of relaxation, and it decreases with number of relaxation periods. Fig. 28 reveals that the growth of strain rate at subsequent intervals of tension causes an increase in the amount of relaxing stress during each interval of relaxation. The stress–strain curves depicted in Figs. 27 and 28 are in qualitative agreement with those reported in Fig. 2 of [25] and Figs. 8 and 9 of [60]. Finally, we analyze the mechanical response of HDPE in cyclic tensile tests. First, we consider tests where a specimen is deformed with a constant strain rate _ ¼ 0:002 s1 up to

Fig. 27. The engineering stress r versus engineering strain  in tensile tests interrupted by intervals of relaxation. Symbols: results of numerical simulation with the strain rates _ ¼ 0:1 (solid line), 0.01 (diamonds), 0.001 (unfilled circles) and 0.0001 s1 (filled circles).

750

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

Fig. 28. The engineering stress r versus engineering strain  in tensile tests interrupted by intervals of relaxation. Symbols: results of numerical simulation. Solid line: loading with the strain rate (79). Circles: loading with the strain rate (80).

the strain 1 = 0.03, unloaded down to some minimum stress rmin, reloaded up to the new strain 2 = 0.06, unloaded down to the same stress rmin as in the first cycle, reloaded up to the strain 3 = 0.09, unloaded down to the stress rmin and reloaded up to the strain 4 = 0.1. Similar cyclic programs were used in mechanical tests on high-density polyethylene in [13,18]. The results of numerical simu-

Fig. 30. The engineering stress r versus strain  in cyclic tensile tests with the strain rate _ ¼ 0:002 s1 , the maximum stress rmax = 18.0 MPa and various minimum stresses rmin (MPa). Symbols: results of numerical simulation. Solid line: rmin = 0.0. Unfilled circles: rmin = 2.0. Filled circles: rmin = 4.0.

lation are plotted in Fig. 29. The stress–strain diagrams in this figure are similar to those reported in Figs. 5 and 15 of [18]. Fig. 29 demonstrates that an increase in rmin results in a decrease in engineering stress at the intervals of tension, but does not affect it at retraction. The stress–strain diagrams in cyclic tensile tests with a stress-controlled program (the strain rate _ ¼ 0:002 s1 , the maximum stress rmax = 18 MPa and various minimum stresses rmin) are presented in Fig. 30. This figure shows that the growth of minimum stress rmin causes a monotonic increase in maximum strain per cycle. Although no observations on semicrystalline polymers are available that confirm this conclusion, it is in qualitative agreement with the experimental data on stainless steel in ratcheting tests [61]. 7. Concluding remarks

Fig. 29. The engineering stress r versus strain  in cyclic tensile tests with the strain rate _ ¼ 0:002 s1 , the maximum strains max = 0.03, 0.06 and 0.09 and various minimum stresses rmin (MPa). Symbols: results of numerical simulation. Solid line: rmin = 0.0. Unfilled circles: rmin = 2.0. Filled circles: rmin = 4.0.

Observations have been reported on two commercial grades of HDPE with different molecular weights in tensile tests with various strain rates below the yield point, standard relaxation tests at various strains, standard creep tests with various stresses and cyclic tests with a strain-controlled deformation program and various maximum strains. Based on the laws of thermodynamics, constitutive equations have been derived for the viscoelastic and viscoplastic responses of a semicrystalline polymer at arbitrary three-dimensional deformations with small strains. A semicrystalline polymer is treated as a two-phase continuum consisting of a crystalline skeleton surrounded by amorphous meso-regions. The crystalline phase is thought of

A.D. Drozdov, J. deC. Christiansen / Computational Materials Science 39 (2007) 729–751

as a viscoplastic medium, whose plastic flow is associated with fine and coarse slip of lamellar stacks. The amorphous phase is modeled as a viscoelasto-plastic medium. Its viscoelastic behavior is attributed to thermally induced rearrangement of strands in a transient network of chains, whereas its viscoplastic response reflects sliding of junctions with respect to their reference positions. The stress–strain relations involve 21 material constants that are found by fitting the experimental data. An advantage of the model is that the adjustable parameters are determined one after another in such a fashion that no more than three quantities are found by matching each interval of a stress–strain diagram. It is demonstrated that the constitutive equations correctly describe the observations, and material parameters are affected by molecular weight of HDPE in a physically plausible way. It is found that an increase in molecular weight leads to a more regular crystalline morphology. This regularity induces (i) growth of the elastic modulus and (ii) slowing down of viscoelastic processes (that are observed in standard creep and relaxation tests). At the same time, enhancement of regularity of the crystalline structure results in more pronounced changes in the rate of viscoplastic flow during the first cycles of retraction and reloading. These changes decay the stronger, the higher is the molecular weight of HDPE. The time- and rate-dependent behavior of a semicrystalline polymer in uniaxial tests with complicated deformation programs is studied by using numerical simulation. The results of numerical analysis are in qualitative agreement with available experimental data. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19]

R.A. Bubeck, Mater. Sci. Eng. Rep. 39 (2002) 1–28. A. Galeski, Prog. Polym. Sci. 28 (2003) 1643–1699. H.E.H. Meijer, L.E. Govaert, Prog. Polym. Sci. 30 (2005) 915–938. R. Seguela, J. Macromol. Sci. C 45 (2005) 263–287. R.A. Sousa, R.L. Reis, A.M. Cunha, M.J. Bevis, J. Appl. Polym. Sci. 89 (2003) 2079–2087. J. Zhang, K. Shen, Y. Gao, Y. Yuan, J. Appl. Polym. Sci. 96 (2005) 818–823. Q. Fu, Y. Men, G. Strobl, Polymer 44 (2003) 1927–1933. J. Kong, X.D. Fan, M. Jia, J. Appl. Polym. Sci. 93 (2004) 2542–2549. Z. Bartczak, Macromolecules 38 (2005) 7702–7713. Z. Bartczak, A. Kozanecki, Polymer 46 (2005) 8210–8221. X. Guan, R. Pitchumani, Polym. Eng. Sci. 44 (2004) 233–251. C. G’Sell, J.M. Hiver, A. Dahoun, A. Souahi, J. Mater. Sci. 27 (1992) 5031–5039. R. Hiss, S. Hobeika, C. Lynn, G. Strobl, Macromolecules 32 (1999) 4390–4403. S. Hobeika, Y. Men, G. Strobl, Macromolecules 33 (2000) 1827– 1833. S. Nikolov, I. Doghri, O. Pierard, L. Zealouk, A. Goldberg, J. Mech. Phys. Solids 50 (2002) 2275–2302. A. Dasari, R.D.K. Misra, Mater. Sci. Eng. A 358 (2003) 356–371. J. Lai, A. Bakker, Polymer 36 (1995) 93–99. C. Zhang, I.D. Moore, Polym. Eng. Sci. 37 (1997) 404–413. J.G.J. Beijer, J.L. Spoormaker, Polymer 41 (2000) 5443–5449.

751

[20] J.F. Mano, R.A. Sousa, R.L. Reis, A.M. Cunha, M.J. Bevis, Polymer 42 (2001) 6187–6198. [21] K. Hong, A. Rastogi, G. Strobl, Macromolecules 37 (2004) 10174– 10179. [22] R. Elleuch, W. Taktak, J. Mater. Eng. Perform. 15 (2006) 111–116. [23] D.S. Li, H. Garmestani, R.G. Alamo, S.R. Kalidindi, Polymer 44 (2003) 5355–5367. [24] K. Hong, A. Rastogi, G. Strobl, Macromolecules 37 (2004) 10165– 10173. [25] K. Mrabet, R. Rahouadj, C. Cunat, Polym. Eng. Sci. 45 (2005) 42–51. [26] B. Na, Q. Zhang, Q. Fu, Y. Men, K. Hong, G. Strobl, Macromolecules 39 (2006) 2584–2591. [27] N.-J. Jo, A. Takahara, T. Kajiyama, Polymer 38 (1997) 5195–5201. [28] F. Addiego, A. Dahoun, C. G’Sell, J.-M. Hiver, Polymer 47 (2006) 4387–4399. [29] B.A.G. Schrauwen, R.P.M. Janssen, L.E. Govaert, H.E.H. Meijer, Macromolecules 37 (2004) 6069–6078. [30] T. Kazmierczak, A. Galeski, A.S. Argon, Polymer 46 (2005) 8926– 8936. [31] Z. Bartczak, A.S. Argon, R.E. Cohen, Polymer 35 (1994) 3427–3441. [32] Z. Bartczak, E. Lezak, Polymer 46 (2005) 6050–6063. [33] C. Zhang, I.D. Moore, Polym. Eng. Sci. 37 (1997) 414–420. [34] J.A.W. van Dommelen, D.M. Parks, M.C. Boyce, W.A.M. Brekelmans, F.P.T. Baaijens, J. Mech. Phys. Solids 51 (2003) 519–541. [35] M.T. Suleiman, B.J. Coree, J. Mater. Civil Eng. 16 (2004) 511–515. [36] T. Kletschkowski, U. Schomburg, A. Bertram, Mech. Mater. 34 (2002) 795–808. [37] N. Brusselle-Dupend, D. Lai, X. Feaugas, M. Guigon, M. Clavel, Polym. Eng. Sci. 43 (2003) 501–518. [38] A.D. Drozdov, R.K. Gupta, Int. J. Eng. Sci. 41 (2003) 2335–2361. [39] A.D. Drozdov, Q. Yuan, Int. J. Solids Struct. 40 (2003) 2321–2342. [40] J.S. Bergstrom, C.M. Rimnac, S.M. Kurtz, Biomaterials 25 (2004) 2171–2178. [41] A.D. Drozdov, J. deC. Christiansen, Polym. Eng. Sci. 44 (2004) 548– 556. [42] V. Oshmyan, S. Patlazhan, Y. Remond, J. Mater. Sci. 39 (2004) 3577– 3586. [43] J.S. Bergstrom, L.B. Hilbert, Mech. Mater. 37 (2005) 899–913. [44] O.U. Colak, Int. J. Plasticity 21 (2005) 145–160. [45] T. Kletschkowski, U. Schomburg, A. Bertram, Mech. Mater. 37 (2005) 687–704. [46] Y. Remond, Compos. Sci. Technol. 65 (2005) 421–428. [47] Z. Xia, X. Shen, F. Ellyin, Polym. Eng. Sci. 45 (2005) 103–113. [48] F. Khan, E. Krempl, Trans. ASME J. Eng. Mater. Technol. 128 (2006) 64–72. [49] F. Tanaka, S.F. Edwards, Macromolecules 25 (1992) 1516–1523. [50] H. Eyring, J. Chem. Phys. 4 (1936) 283–291. [51] S. Fellahi, B.D. Favis, B. Fisa, J. Mater. Sci. 30 (1995) 5522–5530. [52] R.W. Ogden, D.G. Roxburgh, Proc. Roy. Soc. London A 455 (1999) 2861–2877. [53] T. Faravelli, G. Bozzano, C. Scassa, M. Perego, S. Fabini, E. Ranzi, M. Dente, J. Anal. Appl. Pyrolysis 52 (1999) 87–103. [54] K. Kitao, Polym. Eng. Sci. 41 (2001) 1146–1155. [55] R.K. Bayer, T.A. Ezquerra, H.G. Zachmann, F.J. Balta Calleja, J. Martinez Salazar, W. Meins, R.E. Diekow, P. Wiegel, J. Mater. Sci. 23 (1988) 475–480. [56] B. Derrida, Phys. Rev. Lett. 45 (1980) 79–92. [57] Y. Song, H. Yamamoto, N. Nemoto, Macromolecules 37 (2004) 6219–6226. [58] A.D. Drozdov, J. deC. Christiansen, Comput. Mater. Sci., doi:10.1016/j.commatsci.2006.07.014. [59] T. Ariyama, Polym. Eng. Sci. 33 (1993) 18–25. [60] S.G. Bardenhagen, M.G. Stout, G.T. Gray, Mech. Mater. 25 (1997) 235–253. [61] G. Kang, Q. Kan, J. Zhang, Y. Sun, Int. J. Plasticity 22 (2006) 858– 894.