Modeling the thermo-mechanical behavior and constrained recovery performance of cold-programmed amorphous shape-memory polymers

Journal Pre-proof Modeling the thermo-mechanical behavior and constrained recovery performance of cold-programmed amorphous shape-memory polymers Lu Dai, Chuanshuai Tian, Rui Xiao PII:

S0749-6419(19)30643-6

DOI:

https://doi.org/10.1016/j.ijplas.2019.102654

Reference:

INTPLA 102654

To appear in:

International Journal of Plasticity

Received Date: 5 September 2019 Revised Date:

23 December 2019

Accepted Date: 23 December 2019

Please cite this article as: Dai, L., Tian, C., Xiao, R., Modeling the thermo-mechanical behavior and constrained recovery performance of cold-programmed amorphous shape-memory polymers, International Journal of Plasticity (2020), doi: https://doi.org/10.1016/j.ijplas.2019.102654. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Authors’ contribution statement R.X. conceived the project. L.D. and C.T. performed the experiments. L.D. and R.X. performed the simulation. L.D. and R.X. analyzed the data. R.X. wrote the paper. All the authors commented on the paper.

Modeling the thermo-mechanical behavior and constrained recovery performance of cold-programmed amorphous shape-memory polymers Lu Daia , Chuanshuai Tianb , Rui Xiaoc,∗

a School

b Department

c Key

of Civil Engineering and Architecture, Zhejiang Sci-Tech University, Hangzhou, Zhejiang 310018, China

of Engineering Mechanics, College of Mechanics and Materials, Hohai University, Nanjing, Jiangsu 210098, China.

Laboratory of Soft Machines and Smart Devices of Zhejiang Province, Department of Engineering Mechanics, Zhejiang University, Hangzhou 310027, China

Abstract

The low recovery stress has limited utilizing shape-memory polymers (SMPs) as actuators. In this work, we demonstrate that the resultant recovery stress of amorphous SMPs can be significantly increased through cold programming. A three-dimensional model is formulated to describe the thermo-mechanical behavior and shape-memory performance of amorphous polymers in large deformation. The constitutive relationship is derived based on the two-temperature thermodynamic framework employing an effective temperature as a thermodynamic state variable to describe the nonequilibrium structure of amorphous polymers. The model also incorporates the molecular orientation with a relaxation mechanism to describe strain hardening. The model is applied to simulate the stress-strain relationship and stress relaxation behaviors of an amorphous SMP in the programming process and the stress response in the constrained recovery process. The model can accurately reproduce the measured stress response at different temperatures and loading rates. Both experimental and simulation results show that polymers deformed at a larger loading rate exhibit a lower stress in the holding period, though simulation pronouncedly underestimates the magnitude of stress drop in the stress relaxation process. The model also captures the magnitude and location of the peak stress during the constrained recovery process for SMPs programmed at different temperatures and with different applied strains. Keywords: shape-memory polymers; glassy polymers; strain hardening; effective temperature; viscoplasticity; stress relaxation

1. Introduction

The thermally activated shape-memory effect (SME) represents that the programmed temporary shape of polymers can return to the permanent shape when polymers are heated above transition temperatures (He et al. 2015, Eskandari et al. 2018, Lu et al. 2019, Li et al. 2019). Based on the recovery condition, the shape-memory performance can be evaluated in a free recovery test or a constrained recovery test. In the free recovery test, no force is applied while ∗ Corresponding

author. Email address: [email protected] (Rui Xiao )

Preprint submitted to Elsevier

December 31, 2019

the shape is measured with heating. In the constrained recovery test, the shape is fixed during the heating process. Thus, the shape recovery induces a resultant force (Belmonte et al. 2017, Abishera et al. 2018). The recovery force is vital for some medical devices based on shape-memory polymers (SMPs) (Fisher et al. 2015). Compared with various works on characterization the free recovery behaviors of SMPs, only a few works have been performed to investigate the constrained recovery performance of SMPs (Gall et al. 2005, Nguyen et al. 2008, Santiago et al. 2016, Fan and Li 2018). The reported recovery stress is typically several MPa for various polymer systems (Nguyen et al. 2008, Santiago et al. 2016, Abishera et al. 2018), which limits practical applications of SMPs as actuators. The constrained recovery performance also depends on programming temperatures showing an increase in the peak stress at a lower programming temperature (Xiao et al. 2013, Li and Wang 2016). Thus, it is necessary to perform a complete experimental characterization of constrained recovery performance of SMPs with different programming temperatures and strains. In addition to experimental characterization, constitutive modeling has also been widely used to investigate the shapememory performance of polymers (Ge et al. 2016, Mao et al. 2015, Lu et al. 2018, Li et al. 2018, Gu et al. 2019, Mao et al. 2019). As summarized in Nguyen (2013), two main approaches were widely adopted to model the thermally induced shape-memory behavior in amorphous polymers, namely the phase transition approach and the thermoviscoelastic approach. The phase transition approach (Qi et al. 2008) assumes that there coexist an active phase and a frozen phase. The transition between the two phases is the underlying mechanism for shape fixity and shape recovery. Though this approach is phenomenological, it can successfully reproduce many important features in shape-memory tests (Boatti et al. 2016). The thermoviscoelastic approach (Nguyen et al. 2008) can describe the SME through modeling the dramatic change in viscosity across the glass transition region. Yu et al. (2012) have shown that simple viscoelastic models can well capture the temperature memory effect and multiple memory effect in amorphous networks. Xiao et al. (2013) have introduced the effective temperature into the viscoelastic-plastic model to describe the effects of aging on the shape memory performance. However, this model can only be applied in the median strain region. Using the effective temperature to model the thermomechanical behaviors of amorphous solids will be introduced in detail in the late part of this section. Some recent works (Li et al. 2017, Su and Peng 2018) have combined the phase transition approach and the thermoviscoelastic approach to simulate the shape-memory behaviors in different polymer systems. Most amorphous polymers can be potentially used as shape-memory materials, though shape-memory performances still vary with the material systems. Thus, the widely existed models for amorphous polymers can be readily applied or modified to account for shape-memory behaviors (Yu et al. 2017, Jiang et al. 2017, Zhang and To 2016, Engqvist et al. 2016). Here we briefly review several constitutive models for glassy amorphous polymers developed in recent years (Kontou 2016, Johnsen et al. 2019). Krairi et al. (2019) formulated a thermo-mechanical coupled theory to simulate the viscous behaviors of thermoplastic polymers, which was also implemented for finite element analysis to simulate the evolution of temperature field during deformation. Jiang and Jiang (2019) introduced an internal state variable representing the molecular entanglement into an elasto-viscoplastic framework, which can successfully capture the dependence of the stress response of glassy polymers on thermo-mechanical history. Engqvist et al. (2016) employed a micro-structural deformation gradient to describe the orientation of microstructure. The mode captured inhomogeneous deformation behaviors in cold drawing tests of polycarbonate, and more importantly reproduced the fabric orientation measured by the X-ray scattering technique. Chevalier et al. (2018) showed that incorporation of heterogeneous activation of shear transformation zones at mesoscale can reproduce the main macroscale viscoplastic response. These physical based multi-scale models may provide insight on understanding the origin of the complex

2

behaviors in glassy polymers. To model the thermo-mechanical response of amorphous glassy polymers, a series of models have also been developed based on the concept of effective temperature (Semkiv and H¨utter 2016, Das et al. 2018). This thermodynamic framework assumes that the system of amorphous solids is consisted of a kinetic subsystem characterized by temperature and kinetic entropy and a configurational subsystem characterized by effective temperature and configurational entropy. The pioneering work of Langer and coauthors (Langer 2004, Falk and Langer 2011, Bouchbinder and Langer 2009) developed a nonequilibrium thermodynamic framework for amorphous solids and integrated the effective temperature with the theory of shear transformation zones (Falk and Langer 1998). Based on continuum mechanics, Kamrin and Bouchbinder (2014) formulated a finite deformation thermodynamic framework for the elasto-viscoplastic deformation based on effective temperature. Xiao and Nguyen (2015) incorporated multiple relaxation mechanisms into the effective temperature theory and firstly demonstrated the framework can capture the influences of temperature, loading rate, thermal treatment and deformation on the stress response. Semkiv et al. (2017) applied the two-temperature model to study aging in elastomers filled with nanoparticles by assuming there existing a glassy layer around the particles. Das and Roy (2019) showed the constitutive model based on effective temperature can successfully capture the stress response of polyurea at strain rates range from 0.01/s to 6500/s. Though the effective temperature theory has became a powerful tool to model the complex behaviors of amorphous solids, so far it has not been applied to describe the SME in amorphous polymers. Here we extend the model of Xiao and Nguyen (2015) through incorporation the strain hardening mechanism and for the first time demonstrate the effective temperature theory can simultaneously describe the thermo-mechanical response and shape-memory behaviors of amorphous solids under large deformation. The following part is arranged as below. Sec. 2 describes the experimental characterization of an amorphous thermoplastic poly(ethylene terephthalate)-glycol (PETG). Three sets of experiments are performed. The dynamic frequency sweep and differential scanning calorimetry tests are used to obtain the viscoelastic relaxation spectrum and structural relaxation spectrum respectively. Uniaxial compression tests are carried out to obtain the stress response. The constrained recovery tests are used to investigate the recovery stress in the heating process. Sec. 3 first presents a nonequilibrium thermodynamic framework for the theory and then gives an explicit form of the constitutive relations. The followed section shows the procedures of obtaining the model parameters. Sec. 5 evaluates the model performance by comparing simulation results with experimental data. The paper ends with summarizing the main findings of this work.

2. Experiments

The amorphous thermoplastic PETG was obtained from McMaster-Carr Supply Co. The materials are directly used without further chemical treatments.

2.1. Dynamic frequency sweep The DMA Q800 (TA Instruments, USA) was used to measure the dynamic properties of PETG. The size of the film specimen was 15 mm× 5 mm× 1.59 mm. The specimen was first kept at 100o C for 15 mins and then cooled to 30o C at 3o C/min. A dynamic frequency sweep scan was applied with a strain magnitude of 0.2% in the tension mode from 3

30o C to 120o C at a discrete manner of every 3o C. In order to avoid compressing the film specimen, a prestrain with a magnitude of 0.25% was applied in the testing process. The frequencies were chosen as 0.3 Hz, 1 Hz, 3 Hz and 10 Hz.

2.2. Differential scanning calorimetry A DSC 25 (TA Instruments, USA) was employed to perform the differential scanning calorimetry scanning (DSC). A film specimen with a weight of 15 mg was used for the analysis. The specimen was first equilibrated at 100o C for 15 mins to achieve fully thermal rejuvenation and then cooled to 40o C. Four different cooling rates were chosen as 10o C/min, 3o C/min, 1o C/min and 0.3o C/min. After held at 40o C for 3 mins, the specimen was heated to 100o C at 10o C/min.

2.3. Uniaxial compression tests The uniaxial compression tests were carried out on the Instron 3367 Series Testing Machine. The temperature was controlled by a thermal chamber. The cylinder-shaped specimens were used. The diameter was 9.5 mm while the height was chosen as 10 mm. The pictures of experimental setups can be found in Fig. S1 in the supplementary material. The specimens were first annealed at 100o C for 30 mins to achieve thermally rejuvenation. The specimens were then cooled to a lower temperature (25o C, 40o C and 55o C) at 3o C/min and annealed for 30 mins before applying uniaxial compression. The specimens were compressed to an engineering strain of 30% or 60% at a constant engineering strain rate and allowed for 15 mins stress relaxation before unloading to zero stress at the same rate of the loading process. The deformation rates were chosen as 0.0001/s, 0.0003/s and 0.001/s. The relative low strain rates were adopted to guarantee an isothermal deformation condition. For the constrained recovery tests, the programed process was the same as described above. Only specimens with a deformation rate of 0.0003/s were chosen for the shape-memory tests. The compression plateau was first moved to contact the specimens. The specimens were then heated to 110o C at a heating rate of 0.5o C/min. The resultant force was recorded in the fixed-strain recovery process. We have adopted low heating rates to minimize the effects of heat conduction.

3. Theory

3.1. Kinematics We begin the theory by defining the deformation gradient F = ∂ x (X)/∂ X, where X is a point in the reference undeformed configuration and x (X) is a point in the current deformed configuration. A decomposition of the deformation gradient F into multiple pairs of elastic and viscous deformation gradients is applied as F = Fej Fvj ,

j = 1..N.

4

(1)

The right and left Cauchy-Green tensors are further defined as, C = FT F,

e Cej = FeT j F j,

b = FFT ,

bej = Fej FeT j .

(2)

Polymers respond differently to volumetric and shear deformation. Thus, the deformation gradient is split into a volumetric and deviatoric part as, e

F = J −1/3 F,

F j = J e-1/3 Fej , j

C = J −2/3 C,

C j = J e-2/3 Cej , j

b=J

−2/3

b,

e

e bj

(3)

= J e-2/3 bej , j

  where J = det (F) and J ej = det Fej . The velocity gradients are defined as ˙ −1 , L = FF

Lvj = F˙ vj Fv-1 j .

(4)

 
The symmetric parts of velocity gradients are D = F + FT /2 and Dvj = Fvj + FvT /2. j

3.2. Thermodynamic framework The constitutive model should satisfy the requirements of the thermodynamics restrictions. Here we extend the nonequilibrium framework in Xiao and Nguyen (2015) to the large strain region by incorporating the strain hardening effect. The following framework is formulated using the thermodynamic potentials which are related to the volume. The first law of thermodynamics is represented as, 1˙ e˙ = S : C − ∇X · Q. 2

(5)

In the above equation, e represents the internal energy, S reprsents the second Piola-Kirchhoff stress tensor, and Q is ˙ represents the power applied to the system. For the the heat flux flowing into the material systems. The term S : 1 C 2

first law, the heat supply by the radiation is omitted. The entropy inequality requires η˙ + ∇X · H ≥ 0,

(6)

in which we denote the entropy of the system as η and the entropy flux from environments to the system as H. In nonequilibrium thermodynamics, it is assumed that the system of amorphous solid consists of a fast kinetic subsystem and a slow configurational subsystem (Nieuwenhuizen 1998, Bouchbinder and Langer 2009). Xiao and Nguyen (2015) replaces the single configurational subsystem with multiple configurational subsystems based on the fact that structural relaxation exhibits a broad distribution of relaxation mechanisms similar to viscoelastic relaxation. Thus, the corresponding internal energy, heat flux, entropy, and entropy flux are all decomposed into a component of the

5

kinetic subsystem and multiple components of the configurational subsystems, shown as follows, P

P

e = ek + ∑ eci , i

Q = Qk + ∑ Qci , i

P k

η = η +∑

(7)

P

ηic ,

k

H = H +∑

i

Hci .

i

As discussed in Kamrin and Bouchbinder (2014), the entropy flux in the kinetic subsystem and the configurational subsystems is defined as, Hk =

Qk , T

Hci =

Qci . Tei

(8)

The term Tei is the effective temperature, which is a thermodynamic state variable representing the degree of disorder in the ith configurational subsystem. In the two-temperature thermodynamic framework, temperature T and kinetic entropy η k are a pair of thermodynamic conjugate variables, while each effective temperature Tei and the configurational entropy ηic form another pair of thermodynamic conjugate variables. Following the work of Nieuwenhuizen (1998), the Helmholtz free energy density is represented as P

Ψ = ek − T η k + ∑ eci − Tei ηic , | {z } i | {z }

(9)

Ψci

Ψk

where Ψk and Ψci are the Helmholtz free energy density contributed by the kinetic and configurational parts respectively. Substituting Eq. (7) and (8) into (6) yields another form of the second law, P

η˙ k + ∑ η˙ ic + i

P P ∇X · Qci ∇X · Qk 1 1 − 2 Qk · ∇X T + ∑ − ∑ 2 Qci · ∇X Tei ≥ 0. T T Tei i i Tei

(10)

The first law can be expanded by combining Eqs. (5) and (9), which gives P P P ˙ − ∇X · Qk − ∑ ∇X · Qci . ˙ + T η˙ k + ∑ Te η˙ ic + η k T˙ + ∑ ηic T˙e = S : 1 C Ψ i i 2 i i i

(11)

Eq. (11) is then substituted into Eq. (10) to eliminate the kinetic entropy η k , resulting in the following form for the entropy inequality, P

P

∑ (T − Tei ) η˙ ic + ∑ (T − Tei ) i

i

P P ∇X · Qci 1 T 1˙ ˙ ≥ 0. − Qk · ∇X T − ∑ 2 Qci · ∇X Tei + S : C − η k T˙ − ∑ ηic T˙ei − Ψ Tei T 2 i Tei i

(12)

The Helmholtz free energy density is assumed to depend on the total and elastic deformation tensors C and Cej , temperature T , effective temperatures Tei , and the internal variables A j representing orientation tensors for strain

6

hardening responses. This brings the following form for the rate of the Helmholtz free energy density, N P N ˙ + ∑ ∂ Ψ : F˙ vj + ∑ ∂ Ψ : A ˙ j. ˙ = ∂ Ψ T˙ + ∑ ∂ Ψ T˙e + 2 ∂ Ψ : 1 C Ψ i v ∂T ∂C 2 i ∂ Tei j ∂Fj j ∂Aj

The above equation can be rewritten considering

∂Ψ ∂ Fvj

(13)

: F˙ vj = −2Cej ∂∂CΨe : Lvj as, j

P N N ˙ − ∑ 2Cej ∂ Ψ : Lvj + ∑ ∂ Ψ : A ˙ j. ˙ = ∂ Ψ T˙ + ∑ ∂ Ψ T˙e + 2 ∂ Ψ : 1 C Ψ i ∂T ∂C 2 ∂ Cej i ∂ Tei j j ∂Aj

(14)

The second law is then evaluated by substituting Eq. (14) into Eq. (12), P N N ∇X · Qci 1 ∂Ψ T ∂Ψ ˙ v − Qk · ∇X T − ∑ 2 Qci · ∇X Tei − ∑ : A j + ∑ 2Cej e : Lj T T T ∂ A ∂ C e j e i j i j j i i i      P  ∂Ψ 1˙ ∂Ψ ˙ ∂Ψ ˙ − ηk + + S−2 : C T − ∑ ηic + Tei ≥ 0. ∂C 2 ∂T ∂ Tei i P

P

∑ (T − Tei ) η˙ ic + ∑ (T − Tei )

(15)

The equilibrium thermodynamic gives, S=2

∂Ψ , ∂C

ηk = −

∂Ψ , ∂T

ηic = −

∂Ψ . ∂ Tei

(16)

Thus, Eq. 15 is reduced to P

P

∑ (T − Tei ) η˙ ic + ∑ (T − Tei ) i

i

N P N ∇X · Qci 1 ∂Ψ ˙ T ∂Ψ − Qk · ∇X T − ∑ 2 Qci · ∇X Tei − ∑ : A j + ∑ 2Cej : Lvj ≥ 0. (17) Tei T ∂ Cej j ∂Aj i Tei j

We adopt the followed form for the evolution equation of the orientation tensor A j (Frederick and Armstrong 2007, Xiao and Tian 2019), ˙ j = A j Dvj + Dvj A j − φ (A j ) , A

(18)

in which the term φ (A j ) represents the relaxation behavior of the molecular orientation. Xiao and Tian (2019) has demonstrated that the relaxation term is necessary to capture the influences of predeformation on the hardening behavior in the reloading process. Eq. (18) is then substituted to Eq. (17), which gives, P

P

P ∇X · Qci 1 T − Qk · ∇X T − ∑ 2 Qci · ∇X Tei Tei T i Tei ! ∂Ψ ∂Ψ 2Cej − 2A j : Dvj ≥ 0. e ∂Cj ∂Aj

∑ (T − Tei ) η˙ ic + ∑ (T − Tei ) i

i

N

N ∂Ψ +∑ : φ (A j ) + ∑ j ∂Aj j

(19)

The above derivation has used the property that both Cej and A j are symmetric vectors, and the free energy Ψ is a

7

function of the isotropic tensors C, Cej and A j , which leads to A j ∂∂AΨj =

∂Ψ ∂Aj Aj

and Cej ∂∂CΨe = j

∂Ψ e ∂ Cej C j .

The last two terms of the right side of Eq. (19) are defined as the dissipative energy due to the mechanical deformation as, W

N

int

N ∂Ψ ∂Ψ ∂Ψ =∑ : φ (A j ) + ∑ 2Cej e − 2A j ∂ A ∂ A ∂ C j j j j j

! : Dvj .

(20)

Combining Eqs. (11), (14), (16) (18) and (20) yields a new expression for the first law, P

P

i

i

T η˙ k + ∑ Tei η˙ ic = −∇X · Qk − ∑ ∇X · Qci + W int .

(21)

The above equation represents the energy balance of the whole system. Here we make several assumptions to obtain the corresponding first law of kinetic and configurational subsystems. First, in addition to the heat exchange within each subsystem, the kinetic subsystem and each configurational subsystem has a heat exchange denoted as Qkc i , while the interaction between different configurational subsystems is negligible and the related heat exchange between different configurational subsystems is ignored (Xiao and Nguyen 2015). Second, the dissipative energy into  all  is distributed int,c ∂ Ψk ∂ Ψk N N ∂ Ψk P e int,k int int,k = ∑ j ∂ A j : φ (A j ) + ∑ j 2C j ∂ Ce − 2A j ∂ A j : the subsystems represented as W = W + ∑i Wi , where W j  c c c ∂ Ψ ∂ Ψ ∂ Ψ int,c N N v e v i i i = ∑ j ∂ A j : φ (A j ) + ∑ j 2C j ∂ Ce − 2A j ∂ A j : D j (Kamrin and Bouchbinder 2014). Based on these D j and Wi j

assumptions, the first law in each subsystem is represented as P

int,k T η˙ k = −∇X · Qk − ∑ Qkc , i +W

(22a)

i

int,c , Tei η˙ ic = −∇X · Qci + Qkc i + Wi

i = 1...P.

(22b)

Combining Eq. (19) and Eq. (22b) to eliminate the ηic term leads to p P P Qkc 1 T T − Qk · ∇X T − ∑ 2 Qci · ∇X Tei + ∑ (T − Tei ) i + W int,k + ∑ Wi int,c ≥ 0. T T T T ei ei ei i i i

(23)

A strong form of the second law is imposed by requiring the value of each term in Eq. (23) is nonnegative, 1 − Qk · ∇X T ≥ 0, T −

T c Q · ∇X Tei ≥ 0, Te2i i

(T − Tei )

Qkc i ≥ 0, Tei

(24a) (24b)

(24c)

W int,k ≥ 0,

(24d)

T W int,c ≥ 0. Tei i

(24e)

8

Eqs. (24a) and (24b) represent the entropy inequality due to heat flux. A Fourier-type conduction equation can be used to satisfy the above requirement as, qk = −kk ∇x T,

qci = −kic ∇x Tei ,

(25)

where qk = 1/JFQ and qci = 1/JFQci are the heat fluxes contributed by the kinetic and configurational subsystem in the current configuration, and the nonnegative variables of kk and kic are the corresponding thermal conductively of each subsystem (Holzapfel 2000). Eq. (24c) indicates that when the configurational subsystem is out of equilibrium (T 6= Tei ), there exists a heat flux between the kinetic subsystem and the configurational subsystem. This heat flux drives the effective temperature toward temperature, which is the underlying physical mechanism for physical aging. This inequality can be satisfied by adopting the following equation for Qkc i , Qkc i = Ki (T − Tei ) ,

i = 1...P,

(26)

where Ki ≥ 0 is a parameter similar to the heat transfer coefficient. Eqs. (24d) and (24e) require that the dissipative energy flowing into either the kinetic and configurational subsystem cannot be negative, which automatically requires W int ≥ 0. This can be satisfied through ∂Ψ ∂Ψ 2Cej − 2A j ∂ Cej ∂Aj

! : Dvj ≥ 0, (27)

∂Ψ : φ (A j ) ≥ 0. ∂Aj

The first part of the above conditions can be satisfied by using the following relation for Dvj (Xiao and Tian 2019), Dvj =

Meff j 2υ j

,

∂Ψ ∂Ψ − 2A j , ∂ Cej ∂Aj | {z } | {z }

e Meff j = 2C j

neq

Mj neq

where M j

j = 1...N,

(28)

Mback j

is the Mandel stress, Mback is the back stress accounting for the strain hardening response, Meff j j is the

effective Mandel stress for the viscous deformation, and υ j > 0 is the viscosity. Figure 1 shows a rheological model of the constitutive theory, which clearly demonstrates the arrangement of different Mandel stresses. The rheological model also provides a mechanism to explain strain hardening. With the deformation, the development of molecular orientation induces a back stress, which further causes a strain hardening response. As discussed in Xiao and Tian (2019), the following form for φ (A j ) is used to meet the requirement of the nonnegative dissipative energy from chain orientation, φ (A j ) =

2 ∂Ψ Aj A j, ∂ Aj υ ori j

(29)

where υ ori j > 0. Eq. 29 implies a nonnegative dissipative of molecular orientation relaxation. Using this simple

9

Meq e

Mneq 1 (C1 )

Meff 1 Mback 1 (A1 )

e Mneq N (CN)

Meff N back

MN (AN) Figure 1: The rheological representative of the model.

relationship, the dissipative energy contributed by the relaxation of chain orientation can be described as, ∂Ψ 1 : φ (A j ) = Mback : Mback ≥ 0. j j ∂Aj 2υ ori j

(30)

The satisfactory of Eqs. (28) and (29) is not a sufficient condition for guaranteeing W int,k ≥ 0 and Wi int,c ≥ 0. We will revisit this issue once the Helmholtz energy density of the kinetic and configurational subsystem is explicitly defined.

3.3. Constitutive relationships The Helmholtz free energy density needs to be explicitly defined to obtain the constitutive relations. Different forms of Helmholtz free energy density can be found in Reese and Govindjee (1997) and Lion et al. (2017). Here we expand the Helmholtz free energy density in Xiao and Nguyen (2015) by incorporating the free energy density contributed by chain orientation with the following form,    eff     µ neq      x  1
λ y e j j  + κ J 2 − 2 log J − 1 tr C j − 3 + µjori λL2  x + ln − y − ln Ψk = (1 − a) ∑   2 λ sinh x λ sinh y 4 L L j N

+ cg (T − T0 ) − cg T ln

T , T0 (31a)

10

   eff     µ neq      x  1

λ Tei µ eq y e j 2 j  Ψci =aφi ∑ tr C j − 3 + µ ori λ x + ln − y − ln + φ tr C − 3 i j L  λL sinh x λL sinh y T0 2 j  2 N

+ φi ∆c (Tei − T0 ) − φi ∆cTei ln

Tei . T2 (31b)

In the above equation, both Ψk and Ψci have a mechanical part and a thermal part. In the mechanical part, the equilibrium rubbery response is caused by entropy change and is incorporatied in Ψci . The equilibrium shear modulus is denoted as µ eq . The bulk response is included in the kinetic part with κ representing the bulk modulus. The nonequilibrum response is distributed in both kinetic and configurational components with 0 ≤ a ≤ 1 as a partition parameter neq

as the nonequilibrium shear modulus. The eight chain model is used for strain hardening with µ ori j as the q  eff  eff chain stiffness, λL as the stretch limit, λ j = trA j /3, x = L −1 λ j /λL and y = L −1 (1/λL ), where L −1 (·) is and µ j

the inverse Langevin function. For the thermal part, we denote cg as the heat capacity of the kinetic subsystem and ∆c as the total heat capacity as the configurational subsystems. Each structural relaxation process takes a fraction of φi with the restriction that 0 < φi ≤ 1 and ∑Pi φi = 1. Finally, T0 is the reference temperature and T2 is the Kauzmann temperature. Substituting Eqs. (31a) and (31a) into (16) gives the explicit form for stress and entropy as P

S = ∑ φi i

  N  
Tei eq −2/3 1 1 κ neq e−2/3 v−1 µ J Fj I − tr (C) C−1 + ∑ µ j J j I − tr Cej Ce−1 Fv−T + (J 2 − 1)C−1 , (32a) j j T0 3 3 2 j η k = cg ln ηic = ∆cφi ln



T T0

 ,



Tei φi µ eq − tr C − 3 . T2 2T0

(32b) (32c)

As shown in Eq. (32c), the thermal part of configurational entropy is zero when the effective temperature equals to the the Kauzmann temperature. The Cauchy stress can then be obtained through σ = 1J FSFT , which is pushing forward of the second Piola-Kirchhoff stress S, as    

1 P Tei eq 1 1 N neq e 1 1 e
σ = ∑ φi µ b − tr b I + ∑ µ j b j − tr b j I + κ J 2 − 1 I. J i T0 3 J j 3 2J neq

The Mandel stresses M j

(33)

and Mback can be obtained as j   1  e e C j − tr C j I , 3   eff   ori µ
λ λ 1 j j L back −1   Mj = L A j − tr A j I . 3 λ eff λL 3 neq

Mj

neq

= µj

(34)

j

The evolution of Dvj is described by Eq. (28). The viscosity of amorphous polymers depends on temperature, structure 11

as well as the flow stress. To connect the viscosity with the temperature and structure, the Adam-Gibbs model (Adam and Gibbs 1965) is adopted, while the Eyring model (Eyring 1936) is adopted to capture the dependence of viscosity on stress, shown as, υ j = υ ref j exp





B T ∑Pi ηic

  −1 Vs s Vs s sinh , T T

(35)

where υ ref state, B is the activation energy of Adam-Gibbs model, VS is j is the characteristic viscosity at the reference q 1 N N eff the activation volume of the Eyring model, and s = 2 ∑ j Meff j : ∑ j M j is the flow stress. The corresponding stress neq

relaxation time can be defined as τS j = υ j /µ j . From the above equation, it can be seen that the stress relaxation time is small compared with the experimental time scale when the temperature is well above Tg . If deformed in this state, amorphous polymers can quickly return to their original sizes after unloading. Thus, the behavior is a viscoelastic response. The stress relaxation time increases tremendously with decreasing the temperature below Tg . However, the stress relaxation time also depends on the flow stress. Here we use the Eyring model to describe the decrease of the relaxation time with the flow stress. So in uniaxial deformation conditions, amorphous glassy polymers exhibit a clear yield point, where the viscous deformation rate is comparable with the applied deformation rate. With the deformation beyond the yield point, a large fraction of the deformed shape cannot recover after unloading. However, different from plastic deformation in metals, the recovery still occurs. But the recovery rate is extremely low so that no observable recovery occurs in the experimental time scale. Thus, the plastic deformation in amorphous polymers is due to extremely high viscosity, which is also named as the viscoplastic response. Combining Eqs. (18), (29), (31a) and (31b) leads to the evolution equation for the orientation tensor A j , ˙ j = A j Dvj + Dvj A j − A

Mback j υ ori j

A j.

(36)

The viscosity-type parameter υ ori j needs to be determined. Experiments have shown that a more pronounced strain hardening response occurs with decreasing temperature and increasing deformation rates. Thus, the Adam-Gibbs model is not suitable because the effective temperature is higher when polymers deformed at a higher strain rate. This can cause a smaller strain hardening response, which is inconsistent with the experimental observations (Guo 2018). The hybrid function in Xiao and Tian (2019) with a Williams-Landel-Ferry (WLF) function above Tg and an Arrhenius function below Tg is used to describe the effects of temperature on υ ori j . Similar to the viscosity function, a ori stress accelerated term is also incorporated into υ ori j . Thus, the following expression is used for υ j ,

  ori ori −1 V s sinh s , T T
g −C T − Tg log (a (T )) = g1 , T ≥ Tg , C2 + T − Tg A A log (a (T )) = − + , T < Tg , T Tg

υ ori j

V ori sori gori = υ j a (T ) s

gori

(37)

denotes the characteristic viscosity for the relaxation of the orientation tensor at Tg , Vsori is the activation q g g volume parameter, sori = 12 ∑Nj Mback : ∑Nj Mback , C1 and C2 are the WLF constants, and A is a pre-factor of the j j

in which υ j

Arrhenius function.

12

To describe the rejuvenation effect, the total dissipative energy is assumed to flow into both kinetic and configurational subsystems. Here we further require that the partition parameter a to be a constant value. Thus, the dissipative energy flowing into the kinetic and configurational subsystems is a constant fraction of the total dissipative energy. Thus, by adopting Eqs. (28) and (29), the requirement for nonnegative values of W int,k and Wi int,c is automatically satisfied. The total dissipative energy can be represented as, W

N int

=∑

Meff j

:

Dvj +

j

! 1 back back Mj : Mj . 2υ ori j

(38)

To obtain the evolution equation of the effective temperature, Eqs. (25), (26), (32c) and (38) is used to rewrite Eq. (22b), which gives,
∆cφi T˙ei = κic ∇X · C−1 ∇X Tei + Ki (T − Tei ) + aφi W int + φi Hi int , where Hi int =

i = 1...P,

(39)

J −2/3 Tei µ eq 2T0


˙ is a latent heat term. Defining the equilibrium part of the second PiolaI − 13 tr (C) C−1 : C
T e Kirchhoff stress tensor as Seq = ∑Pi φi T0i µ eq J −2/3 I − 13 tr (C) C−1 , we can obtain the following relation ∑Pi φi Hi int = ˙ which represents a thermoelastic coupling effect. For a rubber material, this latent heat term means that Seq : 1 C, 2

deforming a rubber band in an adiabatic condition can change its temperature (Holzapfel 2000). Using the same set of equations, the first law of the kinetic subsystem can be written in the following explicit form,
P cg T˙ = κ k ∇X · C−1 ∇X T − ∑ Ki (T − Tei ) + (1 − a) W int .

(40)

i

The first law of the whole system can be rewritten by summing Eqs. (39) and (40) to give P
P
cg T˙ + ∆c ∑ φi T˙ei = κ k ∇X · C−1 ∇X T + ∑ κic ∇X · C−1 ∇X Tei + W int + H int . i

(41)

i

We further neglect the heat conduction term in the configurational subsystems, though this term may be important to capture the evolution of localization behavior in amorphous solids (Manning et al. 2009). The Eq. (39) is then reduced to

where τRi =

T − Tei a 1 + W int + Hi int , T˙ei = τRi ∆c ∆c ∆cφi Ki

i = 1...P,

(42)

is the structural relaxation time. Without mechanical deformation, the above equation is the same

as the first-order evolution equation of fictive temperature proposed by Tool (1946), which is widely employed to describe the structural relaxation behaviors of volume and enthalpy in amorphous solids (Xiao et al. 2017b). Some other approaches without using the effective temperature or fictive temperature have also been developed to model physical aging. For example, Lion et al. (2011) developed a thermodynamic consistent theory based on internal variables, which captures many important aging behaviors. The term

a int ∆c W

indicates that the dissipative energy can

increase the effective temperature and induce a mechanical rejuvenation effect. Buckley et al. (2004) used the fictive temperature as an internal variable, evolving with the viscous strain rate, to model the rejuvenation effects. However, in

13

the two-temperature thermodynamic framework, the effective temperature is treated as a thermodynamic state variable. The evolution equation is obtained through a strict thermodynamic derivation rather than using a phenomenological approach. This is one of the main advantages of using the two-temperature thermodynamic framework to model the nonequilibrium behaviors of amorphous solids compared with the thermodynamic consistent models based on internal variables. Lastly, the latent heat term is typically much smaller than the viscous dissipative energy when polymers are deformed in the glassy state. The Adam-Gibbs model is also applied for the structural relaxation time, τRi = τRrefi exp



B T ∑Pi ηic

 ,

(43)

where τRrefi denotes the characteristic structural relaxation time in the reference state. In summary, the constitutive relations contain Eqs. (28) and (32)-(43). Eqs. (32)-(34) give the explicit expression for the stress, kinetic entropy and configurational entropy. The evolution equations of viscous deformation velocity Dvj and the orientation tensor A j are shown in Eqs. (28) and (36) respectively. Eq. (42) shows that the effective temperature evolves towards equilibrium due to the aging effect and departs from the equilibrium state caused by mechanical rejuvenation. The model totally contains three relaxation mechanisms: the relaxation of viscous deformation, the relaxation of chain orientation and structural relaxation. Each relaxation is characterized by the corresponding viscosity or relaxation time, with the explicit forms listed in Eqs. (35), (37) and (43). In the above model formulation, we do not consider the thermal strain. The procedures about incorporating the thermal strain can be found in Xiao et al. (2013).

4. Parameter Determination

The model is used to simulate the thermo-mechanical responses of PETG. The model parameters for PETG is listed in Tab. 1. The parameters can be divided into four different groups: the viscoelastic parameters, the viscoplastic parameters, the hardening parameters and the parameters related with structural relaxation. In our previous works (Xiao and Nguyen 2015, Xiao et al. 2017a, Xiao and Tian 2019), we have developed several experimental methods to determine the model parameters. Here we provide a brief description about the procedures. The viscoelastic parameters are determined using the dynamic frequency sweep tests described in Sec. 2.1. Using the principal of time-temperature superposition, the master curve is constructed through shifting the storage modulus measured at various temperatures to the reference temperature T0 = 102o C. The experimentally measured data from 72o C to 102o C corresponding to the glass transition region was used to obtain the master curve. The glass transition temperature was chosen as 72o C, representing the onset temperature of glass transition. The master curve at the reference temperature is plotted in Fig. 2. The shifting process also generates the shift factor shown in Fig. 3. The parameters B, T2 of the Adam-Gibbs model can be obtained through fitting to the measured shift factor while assuming g

g

an equilibrium structure state. The WLF constants C1 and C2 in Eq. (37) can be estimated in the same manner. Figure 3 shows that in the equilibrium state the Adam-Gibbs model is nearly identical to the WLF model, both of which can accurately describe the experimentally measured shift factor. The rubbery Young’s modulus E eq was then obtained from the rubbery plateau, which can further give the value of 14

Table 1: Model parameters for amorphous thermoplastic PETG.

Parameter µ eq (MPa) T 0 (K) 

Physical significance equilibrium shear modulus reference temperature

Values 0.93 375

stress relaxation spectrum

shown in Fig. 4-a

κ (MPa) cg (J/(gK)) ∆c  (J/(gK)) 

bulk modulus heat capacity of the kinetic subsystem total heat capacity of configurational subsystems

1166.7 1.09 0.42

structural relaxation spectrum

shown in Fig. 4-b

Tg (K) B(J/g) T2 (K) g C1 g C2 (K) Vs 4(K/MPa) a λL neq µ ori j /µ j gori g τ j /τS j A(K) Vsori (K/MPa)

glass transition temperature thermal activation energy Kauzmann temperature first WLF constant at Tg second WLF constant at Tg activation parameter for viscous flow partition parameter stretch limit neq the ratio between µ ori j and µ j g gori the ratio between τ j and τ j the Arrhenius parameter for molecular orientation relaxation activation parameter for relaxation of network orientation

345 240.2 323.5 11.8 21.5 125 0.67 1.7 0.0114 100.0 15000 150.0

g

neq

τS j , µ j

g

τRi , φi

the equilibrium shear modulus through µ eq = E eq / (2 (1 + υr )) with the rubbery Poisson’s ratio υr = 0.5. The storage modulus in the glassy plateau region gives the value of the glassy Young’s modulus E g . The bulk modulus can be obtained using κ = E g / (3 (1 − υg )) with the glassy Poisson’s ratio υg = 0.35.

Storage modulus (MPa)

103 102 oC-Exp 96oC-Exp 90oC-Exp 87oC-Exp 84oC-Exp 81oC-Exp 78oC-Exp 75oC-Exp 72oC-Exp Model fit

102

101

100

100

103

106

109

Frequency (Hz)

Figure 2: The experimentally measured and model fitted master curve of PETG at the reference temperature 102o C.

The finite deformation model can be linearized to a small strain model followed the procedure described in Nguyen et al. (2010). For a small strain harmonic deformation, the storage modulus can be calculated from the stress relaxation

15

Shift factor

10

6

10

4

10

2

Exp Adam-Gibbs model Fit WLF model fit

10 0 345

360

375

Temperature (K) Figure 3: The shift factor to construct the master curve of the storage modulus. neq

spectrum τS0j − E j

as (Nguyen et al. 2010),  2 neq E j ω 2 τ 0j 0 Edisc (ω) = E eq + ∑  2 . j 1 + ω2 τ0 j N

(44)

However, it is difficult to directly obtain the discrete relaxation spectrum by fitting the master curve. Instead, the nonequilibrium part of storage modulus (E 0 neq = E 0 − E eq ) can be described by a continuous relaxation spectrum h (τ) as E0

neq

ω 2τ 2 h (τ) d ln τ. 2 2 −∞ 1 + ω τ

Z ∞

(ω) =

(45)

Following the work of Tschoegl (1989), a second order approximation of h (τ) is given as, " h (τ) = E

0 neq

d log E 0 neq 1 − d log ω 2



d log E 0 neq d log ω

2

1 d 2 log E 0 neq − 4.606 d (log ω)2

# √ ω= 2/τ

.

(46)

A polynomial equation was first used to fit the noneqilibrium storage modulus, which was then substituted into Eq. (46) to give an explicit polynomial function of h (τ). The cumulative spectrum, representing the total relaxation modulus above a certain relaxation time, is defined as (Haupt et al. 2000),

Z ∞

Hcont (τ) =

h (z) d ln z.

(47)

ln τ

From the above definition, the cumulative spectrum using discrete relaxation processes has the following form as neq

Hdisc (τ) = ∑ E j hτ − τ ref j i, where hzi=1 when z ≤ 0, and hzi=0 when z > 0.

16

(48)

A power law distribution is assigned for the relaxation time as, τ 0j

=τ

max



τ min τ max

j−1  N−1

, j = 1 : N,

(49)

where τ max denotes the upper bound and τ min denotes the lower bound of the stress relaxation time. Following the discussion in Haupt et al. (2000) and Nguyen et al. (2010), the discrete spectrum is calculated as,     1 Hcont τ1ref + Hcont τ2ref , 2     1 neq ref Ej = Hcont τ ref , 1 < j < N − 1, j+1 − Hcont τ j−1 2 neq

E1

neq

=

(50)

N−1

EN = E neq −

∑ E neq j , j

where E neq = E g − E eq . neq

The obtained stress relaxation spectrum τ 0j − E j

together with the equilibrium rubbery modulus is substituted into

Eq. (44) to calculate the storage modulus. As shown in Fig. 2, the predicted master curve agrees with the measured

neq neq data. The nonequilibrium shear modulus can be calculated as µ j = E j / 2 1 + υg . Figure 4-a plots the stress g

neq

g

relaxation spectrum τS j − µ j . The characteristic stress relaxation time τS j is related with τ 0j through g τS j

= τ 0j exp

B B
− P c P c T0 ∑i ηi (Tei = T0 , F = I) Tg ∑i ηi Tei = Tg , F = I

! .

(51)

10 2

10 1

φi

µ neq (MPa) j

10 -1

10 0

10 -1

10

10

-3

10

0

10

3

10

6

-2

10

-3

10

0

10

g τSj (s)

g τRi (s)

(a)

(b)

3

10

6

Figure 4: a) The stress relaxation spectrum and b) the structural relaxation spectrum of PETG.

We then show the procedures to obtain the structural relaxation spectrum from the caloric DSC tests. The enthalpy response in thermal sanning is measured for specimens without mechanical deformation. Thus, we can simplify the

17

constitutive relations in Sec. 3.3 by only considering structural relaxation with the following equations, T − Tei T˙ei = , τRi 

 B

g

τRi = τRi exp 

T ∑Pi ∆cφi ln

−

Tei T2

B T

Tg ∑Pi ∆cφi ln Tg2

(52)

.

The heat flow can be calculated as Q = cg T˙ + ∆c ∑Pi φi T˙ei . Figure 5 shows the DSC response at different cooling rates. No endothermic overshoot is observed at a cooling rate of 10o C/min, while a large endothermic overshoot appears at a cooling rate of 0.3o C/min. The results also show that at the beginning of the heating process, the heat flow gradually increases from zero to a quasi-steady state. Though this process only takes around 20 seconds, it still has a pronounced effect on determining the structural relaxation spectrum. Thus, here an approximation one-dimensional (1D) model is used to describe the heat diffusion process. The heat conduction is assumed to follow Fick’s law as ∂T ∂ 2T =k 2 , ∂t ∂x

(53)

where k is the effective thermal conductively. The boundary condition is set as T (x = 0) = T (t) and

∂T ∂x

(x = l0 ) = 0,

where l0 is the finite length of the specimen in the 1D model. The specific heat flow can be calculated by a meanaverage of the whole specimen as, R

Q=


cg T˙ + ∆c ∑Pi φi T˙ei dx R . dx

(54)

The parameter cg can be estimated from the specific glassy heat capacity, and the configurational heat capacity ∆c can be obtained from the difference of heat capacities in two plateau regions. The parameter k/l02 can be determined from the DSC response at the beginning stage of heating process, which gives a value of 0.047 s−1 . Similar to determining g

the viscoelastic relaxation spectrum, it is difficult to directly fit the spectrum τRi − φi to the DSC response. Here we assume the enthalpy response can be described by the Kohlrausch-Williams-Watts (KWW) model, shown as, Z T

Te = Ti +

τ

g = τR exp

" 1 − exp −

Ti

Z

t(T )

t(T 0 )

B B − Te T ∆c ln T2 Tg ∆c ln TTg 2

1 dt τ !

β !#

dT 0 ,

Te (t = 0) = Ti , (55)

,

where Ti is the initial starting temperature with an equilibrium structure state and Te can be treated as the weighted g

average effective temperature Te = ∑Pi φi Tei . The KWW model contains two parameters: the elaxation time τR and the breadth of relaxation spectrum β . In Xiao et al. (2017b), we have demonstrated that the continuous model using g

a KWW kernel is equivalent to a model with multiple discrete relaxation processes. The relaxation spectrum τRi − φi g

can be calculated when the two parameters τR and β are provided. For convenience, here we also briefly describe the essential steps.

18

Exp Simu

Heat capacity (Wg -1 K -1 )

1

0.5

0

40

Exp Simu

1.5

-1 -1

Heat capacity (Wg K )

1.5

60

80

1

0.5

0

100

40

60

o

(b)

Exp Simu

2

Heat capacity (Wg -1 K -1 )

1.5

-1 -1

Heat capacity (Wg K )

Exp Simu

1

0.5

40

100

Temperature ( C)

(a)

0

80 o

Temperature ( C)

60

80

1.5

1

0.5

0

100

o

40

60

80

100

Temperature (oC)

Temperature ( C) (c)

(d)

Figure 5: Comparing the measured and predicted DSC response of PETG at cooling rates of (a) 10o C/min, (a) 3o C/min, (c) 1o C/min, and (d) 0.3o C/min.

19

Based on Lindsey and Patterson (1980), the continuous relaxation spectrum of the KWW model is given as, hKWW (τ) = −

 β k+1 τ τR ∞ (−1)k sin (πβ k) Γ (β k + 1) , ∑ k! πτ 2 k=0 τR

(56)

where Γ (·) is the gamma function. The discrete spectrum can be represented using Dirac Delta function δ (x) as,   P hdisc (τ) = ∑ φi δ τ − τRgi .

(57)

i

Similarly, the cumulative spectrum can be defined as, Z τ

HKWW (τ) =

0

hKWW (u) du.

(58)

The cumulative spectrum of the spectrum with discrete processes is represented as, P

Hdisc = ∑ φi i

DD EE g τ − τRi ,

(59)

where hhzii=1 when z ≥ 0, and hhzii=0 when z < 0. We also assume a power-law distribution for the discrete structural relaxation time as τRgi = τRmin



τRmax τRmin

 i−1

P−1

,

, i = 1 : P,

(60)

where τRmax and τRmin are the upper and lower limits of the structural relaxation time. The structural relaxation spectrum is then calculated as,     1 g g HKWW τR1 + HKWW τR2 , 2     1 HKWW τRgi+1 − HKWW τRgi−1 , 1 < i < P − 1, φi = 2

φ1 =

(61)

P−1

φP = 1 −

∑ φi . j

g

In order to obtain the structural relaxation spectrum, different combinations of τR and β were chosen to generate the structural relaxation spectrum, which was then used to simulate the DSC response using Eq. (52). The parameters were then determined through a minimization of the following equation, Sdiff =

2 1  exp c − cfit , N0

(62)

where N0 is the total point chosen, cexp is measured heat capacity and cfit is simulated heat capacity. For each cooling rate, an evenly distributed 300 points from 60o C to 90o C were chosen. The final parameter for the characteristic

20

g

relaxation time τR is 3000 s and the breadth of spectrum β is 0.31. The obtained discrete structural relaxation spectrum is shown in Fig. 4-b. Fig. 5 demonstrates that the discrete model can describe the measured enthalpy response, including dependence of the magnitude and location of endothermic overshoot on the cooling rate. The experimental characterization and theoretical analysis for the enthalpy responses are all based on the pressure constant condition. The theory can also been formulated based on the isochoric condition. The related discussion can be found in Lion et al. (2017). The Eyring parameter Vs controls the yield strength. With all the other parameters fixed, decreasing Vs results in a larger yield strength. Thus, this parameter was determined by comparing the simulated and measured yield strength at various temperatures and deformation rates. Using the parameters listed in Tab. 1, the simulated yield strength is compared with the measured results shown in Fig. 6, showing an acceptable agreement. The results also show that at a higher deformation temperature the slope of the yield strength as a function of strain rate is steeper, though the slope almost does not change in simulation. The partition parameter a controls the strain  softening response. Strain hard gori

, µ ori , the stretch limit λL , the Arrhenius j   gori has the parameter A and the activation parameter Vsori . We further assume that the orientation spectrum τ j , µ ori j   g neq gori g neq ori same shape as the stress relaxation spectrum τS j , µ j , which indicates both ratios between τ j /τS j and µ j /µ j ening is controlled by the molecular orientation relaxation spectrum τ j

neq

are constants. The partition parameter a, the modulus ratio µ ori j /µ j stress-strain curve measured at

25o C

and λL were determined through fitting to the

and a strain rate of 0.001/s, as shown in Fig. 7. To determine the viscoplastic

parameters and hardening parameters, an iteration process was used to obtain a good agreement between experimental gori

and simulation results. The parameter A and the relaxation time ratio τ j

g

/τS j were obtained through strain hardening

response at different temperatures. The parameter Vsori describes the stress-activated relaxation behaviors in network orientation. In Xiao and Tian (2019), we assumed that the relaxation of network orientation only depended on temperature. However, we found that a significant strain recovery occurred in the stress-free state before subject to the constrained recovery, which was not consistent with experimental observations. Thus, following the works of Dupaix and Boyce (2007) and Guo (2018), we also incorporate the stress-activated relaxation in chain orientation. The parameter Vsori was then determined by matching the final recovery strain of specimens deformed at 25o C and with 60%

Yield strength (MPa)

programmed strain between experiments and simulation.

45

25oC-exp 25oC-simu 40oC-exp 40oC-simu 55oC-exp 55oC-simu

35

25 10 -4

10 -3

0.0003 -1

Strain rate (s ) Figure 6: Comparing the experimental and simulation results of the yield strength.

21

Exp Simu

True stress (MPa)

45

30

15

0

0

0.3

0.6

0.9

True strain Figure 7: The stress response of PETG deformed at 25o C and 0.001/s.

5. Results

This paper mainly investigates the constrained recovery behaviors of amorphous SMPs. A typical shape-memory cycle is composed of a programming step and a recovery step. The former is to obtain the temporary shape, which is achieved through uniaxial compression below Tg in this work. Thus, in the following, we first compare the experimental and simulation results on the stress response related with the programming step, including constant strain rate uniaxial compression and stress relaxation. The model is further used to simulate the stress response as a function of the recovery temperature in the fixed-strain recovery process. A comparison between the experimentally measured and theoretical predicted resultant stress is also provided. The resultant stress originates from the shape recovery, which is related with the shape-memory effect. However, it should be noticed that polymers with only a partial recovery ability can also generate a stress during the heating process. Thus, we also performed a free recovery test to demonstrate the shape-memory performance of PETG with the detailed information shown in the appendix part. To simulate the mechanical response of the amorphous thermoplastics in the isothermal uniaxial compression and constrained recovery tests, the highly nonlinear relations were discretized in time using finite difference method and solved by the Newton-Raphson method. A detailed explanation of the numerical procedure is provided in the appendix part. Using the parameters listed in Tab. 1, the model is first used to describe the stress response. Fig. 8 compares the measured and simulated stress response at engineering strain rates of 0.001/s and 0.0001/s. As shown, the model captures the influence of temperature and deformation rate on the measured stress response, including a decrease of yield strength and steady-state flow stress with increasing temperature and decreasing rates. The experimental results show that more strain softening occurs with decreasing temperature and increasing strain rates, which is accurately predicted by the model. As summarized in Medvedev and Caruthers (2016), accurately capturing the effects of strain rate on magnitude of strain softening is an important criteria for evaluating the constitutive models for the glassy amorphous polymers. Medvedev and Caruthers (2016) demonstrated that several models (Boyce et al. 1988, van Breemen et al. 2011) failed to capture this feature. In contrast, the model developed by Anand et al. (2009) and Medvedev and Caruthers (2013) can reproduce the rate-dependent strain softening response. In simulation, we have assumed an isothermal condition by ignoring the heat conduction. This is because a relatively low deformation rate is adopted. The assumption of the isothermal deformation is also consistent with the experimental findings in Arruda

22

et al. (1995) and Ames et al. (2009). 45 Exp Simu

45

Exp Simu

o

25 C

o

25 C o

40 C

True stress (MPa)

True stress (MPa)

40 oC o

55 C 30

15

0

0

0.3

0.6

55 C

15

0

0.9

True strain

o

30

0

0.3

0.6

0.9

True strain

(a)

(b)

Figure 8: The stress response of PETG deformed at engineering strain rates of (a) 0.001/s and (b) 0.0001/s.

The strain hardening response is also accurately captured by the model. The model employs a relaxation mechanism for strain hardening. Thus, at a higher temperature, more relaxation occurs, which leads to a lower hardening modulus. Various models (Srivastava et al. 2010, Bouvard et al. 2013) used a temperature-dependent hardening modulus to describe strain hardening. However, Xiao and Tian (2019) have demonstrated that the approach using a temperaturedependent hardening modulus fails to capture hardening in the reloading process, while the approach containing a relaxation mechanism is able to describe the experimental results. Fig. 8 also shows that some discrepancies between experimental and simulation results exist. For example, the model predicts a larger stress than the measured data in the post-yield region for polymers deformed at 55o C. We next employ the model to simulate the stress response during the programming process. The temporary shape was obtained through compressing the specimens to a certain strain and then held for 15 mins for stress relaxation. Fig. 9 shows the corresponding stress response from experiments and simulation during the loading, holding and unloading processes. In general, the model accurately predicts all the characteristics observed in experiments. However, the model significantly underestimates stress relaxation. For all the conditions, the total stress drop is between 14.5 MPa and 15.9 MPa from experiments, while the stress drop predicted by the model is from 8.7 MPa to 9.9 MPa. Stress relaxation depends on relaxation time and the breadth of relaxation spectrum. In Xiao et al. (2017c), we show that the breadth of spectrum has a limited influence on the viscoplastic stress response during a constant strain rate loading process. In contrast, a narrower relaxation spectrum can induce more pronounced stress relaxation in the holding period. Some recent experimental works (Lee et al. 2009, Hebert et al. 2015) demonstrated that the breadth of the relaxation spectrum decreased with the deformation, which was attributed to a decrease in the dynamic heterogeneity. This is a big challenge for nearly all the current constitutive models for glassy polymers in the literature, which assume that the relaxation spectrum remains unchanged. Mathiesen et al. (2014) developed a viscoplastic model incorporating a single relaxation process for intermolecular relaxation and a single relaxation process for network orientation. This model can capture the stress relaxation response of poly(methyl methacrylate) spanning the transition region. However, using a single relaxation process cannot accurately predict the shape-memory behavior which closely relates with the broad distribution of relaxation mechanisms. The stochastic constitutive model developed by Medvedev and Caruthers 23

(2013, 2015) included the nanoscale fluctuations into the constitutive relations and the model assumed the macroscopic response was an assemble average of the whole systems. This model successfully captured the change in the shape of the viscoelastic relaxation spectrum during deformation, which may have potential to accurately simulate the stress relaxation response. However, strain hardening is not incorporated into the framework. The computational complexity also limits its applications. 45

45 Exp Simu

Exp Simu

o

o

25 C

True stress (MPa)

True stress (MPa)

25 C o

40 C

30

55 oC 15

0

0

500

1000

1500

55 oC 15

0

2000

Time (s)

40 oC

30

0

500

1000

1500

2000

Time (s)

(a)

(b)

Figure 9: Comparing the measured and simulated stress response in the programming process.

Though the model fails to capture the magnitude of stress relaxation, it successfully predicts the effects of strain rates on stress relaxation. Fig. 10 shows the relaxation response of PETG from both experimental and simulation results, which indicates that more relaxation occurs for a larger loading rate for both specimens with programmed strains of 30% and 60%. The constitutive model can provide a physical mechanism to explain this observation. Both aging and rejuvenation effects coherently occur in amorphous polymers. With aging, the polymer structure approaches a more sluggish equilibrium state, resulting in an increase in viscosity. Upon deformation, the dissipative energy rejuvenates the polymers and causes a decrease in viscosity, which leads to the observed strain softening phenomenon (Jiang et al. 2015). When deformed at a higher strain rate, more energy is generated to rejuvenate the configurational subsystem and less aging occurs, both of which can produce a more mobile structure state. As shown in Fig. 11, a higher weighted average effective temperature Te , can be observed with a higher loading rate in both loading and holding processes. This more mobile state can lead to a larger relaxation rate in the holding period. Fig. 11 (b) also shows that during the relaxation process the effective temperature decreases, which is caused by that the aging effect surpasses the rejuvenation effect when deformation ceases. The programmed polymers were then subject to the fixed-strain recovery tests. The resulting stress response is shown in Fig. 12. In all conditions, the stress first increases with temperature until reaching a peak value and then decreases to zero at around 100o C. Decreasing the programming temperature leads to a larger peak stress. The programming strain also affects the magnitude of the peak stress. For example, when deformed at 25o C, the peak stress increases from 12.8 MPa to 16.6 MPa when the programmed strain changes from 30% to 60%. This is because more energy is stored in polymers when deformed at a lower temperature or to a larger strain. Several works (Nguyen et al. 2008, Abishera et al. 2018, Gall et al. 2005) have also investigated the constrained recovery tests. The maximum recovery stress in these works is around 3 MPa. This is because the temporary shapes of SMPs in the above works are obtained 24

35

45 0.001/s 0.0003/s 0.0001/s

True stress (MPa)

True stress (MPa)

0.001/s 0.0003/s 0.0001/s

25

15

0

300

600

35

25

900

0

300

Time (s)

600

(a)

(b)

35

45 0.001/s 0.0003/s 0.0001/s

True stress (MPa)

0.001/s 0.0003/s 0.0001/s

True stress (MPa)

900

Time (s)

30

40

35

25

0

300

600

30

900

0

300

Time (s)

600

(c)

(d)

35

45 0.001/s Exp 0.001/s Simu 0.0001/s Exp 0.0001/s Simu

True stress (MPa)

True stress (MPa)

0.001/s Exp 0.001/s Simu 0.0001/s Exp 0.0001/s Simu

25

15

900

Time (s)

0

300

600

35

25

900

Time (s)

0

300

600

900

Time (s)

(e)

(f)

Figure 10: The influence of strain rates on the stress relaxation behavior of PETG deformed at 25o C: measured response with strains of (a) 30% and (b) 60%, and simulated response with strains of (c) 30% and (d) 60%, and a direct comparison of measured and predicted stress relaxation response with strains of (e) 30% and (f) 60%

.

25

0.001/s 0.0003/s 0.0001/s

353

0.001/s 0.0003/s 0.0001/s

353

Te (K)

Te (K)

350

347

351

344

341

0

0.3

0.6

349

0.9

True strain

0

300

600

900

Time (s)

(a)

(b)

Figure 11: The evolution of the effective temperature (a) in the loading process as a function of strain and (b) in the relaxation process as a function of time.

by deforming the polymers above or around Tg . Here we demonstrate that the recovery stress can be significantly increased to more than 10 MPa by cold programming, which can potentially promote the broad applications of SMPs. The location of the peak stress shifts to a higher temperature with increasing programming temperature and strain. The curves of the recovery stress overlap in the high temperature region. All these important features have been successfully captured by the model, including the change of magnitude and location of the peak stress with temperature and programming strain. There are also several inconsistencies between experiments and simulation. For example, the model significantly overestimates the stress response of polymers programmed at 55o C and with a strain of 30%. This is because the model fails to accurately predict the stress relaxation behaviors especially at 55o C. The initial slope of the recovery stress predicted in the model is also stiffer for the programmed strain of 60%. The predicted final stress is not zero, which is inconsistent with the experimental observations. This is because an equilibrium rubbery modulus is assumed, while for thermoplastics the plateau region is only a quasi-equilibrium state associating with the relaxation processes with a large relaxation time. These relaxation processes cannot be captured by the dynamic frequency sweep tests (Collins et al. 2016). In Xiao et al. (2016), we showed that the relaxation behaviors in the glass transition of thermoplastics can be captured by the stress relaxation tests. However, PETG exhibits a strong volumetric relaxation in the glass transition region, which prohibits using the small strain stress relaxation tests to get the viscoelastic relaxation spectrum.

6. General Remarks

We have developed a constitutive model to describe the constrained recovery performance of amorphous thermoplastic PETG with the temporary shape programmed below Tg . The uniaxial compression tests were used to verify the theoretic framework. For uniaxial compression tests, we chose a cylinder-shaped specimen with a ratio of height and diameter close to 1 to minimize the buckling and bulging effects. To further reduce the friction, the lubricant oil was used between the specimens and compression plates. No observable bulge appeared during all the tests. Though the uniaxial compression tests are most widely employed to demonstrate the prediction ability for viscoplastic models (Dupaix and Boyce 2007, Ames et al. 2009, Bouvard et al. 2013), it is still intriguing to show the performance of 26

18 25oC-Exp 25oC-Simu 40oC-Exp 40oC-Simu 55oC-Exp 55oC-Simu

9

6

3

0 20

25oC-Exp 25oC-Simu 40oC-Exp 40oC-Simu 55oC-Exp 55oC-Simu

15

True stress (MPa)

True stress (MPa)

12

12 9 6 3

50

80

0 20

110

o

50

80

110

Temperature (oC)

Temperature ( C) (a)

(b)

Figure 12: The resultant stress during the constrained recovery process of specimens programmed with strains of (a)30% and (b) 60%.

the model to describe the tension deformation. However, due to strain softening, localization can easily occur during the uniaxial tension tests for PETG as shown in Zhang et al. (2018). Thus, it is difficult to obtain the stress-strain relationship in the tension tests. However, the recent development of the experimental technique, especially digital image correlation (DIC), has provided a useful tool to measure the inhomogeneous strain field. Thus, here we use the DIC technique to extract the stress-strain relationship in uniaxial tension tests. A detailed experimental description can be found in the supplementary material. Fig. 13-a shows the change of engineering strain in the loading direction of the middle part of the specimens with time, which is also used in simulation. Fig. 13-b compares the simulated and measured stress response in the uniaxial tension tests. In general, the simulations reasonably capture the experimental results with some discrepancies. First, the simulation overestimates the yield strength. This is because the model does not consider the effect of hydrostatic pressure on the viscosity (Arruda et al. 1995). Second, the model also does not fully capture the hardening responses in the very large strain region. This is because a dissipative mechanism has been adopted for the strain hardening term, which is necessary to capture the strain recovery in the unloading process. In very large strain, the dissipative should be suppressed due to chain alignments, which can be incorporated by using the approach proposed by Dupaix and Boyce (2007). We have shown that the model can describe the stress response in the uniaxial loading conditions including compression and tension. Here we also perform the simulation to investigate the stress-strain relationships in plane strain and biaxial loading conditions to demonstrate that the model has the ability to predict the mechanical behaviors in the three-dimensional deformation. For plane strain, the deformation is applied at an engineering strain rate of 0.001/s. As shown in Fig. 14, the yield strength and strain hardening significantly increase compared with the uniaxial compression tests at the same loading rates, which is consistent with the experimental results in Dupaix and Boyce (2007). For biaxial tension tests, the deformation is applied in the e1 direction at an engineering strain rate 0.001/s while the loading rate in the e3 direction is varied from 0.00025/s to 0.001/s. The stress responses in the two loading directions are plotted in Fig. 15. Fig. 15-a shows that a larger loading rate in the e3 direction results in a lower yield strength and a larger hardening response in the e1 direction. The lower yield strength is caused by more rejuvenation, while a larger stretch in the e3 direction also induces more chain alignments and leads to a stronger hardening response. Due to the limited data available in the literature on biaxial loading conditions, these findings need to be verified in the future. The model has not been implemented for finite element analysis, which does not allow us to evaluate the performance

27

1.8 Exp Simu

True stress (MPa)

Engineering strain

60

1.2

0.6

0

0

50

100

150

40

20

0

200

0

0.3

Time (s)

0.6

0.9

True strain

(a)

(b)

Figure 13: (a) The engineering strain as a function of time in the middle part of the tensile specimen and (b) comparison between the simulated and measured stress response in uniaxial tension tests.

of the model to predict complex deformations at the current stage. This is our on-going work. 75

True stress (MPa)

Uniaxial compression Plane strain

50

25

0

0

0.3

0.6

0.9

True strain Figure 14: Comparison between the predicted stress response in plane strain tests and uniaxial compression tests.

7. Conclusion

We have developed a constitutive model to describe the thermo-mechanical response of amorphous SMPs, which incorporates the mechanisms of stress relaxation, structural relaxation and strain hardening. The effective temperature is introduced as a thermodynamic state variable to represent the state of the nonequilibrium structure, which is coupled with viscoplastic deformation to account for the rejuvenation-induced strain softening. Distinct from the models using a non-dissipative orientation process for strain hardening, the framework includes the relaxation mechanism into network orientation to capture the temperature and rate dependent strain hardening in large deformation. A nonequilibriumm thermodynamic framework is provided to derive the constitutive relations. A series of thermo-

28

60

60 0.001/s 0.0005/s 0.00025/s

True stress σ3 (MPa)

40

1

True stress σ (MPa)

0.001/s 0.0005/s 0.00025/s

20

0

0

0.25

0.5

40

20

0

0.75

0

0.25

True strain ǫ 1

0.5

0.75

True strain ǫ 3

(a)

(b)

Figure 15: The stress response in biaxial tension tests with an engineering strain rate 0.001/s in the e1 direction and an engineering strain rate 0.001/s, 0.0005/s or 0.00025/s in the e3 direction: (a) the stress response in the e1 direction and (b) the stress response in the e3 direction.

mechanical tests are designed to obtain the model parameters related with the stress relaxation spectrum, structural relaxation spectrum, viscoplastic deformation and strain hardening behaviors. The model is then used to simulate the stress-strain relationship and stress relaxation response in uniaxial compression tests. The model provides a reasonable description of the stress response, but predicts a less pronounced stress relaxation compared with experimental observations. One possible reason for the underestimation is that the deformation reduces the dynamic heterogeneity while most of the current constitutive modeling approaches assume that the breadth of the relaxation spectrum remains unchanged during the deformation. Experimental results also show that amorphous polymers with a larger deformation rate exhibit a lower stress during the holding period, which is accurately predicted by the model. The model also captures the effects of deformation temperatures and programming strains on the stress response in the constraint recovery process including the magnitude and localization of the peak stress. Both experiments and simulations show that a large recovery stress can be achieved by cold programming. Thus, the model developed in this work can be a useful tool to promote the applications of amorphous SMPs. We are currently implementing the above thermo-mechanical coupled theory into a finite element suite to investigate the responses of amorphous polymers at high loading rates. The finite element model will also be employed to characterize the onset and propagation of the strain localization behaviors in glassy polymers.

Acknowledgments

R. Xiao acknowledges the funding support from the National Natural Science Foundation of China (Grant Nos. 11872170,11828201) and the One-hundred Talents Program of Zhejiang University.

29

Appendix A. Free Recovery Test

In order to evaluate the shape-memory performance of PETG, a free recovery test was performed on the DMA. A film specimen with a dimension of 20 mm× 5 mm× 1.59 mm was used. The specimen was first stretched to an engineering strain of 0.3 at 80o C and an engineering strain rate of 0.003/s. The specimen was then cooled to 25o C at 5o C/min before unloaded to zero force. To achieve shape recovery, the specimen was heated to 100o C at 1o C/min. Fig. 16 shows that the simulation exhibits a good agreement with the measured strain recovery, which further validates the model. A discrepancy at the high temperature region exists because an equilibrium rubbery modulus is assumed for the thermoplastic PETG. The results also show that the final recovery ratio of PETG is above 95%, demonstrating a good shape memory behavior of this engineering polymer.

Engineering strain

0.3

Exp Simu

0.2

0.1

0 25

50

75

100

Temperature (oC) Figure 16: Comparison between the measured and predicted free recovery performance of PETG.

Appendix B. Numerical Algorithm

For the uniaxial deformation, the deformation gradient can be written as F = λ1 e1 ⊗ e1 + λ2 e2 ⊗ e2 + λ3 e3 ⊗ e3 , while we assume that the deformation is applied in the e2 direction. The model can be implemented in two modes, namely a strain-control mode and a stress control model. In the strain control mode, the value of λ2 is given, while in the stress mode, the value of σ2 is provided. In the following, we use the strain-control mode to explain the implementation process.     We first define the Hencky strain ε = 12 ln (b), ε ej = 12 ln bej and ε vj = 12 ln bvj . In the principal direction, the Hencky e and ε v = ln λ v , where j = 1 : N and k = 1 : 3. The log strain strain can be represented as εk = ln (λk ), ε ejk = ln λ jk jk jk

of each orientational tensor is also defined in the similar way as ε aj = 21 ln (A j ), which has three components in the principal directions. The objective of the numerical algorithm is to obtain the values of all the variables at the n + 1 time step based all the values provided at the n time step. The variables Tei , ε ej and ε aj can be obtained through solving Eqs. (42), (28) and (36). The finite difference method together with the predictor/corrector scheme (Reese and Govindjee 1998) is used,

30

which gives: r1i = (Tei )n+1 − (Tei )n −

T n+1 − (Tei )n+1 (τRi )n+1

∆t −

a int 1 W ∆t − Hi int ∆t = 0 ∆c ∆c

i = 1 : P,

 n+1
tr eff M − ε ej = 0, j = 1 : N, j n+1 2 (υ j )  n+1
n
n+1
n+1
tr ∆t back M − ε aj = 0, j = 1 : N, r3 j = ε aj + ε ej − ε ej +  n+1 j 2 υ ori j r2 j = ε ej


n+1

+

∆t

(63a)

(63b) (63c)

 tr  n where ε ej = (ε j )n+1 − ε vj is the trial solution of the elastic Hencky strain. Eq. (63a) consists of P equations, while both Eqs. (63b) and (63c) contain 3N equations. These highly nonlinear equations are solved using the NewtonRaphson method. The value of ε1 and ε3 at the n + 1 time step can be determined by solving σ1 = 0 and σ3 = 0. Using Taylor expansion, this relationship can be written as: e+ σ

N e ∂ σe ∂σ ∆e ε + ∑ e ∆ε ej = 0, ∂ ε ∂e ε j j

(64)

where σe = [σ1 , σ3 ]T and e ε = [ε1 , ε3 ]T . To find the value of e ε , one central task is building the relationship between ∆ε ej and ∆e ε , which can be obtained through Eqs. (63b) and (63c). In brief, the Eqs. (63b) and (63c) are treated by the variation method, which gives, N ∂e r ∂e r ∆ε + ∑ in ∆ε in j = 0, ∂ε j ∂εj

(65)

h iT e a . From the above equation, the incremental ∆ε e can be represented as a function where e r = [r2, r3]T and ε in j = εj ,εj j of ∆ε. The interested readers can find more information about the implementation of finite deformation viscous model in Reese and Govindjee (1998).

References Abishera, R., Velmurugan, R., Nagendra Gopal, K., 2018. Free, partial, and fully constrained recovery analysis of cold-programmed shape memory epoxy/carbon nanotube nanocomposites: Experiments and predictions. Journal of Intelligent Material Systems and Structures 29 (10), 2164–2176. Adam, G., Gibbs, J. H., 1965. On the temperature dependence of cooperative relaxation properties in the glass-forming liquids. J. Chem. Phys. 43, 139–146. Ames, N. M., Srivastava, V., Chester, S. A., Anand, L., 2009. A thermo-mechanically coupled theory for large deformations of amorphous polymers. part ii: Applications. International Journal of Plasticity 25 (8), 1495–1539. Anand, L., Ames, N. M., Srivastava, V., Chester, S. A., 2009. A thermo-mechanically coupled theory for large deformations of amorphous polymers. part i: Formulation. International Journal of Plasticity 25 (8), 1474–1494. Arruda, E. M., Boyce, M. C., 1993. A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials. J. Mech. Phys. Solids 41, 389–412. 31

Arruda, E. M., Boyce, M. C., Jayachandran, R., 1995. Effects of strain rate, temperature and thermomechanical coupling on the finite strain deformation of glassy polymers. Mechanics of Materials 19 (2-3), 193–212. Belmonte, A., Russo, C., Ambrogi, V., Fern´andez-Francos, X., De la Flor, S., 2017. Epoxy-based shape-memory actuators obtained via dual-curing of off-stoichiometric “thiol–epoxy” mixtures. Polymers 9 (3), 113. Boatti, E., Scalet, G., Auricchio, F., 2016. A three-dimensional finite-strain phenomenological model for shapememory polymers: Formulation, numerical simulations, and comparison with experimental data. International Journal of Plasticity 83, 153–177. Bouchbinder, E., Langer, J. S., 2009. Nonequilibrium thermodynamics of driven amorphous materials. ii. effectivetemperature theory. Phys. Rev. E. 80, 031132. Bouvard, J.-L., Francis, D. K., Tschopp, M. A., Marin, E., Bammann, D., Horstemeyer, M., 2013. An internal state variable material model for predicting the time, thermomechanical, and stress state dependence of amorphous glassy polymers under large deformation. International Journal of Plasticity 42, 168–193. Boyce, M. C., Parks, D. M., Argon, A. S., 1988. Large inelastic deformation of glassy polymers. part i: rate dependent constitutive model. Mechanics of Materials 7, 15 – 33. Buckley, C., Dooling, P., Harding, J., Ruiz, C., 2004. Deformation of thermosetting resins at impact rates of strain. part 2: constitutive model with rejuvenation. Journal of the Mechanics and Physics of Solids 52 (10), 2355–2377. Chevalier, J., Brassart, L., Lani, F., Bailly, C., Pardoen, T., Morelle, X. P., 2018. Unveiling the nanoscale heterogeneity controlled deformation of thermosets. Journal of the Mechanics and Physics of Solids 121, 432–446. Collins, D. A., Yakacki, C. M., Lightbody, D., Patel, R. R., Frick, C. P., 2016. Shape-memory behavior of high-strength amorphous thermoplastic poly (para-phenylene). Journal of Applied Polymer Science 133 (3), 42903. Das, S., Chowdhury, S. R., Roy, D., 2018. A constitutive model for thermoplastics based on two temperatures. European Journal of Mechanics-A/Solids 72, 440–451. Das, S., Roy, D., 2019. A constitutive model for block-copolymers based on effective temperature. International Journal of Mechanical Sciences, 105082. Dupaix, R. B., Boyce, M. C., 2007. Constitutive modeling of the finite strain behavior of amorphous polymers in and above the glass transition. Mechanics of Materials 39, 39–52. Engqvist, J., Wallin, M., Ristinmaa, M., Hall, S. A., Plivelic, T. S., 2016. Modelling multi-scale deformation of amorphous glassy polymers with experimentally motivated evolution of the microstructure. Journal of the Mechanics and Physics of Solids 96, 497–510. Eskandari, A. H., Baghani, M., Sohrabpour, S., 2018. A time-dependent finite element formulation for thick shape memory polymer beams considering shear effects. International Journal of Applied Mechanics 10 (04), 1850043. Eyring, H., 1936. Viscosity, Plasticity, and Diffusion as Examples of Absolute Reaction Rates. J. Chem. Phys. 4, 283–291. Falk, M. L., Langer, J. S., 1998. Dynamics of viscoplastic deformation in amorphous solids. Physical Review E 57 (6), 7192.

32

Falk, M. L., Langer, J. S., 2011. Deformation and failure of amorphous, solidlike materials. Annu. Rev. Condens. Matter Phys. 2 (1), 353–373. Fan, J., Li, G., 2018. High enthalpy storage thermoset network with giant stress and energy output in rubbery state. Nature communications 9 (1), 642. Fisher, J. G., Sparks, E. A., Khan, F. A., Dionigi, B., Wu, H., Brazzo III, J., Fauza, D., Modi, B., Safranski, D. L., Jaksic, T., 2015. Extraluminal distraction enterogenesis using shape-memory polymer. Journal of pediatric surgery 50 (6), 938–942. Frederick, C. O., Armstrong, P., 2007. A mathematical representation of the multiaxial bauschinger effect. Materials at High Temperatures 24 (1), 1–26. Gall, K., Yakacki, C. M., Liu, Y., adn N. Willet, R. S., Anseth, K. S., 2005. Thermomechanics of the shape memory effect in polymers for biomedical applications. Journal of Biomedical Materials Research: Part A 73, 339–348. Ge, Q., Serjouei, A., Qi, H. J., Dunn, M. L., 2016. Thermomechanics of printed anisotropic shape memory elastomeric composites. International Journal of Solids and Structures 102, 186–199. Gu, J., Leng, J., Sun, H., Zeng, H., Cai, Z., 2019. Thermomechanical constitutive modeling of fiber reinforced shape memory polymer composites based on thermodynamics with internal state variables. Mechanics of Materials 130, 9–19. Guo, J., 2018. Modeling the mechanical behavior of soft active materials. Ph.D. thesis, Johns Hopkins University. Haupt, P., Lion, A., Backhaus, E., 2000. On the dynamic behaviour of polymers under finite strains: constitutive modelling and identification of parameters. Int. J. Solids Struct. 37 (26), 3633–3646. He, Y., Guo, S., Liu, Z., Liew, K., 2015. Pattern transformation of thermo-responsive shape memory polymer periodic cellular structures. International Journal of Solids and Structures 71, 194–205. Hebert, K., Bending, B., Ricci, J., Ediger, M., 2015. Effect of temperature on postyield segmental dynamics of poly (methyl methacrylate) glasses: Thermally activated transitions are important. Macromolecules 48 (18), 6736–6744. Holzapfel, G. A., 2000. Nonlinear solid mechanics: a continuum approach for engineers. John Wiley and Sons, LTD, Chichester. Jiang, H., Jiang, C., 2019. Finite deformation constitutive model for macro-yield behavior of amorphous glassy polymers with a molecular entanglement-based internal-state variable. International Journal of Mechanical Sciences, 105064. Jiang, H., Zhang, J., Yang, Z., Jiang, C., Kang, G., 2017. Modeling of competition between shear yielding and crazing in amorphous polymers’ scratch. International Journal of Solids and Structures 124, 215–228. Jiang, M., Wilde, G., Dai, L., 2015. Origin of stress overshoot in amorphous solids. Mechanics of Materials 81, 72–83. Johnsen, J., Clausen, A. H., Grytten, F., Benallal, A., Hopperstad, O. S., 2019. A thermo-elasto-viscoplastic constitutive model for polymers. Journal of the Mechanics and Physics of Solids 124, 681–701. Kamrin, K., Bouchbinder, E., 2014. Two-temperature continuum thermomechanics of deforming amorphous solids. J. Mech. Phys. Solids 73, 269–288. 33

Kontou, E., 2016. Lower and higher strain regime modeling of cyclic viscoplastic response of an amorphous glassy polymer. International Journal of Solids and Structures 97, 489–495. Krairi, A., Doghri, I., Schalnat, J., Robert, G., Van Paepegem, W., 2019. Thermo-mechanical coupling of a viscoelastic-viscoplastic model for thermoplastic polymers: Thermodynamical derivation and experimental assessment. International Journal of Plasticity 115, 154–177. Langer, J., 2004. Dynamics of shear-transformation zones in amorphous plasticity: Formulation in terms of an effective disorder temperature. Physical Review E 70 (4), 041502. Lee, H.-N., Riggleman, R. A., de Pablo, J. J., Ediger, M., 2009. Deformation-induced mobility in polymer glasses during multistep creep experiments and simulations. Macromolecules 42 (12), 4328–4336. Li, G., Wang, A., 2016. Cold, warm, and hot programming of shape memory polymers. Journal of Polymer Science Part B: Polymer Physics 54 (14), 1319–1339. Li, J., Kan, Q., Chen, K., Liang, Z., Kang, G., 2019. In situ observation on rate-dependent strain localization of thermo-induced shape memory polyurethane. Polymers 11 (6), 982. Li, J., Kan, Q.-h., Zhang, Z.-b., Kang, G.-z., Yan, W., 2018. Thermo-mechanically coupled thermo-elasto-viscoplastic modeling of thermo-induced shape memory polyurethane at finite deformation. Acta Mechanica Solida Sinica 31 (2), 141–160. Li, Y., He, Y., Liu, Z., 2017. A viscoelastic constitutive model for shape memory polymers based on multiplicative decompositions of the deformation gradient. International Journal of Plasticity 91, 300–317. Lindsey, G. P., Patterson, G. D., 1980. Detailed comparison of the williams-watts and cole-davidson functions. J. Chem. Phys. 73, 3348–3357. Lion, A., Mittermeier, C., Johlitz, M., 2017. Heat capacities and volumetric changes in the glass transition range: a constitutive approach based on the standard linear solid. Continuum Mechanics and Thermodynamics 29 (5), 1061–1079. Lion, A., Peters, J., Kolmeder, S., 2011. Simulation of temperature history-dependent phenomena of glass-forming materials based on thermodynamics with internal state variables. Thermochimica acta 522 (1-2), 182–193. Lu, H., Wang, X., Xing, Z., Fu, Y.-Q., 2019. A cooperative domain model for multiple phase transitions and complex conformational relaxations in polymers with shape memory effect. Journal of Physics D: Applied Physics 52 (24), 245301. Lu, H., Wang, X., Yao, Y., Fu, Y. Q., 2018. A “frozen volume” transition model and working mechanism for the shape memory effect in amorphous polymers. Smart Materials and Structures 27 (6), 065023. Manning, M., Daub, E., Langer, J., Carlson, J., 2009. Rate-dependent shear bands in a shear-transformation-zone model of amorphous solids. Physical Review E 79 (1), 016110. Mao, Y., Chen, F., Hou, S., Qi, H. J., Yu, K., 2019. A viscoelastic model for hydrothermally activated malleable covalent network polymer and its application in shape memory analysis. Journal of the Mechanics and Physics of Solids 127, 239–265.

34

Mao, Y., Robertson, J. M., Mu, X., Mather, P. T., Qi, H. J., 2015. Thermoviscoplastic behaviors of anisotropic shape memory elastomeric composites for cold programmed non-affine shape change. Journal of the Mechanics and Physics of Solids 85, 219–244. Mathiesen, D., Vogtmann, D., Dupaix, R. B., 2014. Characterization and constitutive modeling of stress-relaxation behavior of poly (methyl methacrylate)(pmma) across the glass transition temperature. Mechanics of Materials 71, 74–84. Medvedev, G. A., Caruthers, J. M., 2013. Development of a stochastic constitutive model for prediction of postyield softening in glassy polymers. Journal of Rheology 57 (3), 949–1002. Medvedev, G. A., Caruthers, J. M., 2015. Stochastic model prediction of nonlinear creep in glassy polymers. Polymer 74, 235–253. Medvedev, G. A., Caruthers, J. M., 2016. A comparison of constitutive descriptions of the thermo-mechanical behavior of polymeric glasses. In: Polymer Glasses. CRC Press, pp. 467–552. Nguyen, T., Qi, H. J., Castro, F., Long, K. N., 2008. A thermoviscoelastic model for amorphous shape memory polymers: Incorporating structural and stress relaxation. J. Mech. Phys. Solids 56, 2792–2814. Nguyen, T. D., 2013. Modeling shape-memory behavior of polymers. Polymer Reviews 53 (1), 130–152. Nguyen, T. D., Yakacki, C. M., Brahmbhatt, P. D., Chambers, M. L., 2010. Modeling the Relaxation Mechanisms of Amorphous Shape Memory Polymers. Advanced Materials 22 (31), 3411–3423. Nieuwenhuizen, T. M., 1998. Thermodynamics of the glassy state: Effective temperatureas an additional system parameter. Phys. Rev. Lett. 80, 5580–5583. Qi, H. J., Nguyen, T. D., Castro, F., Yakacki, C., Shandas, R., 2008. Finite deformation thermo-mechanical behavior of thermally-induced shape memory polymers. J. Mech. Phys. Solids 56, 1730–1751. Reese, S., Govindjee, S., 1997. Theoretical and numerical aspects in the thermo-viscoelastic material behaviour of rubber-like polymers. Mechanics of Time-Dependent Materials 1 (4), 357–396. Reese, S., Govindjee, S., 1998. A theory of finite viscoelasticity and numerical aspects. Int. J. Solids Struct. 35, 3455–82. Santiago, D., Fabregat-Sanjuan, A., Ferrando, F., De la Flor, S., 2016. Recovery stress and work output in hyperbranched poly (ethyleneimine)-modified shape-memory epoxy polymers. Journal of Polymer Science Part B: Polymer Physics 54 (10), 1002–1013. Semkiv, M., Anderson, P. D., H¨utter, M., 2017. Two-scale model for the effect of physical aging in elastomers filled with hard nanoparticles. Journal of Computational Physics 350, 184–206. Semkiv, M., H¨utter, M., 2016. Modeling aging and mechanical rejuvenation of amorphous solids. Journal of NonEquilibrium Thermodynamics 41 (2), 79–88. Srivastava, V., Chester, S. A., Ames, N. M., Anand, L., 2010. A thermo-mechanically-coupled large-deformation theory for amorphous polymers in a temperature range which spans their glass transition. International Journal of Plasticity 26 (8), 1138–1182.

35

Su, X., Peng, X., 2018. A 3d finite strain viscoelastic constitutive model for thermally induced shape memory polymers based on energy decomposition. International Journal of Plasticity 110, 166–182. Tool, A. Q., 1946. Viscosity and extraordinary heat effects in glass. J. Am. Ceram. Soc. 29, 240–253. Tschoegl, N. W., 1989. The Phenomenological Theory of Linear Viscoelastic Material Behaviour. Springer Verlag. van Breemen, L., Klompen, E. T., Govaert, L., Meijer, H., 2011. Extending the egp constitutive model for polymer glasses to multiple relaxation times. J. Mech. Phys. Solids 59, 2191–2207. Xiao, R., Choi, J., Lakhera, N., Yakacki, C., Frick, C., Nguyen, T., 2013. Modeling the glass transition of amorphous networks for shape-memory behavior. Journal of the Mechanics and Physics of Solids 61 (7), 1612–1635. Xiao, R., Ghazaryan, G., Tervoort, T. A., Nguyen, T. D., 2017a. Modeling energy storage and structural evolution during finite viscoplastic deformation of glassy polymers. Physical Review E 95 (6), 063001. Xiao, R., Guo, J., Tian, C., Chen, W., 2017b. Modeling enthalpy relaxation using the mittag-leffler function. Journal of Non-Crystalline Solids 465, 17–25. Xiao, R., Nguyen, T. D., 2015. An effective temperature theory for the nonequilibrium behavior of amorphous polymers. Journal of the Mechanics and Physics of Solids 82, 62–81. Xiao, R., Sun, H., Chen, W., 2017c. A finite deformation fractional viscoplastic model for the glass transition behavior of amorphous polymers. International Journal of Non-Linear Mechanics 93, 7–14. Xiao, R., Tian, C., 2019. A constitutive model for strain hardening behavior of predeformed amorphous polymers: Incorporating dissipative dynamics of molecular orientation. Journal of the Mechanics and Physics of Solids 125, 472–487. Xiao, R., Yakacki, C. M., Guo, J., Frick, C. P., Nguyen, T. D., 2016. A predictive parameter for the shape memory behavior of thermoplastic polymers. Journal of Polymer Science Part B: Polymer Physics 54 (14), 1405–1414. Yu, C., Kang, G., Chen, K., 2017. A hygro-thermo-mechanical coupled cyclic constitutive model for polymers with considering glass transition. International Journal of Plasticity 89, 29–65. Yu, K., Xie, T., Leng, J., Ding, Y., Qi, H. J., 2012. Mechanisms of multi-shape memory effects and associated energy release in shape memory polymers. Soft Matter 8 (20), 5687–5695. Zhang, P., To, A. C., 2016. Transversely isotropic hyperelastic-viscoplastic model for glassy polymers with application to additive manufactured photopolymers. International Journal of Plasticity 80, 56–74. Zhang, R., Bai, P., Lei, D., Xiao, R., 2018. Aging-dependent strain localization in amorphous glassy polymers: From necking to shear banding. International Journal of Solids and Structures 146, 203–213.

36

Highlights • • • •

A viscoelastic-viscoplastic model is formulated based on the two-temperature thermodynamic framework. The model incorporates the molecular orientation with a relaxation mechanism for strain hardening. The simulation accurately captures the measured stress response at various temperatures and rates. The theory reproduces the main features of the constrained recovery performance of cold-programmed shape-memory polymers.

Declaration of interests

RThe authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: