Abstract

Chemically cross-linked thermoset shape memory polymers (TSMPs) are an important branch of smart materials due to their potentially wide applications in deplorable structures, soft robots, damage self-healing, and 4D printing. Further development and design of TSMP structures call for constitutive models. Although the Arruda–Boyce eight-chain model has been very successful and widely used for entropy-driven TSMPs, recent studies found that some new TSMPs, such as those using enthalpy as the primary driving force, show unit cells different from the eight-chain model. Considering that these new epoxy-based TSMP networks consist of a plenty of four-chain features, this study proposes a four-chain tetrahedron structure as the unit cell of the network to construct the constitutive model. In this model, Gibbs free energy is used to formulate the thermodynamic driving force. Then, by introducing a transition of the molecule deformation mechanism from that dominated by bond stretch to that dominated by bond angle opening, the traditional Langevin chain model is modified. It is found that this model can well capture the dramatic modulus change for the new TSMP in the thermomechanical experiments. Moreover, it shows that the original Treloar four-chain model and Arruda–Boyce eight-chain model underestimate the driving force for the enthalpy-driven TSMPs, and thus cannot well capture the thermomechanical behaviors. It is also found that under certain conditions, our four-chain model produces the same Cauchy stress as the eight-chain model does. This study may help researchers understand the thermomechanical response and design a special category of TSMPs with high recovery stress.

1 Introduction

Shape memory polymers (SMPs) have attracted growing attention for decades because they have many potential applications in engineering structures and devices [15]. SMPs can be broadly divided into two categories: one-way SMPs and two-way SMPs, including one-way multi-shape SMPs and two-way multi-shape SMPs. One-way SMPs are capable of deforming to a temporal shape and recovering to its original shape upon stimuli, i.e., one-way shape memory effect (SME) [613]; two-way SMPs are capable of switching reversibly between two different shapes upon stimuli, i.e., two-way SME [1418]. Several approaches have been used to trigger the shape memory effect by external stimuli such as heat [19], light [20], and moisture [21]. To date, SMPs have found wide applications in artificial muscles [22], self-deployable aerospace structures [23,24], soft robots [25], self-healing materials [26], fire dampers [27], 4D printing [28], and so on.

In order to design SMPs with good SME such as high shape recovery ratio and high recovery stress, a plenty of models have been developed to unravel the underlying mechanism. Basically, these models can be divided into several types: the storage strain-based macroscopic model [9,12,29,30], macroscopic rheological model [3134], viscoelasticity and viscoplasticity model [8,3537], molecular dynamics simulation [38], and multi-branch model [39] by combing statistical mechanics based Arruda–Boyce eight-chain model [40] and rheological model. Among them, the storage strain-based model, rheological model, and viscoelasticity and viscoplasticity model are macroscopic ones, which cannot provide a bridge between the compositional and topological features of the SMPs and their SMEs. In other words, they fail to explain the physical mechanisms in multiple length scales. The molecular dynamic model can provide in-depth understanding in the compositional and topological scale but is impractical for engineering design. As for the multi-branch model, it connects microscopic statistical mechanics to macroscopic deformation based on physical understanding and has gained wide recognitions [39,41].

In the Arruda–Boyce eight-chain model, a basic assumption is that the microscopic structures of cross-linked polymers possess diagonal square unit cells. However, for thermoset shape memory polymers (TSMPs), the unit cell varies from one TSMP to another. Therefore, it should be realized that the eight-chain model is only suitable for describing the thermomechanical response of a certain type of TSMPs, which are composed of molecular chains crosslinked to the diagonal square structure, instead of all types of TSMPs. For example, Fan and Li [42] synthesized a new TSMP made of commercially available epoxy resin (EPON 826, a bisphenol A based epoxy resin) and hardener (isophorone diamine (IPD)). Instead of the widely used entropy-driven mechanism, the driving force is primarily through enthalpy change, which leads to about ten times higher recovery stress in the rubbery state and in the bulk form. The enthalpy-driven mechanism is proved through characterizations by Raman spectroscopy and near-edge X-ray absorption fine structure spectroscopy. The characterizations show that the cross-linked network has high steric hindrance, leading to energy storage and release primarily through enthalpy change [42]. This polymer actually belongs to a class of TSMPs. In addition to the EPON 826 cured by IPD, EPON 826 cured by 1,3-bis(aminomethyl)cyclohexane (BACH) [42], and ultraviolet curing Bisphenol A glycerolate dimethacrylate (BPAGMA) [28], all demonstrated very high recovery stress at the rubbery state and in the bulk form. It is found that the unit cell of the new high recovery stress TSMP possesses a tetrahedron structure [42], see Fig. 1. This perceived unit cell was also echoed by molecular dynamic simulation for a similar EPON-IPD system (see Fig. 4 in Ref. [43]).

Fig. 1
(a) Planar schematic for molecular configuration of the EPON-IPD network; (b) unit cell in the EPON-IPD network; (c) monomer of IPD; (d) partial structure of the epoxy resin (EPON 826);  represents the extend structure. For each dotted circle in (a), it consists of one IPD at the center, connected by four epoxy arms, which are the EPON epoxy monomer. In 3D, this planar schematic is represented by the tetrahedron unit cell in (b). From (a), it is seen that every epoxy monomer arm is shared by two unit cells, growing into the 3D space network.
Fig. 1
(a) Planar schematic for molecular configuration of the EPON-IPD network; (b) unit cell in the EPON-IPD network; (c) monomer of IPD; (d) partial structure of the epoxy resin (EPON 826);  represents the extend structure. For each dotted circle in (a), it consists of one IPD at the center, connected by four epoxy arms, which are the EPON epoxy monomer. In 3D, this planar schematic is represented by the tetrahedron unit cell in (b). From (a), it is seen that every epoxy monomer arm is shared by two unit cells, growing into the 3D space network.
Close modal

Historically, a four-chain unit cell was firstly modeled by Treloar in 1946 for rubber elasticity [44] based on the assumption of the Flory–Rehner model which involved the calculation of the entropy formation for a network of molecules of equal contour or chain length. In this formulation, the driving force for rubbery elasticity was entropy. As shown later, this entropy-driven four-chain model cannot capture the thermomechanical behaviors of the enthalpy-driven TSMPs because of the underestimation of the driving force. The same problem exists for the Arruda–Boyce eight-chain model. Meanwhile, we also notice that enthalpy has been considered as one of the driving forces for the chemo-responsive shape recovery effect and has been used to constitute the free energy forms in a couple of studies [45,46]. Therefore, the formulation of the four-chain model must correctly represent the driving force, which is, in the current study, enthalpy. For this purpose, we present a Gibbs free energy driven constitutive model and compare it to the previous constitutive models with entropy as the driving forces, aiming at exploring and studying the difference under the different driving forces.

It is noted that there exists a longstanding controversy if the moduli parameters in the strain energy formulation are constants or not. So far, the moduli parameters in most models are treated as constants and are fitted by experiments (see Table 1; A10, A01, A30, A40, and A50 are constants). However, there are various molecule deformation mechanisms and energy barriers for polymers, and all of them could contribute to modulus variation. For example, according to Onck et al. [47], during the filamentous protein network deformation, there is a transition from a bending-dominated response at small deformation to a stretching-dominated response at large deformation, and thus different moduli are found corresponding to these two stages of deformation. Besides, according to Fan and Li [42], there is a transition from bond rotation in low energy state (at small deformation) to bond stretch in high energy state (at large deformation) for an EPON-IPD specimen, which also naturally causes a modulus change. Therefore, the modulus of TSMPs should be treated as a function of deformation, instead of a constant, which better represents the deformation mechanisms. Based on this understanding, the traditional Langevin chain model is modified by considering the bond angle change and bond stretch, which will be described in detail later.

Table 1

Different descriptions of strain energy of rubber elasticity

Model nameFormula
Neo-Hookean model [48]W = A10(I1 − 3)
Mooney–Rivlin model [49]W = A10(I1 − 3) + A01(I2 − 3)
Yeoh model [50]W=A10(I13)+A01(I23)2+A30(I13)3
Gent model [51]W=A402Jmlog(1I13Jm)
Ogden model [49]W=i=1NA50αi(λ1αi+λ2αi+λ3αi3)
Model nameFormula
Neo-Hookean model [48]W = A10(I1 − 3)
Mooney–Rivlin model [49]W = A10(I1 − 3) + A01(I2 − 3)
Yeoh model [50]W=A10(I13)+A01(I23)2+A30(I13)3
Gent model [51]W=A402Jmlog(1I13Jm)
Ogden model [49]W=i=1NA50αi(λ1αi+λ2αi+λ3αi3)

The aim of this study is to develop a three-dimensional model which is able to predict the thermomechanical behaviors of enthalpy-driven TSMPs with four-chain unit cells. Section 2 introduces the constitutive model on the basis of the basic thermodynamics laws. In Sec. 3, the constitutive model is verified by an EPON-IPD specimen under thermomechanical experiments. Next, a comparison between the classical Arruda–Boyce eight-chain model, the original Treloar four-chain model, and our new four-chain model is elaborated in Sec. 4. Finally, in Sec. 5, some important conclusions are drawn.

2 Constitutive Model

2.1 The Geometry of the Four-Chain Model.

As shown in Fig. 2(a), we set up the origin of the Cartesian coordinate system at the body center of the tetrahedron. In the tetrahedron model, let the original length between the origin and the four vertexes (A, B, C, and D) be r0. Four chains are assumed to originate from vertexes and crosslinked at the body center O. Also, the deformation of the four chains will not be affected by the other unit cells of the network. It should be noticed that this is also an average assumption (in reality, the origin can be at any location in the tetrahedron). Therefore, the coordinates for A, B, C, and D are
A=(0,0,1)r0,B=(0,223,13)r0,C=(23,23,13)r0,D=(23,23,13)r0
(1)
where r0 is the length of a single representative polymer chain (see Fig. 2(b)). The stretch in the three principal directions are λ1, λ2, and λ3, respectively. With a rigid body rotation Q, the principal stretch λ1, λ2, and λ3 in the new Cartesian coordinate system along the direction OX′, OY′, OZ′ are related to the original Cartesian coordinate system
x=QxwhereQ=[l1l2l3m1m2m3n1n2n3]
(2)
where l1, l2, l3, m1, m2, m3, n1,n2, and n3 are the direction cosine between the original coordinate and the principal stretch. The coordinates for A, B, C, and D in the new Cartesian coordinate can be obtained from Eq. (2) as
A=r0(l3,m3,n3)B=r0(223l213l3,223m213m3,223n213n3)C=r0(23l123l213l3,23m123m213m3,23n123n213n3)D=r0(23l123l213l3,23m123m213m3,23n123n213n3)
(3)
Fig. 2
(a) Schematic of the four-chain shape memory polymer model. Four chains originated from vertexes A, B, C, and D are cross-linked at the body center O. (b) Geometry of a single representative polymer chain (r0: original chain length of a representative chain).
Fig. 2
(a) Schematic of the four-chain shape memory polymer model. Four chains originated from vertexes A, B, C, and D are cross-linked at the body center O. (b) Geometry of a single representative polymer chain (r0: original chain length of a representative chain).
Close modal
The average stretch square of the four chains is
r¯2=OA¯2+OB¯2+OC¯2+OD¯24
(4)
where
OA¯2=λ12r02l32+λ22r02m32+λ32r02n32
(5)
OB¯2=λ12r02(223l213l3)2+λ22r02(223m213m3)2+λ32r02(223n213n3)2
(6)
OC¯2=λ12r02(23l123l213l3)2+λ22r02(23m123m213m3)2+λ32r02(23n123n213n3)2
(7)
OD¯2=λ12r02(23l123l213l3)2+λ22r02(23m123m213m3)2+λ32r02(23n123n213n3)2
(8)
The summation of the direction cosine in every direction of the coordinate is
l12+l22+l32=1,m12+m22+m32=1,n12+n22+n32=1
(9)
After some complex and tedious calculations by combining Eqs. (S30)–(S36) in Supplemental S3 available in the Supplemental Materials on the ASME Digital Collection, we obtain the square of the average chain length
r¯2=r02(λ12+λ22+λ32)3
(10)
Therefore, the average single chain stretch is
λ¯mic=r¯r0=13λ12+λ22+λ32=13I1(C)
(11)
Based on basic thermodynamics, the Cauchy stress for enthalpy-driven TSMP can be further written as (see the details in Supplemental Material S1)
σ=1J([2(ψI1(I13I1F1FT))F1J(σ:EF+E:σF)](pF))FT+κ(J1)F1FTFFT3κα(θθ0)I
(12)
Furthermore, by using the chain rule for the average single chain stretch, the Cauchy stress can be written as
σ=1J([2(ψλ¯micλ¯micI1)(I13I1F1FT)F1J(σ:EF+E:σF)](pF))FT+κ(J1)F1FTFFT3κα1(θθ0)I
(13)

Since the derivative λ¯mic/I1 can be obtained from Eq. (11), we still need to solve ψ/λ¯mic, otherwise the Cauchy stress cannot be derived. In order to obtain the explicit form of the Cauchy stress, the Helmholtz free energy form and the transition of the molecule deformation mechanism are derived in Sec. 2.2.

2.2 Helmholtz Free Energy Form and Molecule Deformation Mechanism Transition.

The infinitesimal probability dP for the conformation of a chain falls in the geometrical length interval [λ¯mic, λ¯mic+dλ¯mic] is [52]
dP(λ¯mic)=p(λ¯mic)dλ¯mic
(14)
The probability density for the stretch of a representative polymer chain is [53]
p(λ¯mic)=γexp{n(λ¯micnβ+lnβsinhβ)}
(15)
where rv is a normalization constant, n is the segment number of a single representative polymer chain, and the parameter β can be determined by the microscopic stretch of the Langevin function as
L(β)=cothβ1β=λ¯micn
(16)
Combining Eqs. (14)(16), the Helmholtz free energy density for the polymer with the random walk chain whose conformation incorporates the non-Gaussian statistics can be simply written as [54] follows:
ψ=μn(λ¯micnβ+lnβsinhβ)
(17)
where µ is the shear modulus of the material. In previous works [39,40,52], the material parameter µ is treated as a constant. From the microscopic level, a constant parameter µ means that the constraints on the molecule chain are constant, which suggests that only one type of microscopic deformation mechanism exists in the whole deformation process. However, as indicated above, some researchers [42,55,56] pointed out that there are different microscopic deformation mechanisms responsible for the macro-stretch and thus, correspondingly, the moduli should depend on deformation. According to Lyons [56], the modulus is determined by both the elongation of the repeating unit due to stretching of valence bonds and bond angle opening. Therefore, the material parameter µ is a deformation-dependent variable instead of a constant. As shown in Figs. 3(a) and 3(b), Young’s moduli for these two mechanisms (bond length change and bond angle opening) are defined, respectively, as [56]
Es=l0b0i=1n11kicos2θi,Ea=l0b0i=1n11disin2θi
(18)
where θi is the angle between the bond direction and the chain axis, which is assumed as a constant [56]; l0 and b0 are the repeating distance along the chain and the average cross-sectional area assignable to each chain molecule (or the area of the projected basal plane); ki and di are the interatomic force constant and angular force constant; n1 is the number of bond type in a segment.
Fig. 3
(a) Schematic definition of a zigzag shape of the valence and axial angle (αi is the valence angle.) (b) Schematic of bond stretch and bond angle opening upon loading by P.
Fig. 3
(a) Schematic definition of a zigzag shape of the valence and axial angle (αi is the valence angle.) (b) Schematic of bond stretch and bond angle opening upon loading by P.
Close modal
Because the energy barrier for bond angle opening is greater than that for bond stretching, we assume that the bond stretching occurs with a greater probability at the beginning or at small strain, and then bond angel opening gradually becomes the dominating mechanism with the increase in the average microscopic stretch. Based on this understanding, it is reasonable for us to define µ as
μ=12(1+ν)(fs(λ¯mic)Es+fa(λ¯mic)Ea)withfa(λ¯mic)=1fs(λ¯mic)
(19)
where fs and fa are the probabilities for the bond stretch and bond angle opening, respectively. fs is assumed to accord with quasi-Gaussian distribution, which has been commonly adopted in polymer science [29,57,58]
ps(λ¯mic)={p0exp[(λ¯micλ¯m)2/2Σ2],λ¯mic>00,λ¯mic<0
(20)
in which λ¯m is the mean of tube stretch, p0 is the normalization factor, and Σ is the standard deviation. Therefore, the fraction of possibility for modulus from chain stretch at certain microscopic stretch λ¯mic is
fs(λ¯mic)=λ¯micp(λ¯)dλ¯
(21)
The derivative of the strain energy ψ for representative polymer chains with respect to the average microscopic stretch reads
ψλ¯mic=μ(λ¯micnβ+lnβsinhβ)r¯r¯λ¯mic+μλ¯mic(λ¯micnβ+lnβsinhβ)=μnL1(λ¯micn)+μλ¯mic(λ¯micnβ+lnβsinhβ)
(22)
where the first term in Eq. (22) on the right-hand side is obtained from the freely jointed chain model [59], and the derivative of the shear modulus with respect to the microscopic stretch in the second term is
μλ¯mic=p0exp[(3λ¯mic3λ¯m)/18Σ2]18λ¯micb0n1sin2θi(3λ¯micλ¯mic3+2λ¯mic2)lodi+p0exp[(3λ¯mic3λ¯m)/18Σ2]18λ¯micb0n1cos2θi(3λ¯micλ¯mic3+2λ¯mic2)loki
(23)
Combining the expression of Cauchy stress (Eq. (13)) with average single chain stretch (Eq. (11)) and the derivative (Eq. (22)), the Cauchy stress can be further written as
σ=2MJ1(I13I1F1FT)FFT1J2(σ:EF+E:σF)FT1J(pF)FT+κ(J1)F1FTFFT3κα(θθ0)I
(24)
where I is the second-order identity tensor and M reads
M=μnL1(λ¯micn)(123I11/2)+(p0l018λ¯micb0n1(2λ¯mic2λ¯mic2)exp(3(λ¯micλ¯m)18Σ2)×(disin2θi+kicos2θi))(λ¯micnβ+lnβsinhβ)
(25)

Since the constitutive model for a polymer under mechanical loading has been derived, the constitutive model for TSMPs can be further developed.

2.3 Constitutive Model for Thermoset Shape Memory Polymers.

Similar to the phase model presented by Qi et al. [39], the total stress can be decomposed into two parts
σ=fgσg+frσr
(26)
where fg and fr are the volume fractions of the glassy phase and rubbery phase, respectively, which are defined as
fg=111+cf(ThT)n2
(27)
fr=1fg=11+cf(ThT)n2
(28)
where cf and n2 are the parameters for the phase evolution law, respectively. It is worth noting that the initial glassy phase is not considered in this constitutive model because as Qi et al. [39] indicated, the pre-deformation step is performed at θ >> θg and thus the fraction of initial glassy phase is about zero. σg and σr are the stress in the glassy phase and stress in the rubbery phase, respectively. By considering the TSMPs as an incompressible material (µr = 0.5), the stress of the rubbery phase is simplified as
σr=2[Mr(I13I1F1FT)]FFT1J(σrEF+EσrF)FT(prF)FT3κrαr(θθ0)I
(29)
where
Mr=μnrL1(λ¯micnr)(123I11/2)+(p0l018λ¯micb0n1(2λ2λ2)exp(3(λ¯micλ¯m)18Σ2)×(dirsin2θi+kircos2θi))(λ¯micnβ+lnβsinhβ)
(30)

Under relaxation, the stress is σrel = σr · exp(−t/τ), where τ is the relaxation time of the TSMP.

Generally speaking, the response to stress for glassy polymers can be divided into four steps: pure elastic deformation, yielding, strain softening and plastic flow, and strain hardening. The pure elastic deformation, yielding, and strain softening can be represented by the viscoplastic branch in Fig. 4, which was described in detail by Boyce et al. [60,61] and has been widely used later. According to Boyce et al. [60], the molecular chains in the glassy state are initially frozen and coordinated segmental rotation is prohibited, which corresponds to the elastic step. Once the stretch reaches a certain level, coordinated segmental rotation becomes possible, leading to yielding. After that, the polymer exhibits strain softening and plastic flow. Finally, the taut network is further stretched with an increase in stress, leading to strain hardening and ultimate breaking. Basically, the first three steps (elastic, yielding, and strain softening/plastic flow) can be represented by the viscoplastic branch. For the strain hardening, it can be represented by hyperelasticity. However, as discussed above, this deformation mechanism does not take effect until the stretch comes to a certain level. Therefore, in Fig. 4, we added another branch to represent the hyperelasticity. However, this branch needs to be switched off in the viscoplastic stage and needs to be turned on once the viscoplastic stage is over, see Fig. 4. Also, as summarized by van Melick et al. [62], in all the models for glassy polymers, including “BPA model” [63], “full chain model” [64], and “Neo-Hookean (Gaussian) relation model” [62], the strain hardening is modeled as a rubber-elastic model. Additionally, some recent studies also indicated that strain hardening can be modeled by the viscoplastic part [65]. Therefore, both hyperelastic and viscoplastic methods can be used to model strain hardening. In this study, we adopt the hyperelastic strain hardening model because it can be easily implemented, i.e., this method does not require us to introduce the back stress tensor to compute the yield surface shifting in every time-step, and thus saves computational resources. However, we do believe that the viscoplastic method may provide some new insights to our model, which deserves investigation in future studies. Motived by these observations and studies, the stress of the glassy phase σg is divided into two parts, i.e., the viscoplastic stress σgp and the switchable hyperelastic (rubbery) stress σgr (see Fig. 4)
σg=σgp+σgr
(31)
Fig. 4
One-dimensional rheological analogy for the TSMP at the glassy state (the switch is turned on once the critical equivalent strain is achieved)
Fig. 4
One-dimensional rheological analogy for the TSMP at the glassy state (the switch is turned on once the critical equivalent strain is achieved)
Close modal
The stress for the switchable hyperelastic part follows the rubbery elasticity as
σgr=H(EeqEc)(2[Mg(I13I1Fgr1FgrT)]FgrFgrT(σgEgrFgr+EσgrFgr)FgrT(pgFgr)FgrT3κgαg(θθ0)I)
(32)
where H(x) is the Heaviside function as the “switch” and can be written as
H(x)={0,x<01,x>0
(33)
and Fgr and Egr are defined as
Fgr=FFc,Egr=EEc
(34)
in which Fc and Ec are the deformation gradient and strain tensor corresponding to the critical equivalent strain Ec, respectively. Eeq is the equivalent strain and reads
Eeq=23E:E
(35)
It is worth noting that the critical equivalent strain Ec is dependent on temperature and strain rate. Mg in Eq. (32) reads
Mg=μngL1(λ¯micrng)(123(I1r)1/2)+(p0l018λ¯micb0n1(2λ¯micr2(λ¯micr)2)exp(3(λ¯micrλ¯mr)18Σ2)×(digsin2θi+kigcos2θi))(λ¯micrnβ+lnβsinhβ)
(36)
where
λ¯micr=λ¯micλ¯cmic
(37)
in which λ¯cmic is the microscopic stretch corresponding to the critical equivalent strain and I1r can be evaluated from the stretch increment when the critical equivalent strain occurs. It is assumed that the glassy phase and rubbery phase share the same energy transition parameters except for the force constants. By using the multiplicative decomposition scheme, the elastic deformation gradient for the glassy phase can be determined from the total deformation gradient and plastic gradient
Fe=FFp1
(38)
The viscoplastic stress is defined as [39]
σgp=1Je[Ce:Ee3αgκg(θθ0)I]
(39)
where Je is the Jacobin for the elastic part, Ce is the isotropic elasticity tensor, I is the second-order identity tensor, and κg is the bulk modulus for the glassy phase.
The spatial velocity gradient is
l=F˙F1=F˙e(Fe)1+Felp(Fe)1
(40)
Without losing generality, the spin rate of the spatial velocity gradient is assumed to be zero [66]. Then, the spatial velocity gradient of the plastic part can be represented by the stretch rate, which reads
lp=Dv=γ˙p2τ¯σg
(41)
where σgp is the deviatoric part of stress for the viscoplastic phase, and τ¯ is the effective equivalent shear stress and is defined as [60]
τ¯=[12σgp:σgp]1/2
(42)
γ˙p is the plastic shear strain rate and is defined as
γ˙p=γ˙0exp(ΔGkT(1(τ¯s)))
(43)
where γ˙0 is the pre-exponential factor proportional to the attempt frequency [66], ΔG is the zero stress level activation energy, and s represents the current athermal deformation resistance of the material, indicating the current state of the structure [60]. In order to describe the material softening evolution, the athermal shear rate is defined as [39,60,61]
s˙=h0(1s/ss)γ˙p
(44)
where ss is the saturation value of the athermal shear deformation resistance. The initial condition is
s=s0whenγp=0
(45)
in which s0 is the initial athermal shear deformation resistance. Besides, the deformation gradient for the temperature-dependent glassy phase does not inherit the deformation of the rubbery phase and will behave as an undeformed material. Hence, during cooling, the deformation increment of the glassy phase is given by a piecewise function [39]
ΔFgn+1={Fn+1(Fn)1,Δθ0I,Δθ=0
(46)
where Fn and Fn+1 are the overall deformation gradient at the instant tn and tn+1, and the total deformation gradient evolves as
Fgn+1=ΔFgn+1Fgn
(47)
in which, the Fg will be used to replace F in the glassy phase.
In the constrained stress-recovery process, it is assumed that the conversion efficiency from the released Gibbs free energy to the stored Gibbs free energy in the programming process is ηe, i.e., gre = ηegr, which will work as a driving force in the recovery process. Under the boundary condition F = F0, the Cauchy stress for the rubbery phase at the rubbery state can be written as
σr=(2ηe[Mr(I13I1F1FT)]FFT1Jηe(σr:EF+E:σrF):FT(pF)FT3ηeκgα(θθ0)I)|F=F0
(48)

3 Model Validation

In this section, by assuming the TSMP as an isotropic material, the model is validated by a TSMP (EPON-IPD) under the uniaxial compression test [42]. Specifically, the modeling for the TSMP in the glassy phase, rubbery phase and recovery stress under different constraint boundary conditions are conducted by using the computation process in Supplemental Material S8 and the material parameters in Table 2. Because we only compare the simulation with experimental data at 100% rubbery state (fg = 0) or 100% glassy state (fg = 1), the thermal expansion effect and phase transition parameters are not necessary and thus are not listed in Table 2. The parameter calibration and sensitivity analysis are presented in Supplemental Materials S5 and S6.

Table 2

Model parameters

DescriptionParameterValue
Common parameters for both rubbery and glassy phases
Type of bonds for the single segmentn1 (−)4
Average bond lengthl0 (Å)1.5
Average bond angle between bond direction and chain axisθi (deg)35
Cross-sectional area assigned to each chain moleculeb02)2.32
Exclusive parameters for the rubbery phase
Mean of average microscopic stretchλ¯rm (−)1.04
Standard deviation of microscopic stretchΣ (−)0.4
Interatomic force constantkir (dyn/cm)0.160
Angular force constantdir (dyn/cm)0.467
Conversion efficiency for Gibbs free energyη0.215
Segment number of the single representative polymer chainnr (−)1.5
Exclusive parameters for glassy phase
Mean of average microscopic stretchλ¯gm (−)1.01
Standard deviationΣ (−)0.3
Interatomic force constant for equilibrium time-independent behaviorkig (dyn/cm)0.604
Angular force constant for equilibrium time-independent behaviordig (dyn/cm)0.982
Elastic modulus for non-equilibrium time-dependent behaviorEg (MPa)1350
Zero level activation energyΔG (Pa)0.4 × 10−19
Pre-exponential factorγ˙0 (s−1)0.05
Initial value of athermal shear strengths0 (MPa)69
Saturation value of athermal shear strengthss63
Segment number of the single representative polymer chainng1.25
Parameter for material softeningh0 (MPa)2.405 × 103
Critical equivalent strainεc0.27
Other parameters
Relaxation timeτ (min)11.01
DescriptionParameterValue
Common parameters for both rubbery and glassy phases
Type of bonds for the single segmentn1 (−)4
Average bond lengthl0 (Å)1.5
Average bond angle between bond direction and chain axisθi (deg)35
Cross-sectional area assigned to each chain moleculeb02)2.32
Exclusive parameters for the rubbery phase
Mean of average microscopic stretchλ¯rm (−)1.04
Standard deviation of microscopic stretchΣ (−)0.4
Interatomic force constantkir (dyn/cm)0.160
Angular force constantdir (dyn/cm)0.467
Conversion efficiency for Gibbs free energyη0.215
Segment number of the single representative polymer chainnr (−)1.5
Exclusive parameters for glassy phase
Mean of average microscopic stretchλ¯gm (−)1.01
Standard deviationΣ (−)0.3
Interatomic force constant for equilibrium time-independent behaviorkig (dyn/cm)0.604
Angular force constant for equilibrium time-independent behaviordig (dyn/cm)0.982
Elastic modulus for non-equilibrium time-dependent behaviorEg (MPa)1350
Zero level activation energyΔG (Pa)0.4 × 10−19
Pre-exponential factorγ˙0 (s−1)0.05
Initial value of athermal shear strengths0 (MPa)69
Saturation value of athermal shear strengthss63
Segment number of the single representative polymer chainng1.25
Parameter for material softeningh0 (MPa)2.405 × 103
Critical equivalent strainεc0.27
Other parameters
Relaxation timeτ (min)11.01

It can be clearly seen from Figs. 57 that the model results agree with the experimental data well for most of the deformation. Specifically, for the rubbery state during stepwise programming (Fig. 6), the model can capture the dramatic modulus change in the whole deformation process; for the glassy state (Fig. 5), the model can well capture some critical characteristics, including the initial elastic deformation, yield, and strain softening as well as the post-yield strain hardening. However, there are some small discrepancies especially under small deformation in the stepwise relaxation process (Fig. 6). One reason may be that we only consider constant relaxation time in this study, but in reality, the relaxation time τ could be within a certain range due to the statistical nature of the polymer network. For the recovery stress-recovery strain curve in Fig. 7, the model also captures the test results. It is worth noting that the stresses in Fig. 7 were measured after the relaxation process under different constrained conditions, which accords with the general definition for stress relaxation. Therefore, calling it recovery stress after relaxation may be more accurate than calling it recovery stress. It is also noted that in addition to change the molecular structures such as the EPON-IPD network, the constrained recovery stress is also controlled by programming temperature and pre-strain [37], which is a useful strategy to make TSMPs with higher recovery stress.

Fig. 5
Comparison of Cauchy stress and Cauchy strain for the TSMP in the glassy state under uniaxial compression
Fig. 5
Comparison of Cauchy stress and Cauchy strain for the TSMP in the glassy state under uniaxial compression
Close modal
Fig. 6
Comparison of Cauchy stress and Cauchy strain for the TSMP in the rubbery state during the uniaxial compression programming process with stepwise relaxation (in each step, the specimen was compressed Δε = 2% and then relaxed while keep strain constant for 4 min at 170 °C)
Fig. 6
Comparison of Cauchy stress and Cauchy strain for the TSMP in the rubbery state during the uniaxial compression programming process with stepwise relaxation (in each step, the specimen was compressed Δε = 2% and then relaxed while keep strain constant for 4 min at 170 °C)
Close modal
Fig. 7
Comparison of recovery Cauchy stress and recovery Cauchy strain for the TSMP in the rubbery state
Fig. 7
Comparison of recovery Cauchy stress and recovery Cauchy strain for the TSMP in the rubbery state
Close modal

From the sensitivity analysis, the rankings for the sensitivity coefficient corresponding to +10% and +20% parameter variation can be obtained in Table S1 (see details in Supplemental Material S6). By comparing the absolute values of the sensitivity coefficient, we obtained the rankings Xss>Xεc>Xng>Xθi>Xλ¯m>Xdig>Xl0>Xn1>Xb0>XΣ>Xγ˙0>XΔG>XEg>Xkig>Xh0>Xs0 and Xss>Xεc>Xng>Xθi>Xλ¯m>Xdig>Xl0>Xn1>Xb0>XΣ>Xγ˙0>XΔG>XEg>Xkig>Xh0>Xs0 for +10% and +20% parameter variations, respectively.

4 Comparison of Constitutive Behavior Among the New Four-Chain Model, the Original Four-Chain Model, and the Eight-Chain Model

It seems that the Cauchy stress in Eq. (24) is relatively complex. However, if we ignore the pressure and microscopic deformation transition as well as adopt Helmholtz free energy as the driving force, and taking a closer look at the Cauchy stress expression in Eq. (24) again, we find that
σ=μ3J1nλ¯micL1(λ¯micn)(I13I1F1FT)FFT+κ(J1)F1FTFFT3κα1(θθ0)I
(49)
If the deformation is symmetric, then the Cauchy stress for TSMP at the rubbery state can be simplified as
σr=μr3nλ¯micL1(λ¯micn)(B13I1I)3κα1(θθ0)I
(50)
where B is the left Cauchy Green tensor. It is found that Eq. (50) is equal to the Cauchy stress of the Arruda–Boyce eight-chain model. In other words, the four-chain model can be reduced to the eight-chain model by ignoring the pressure term and the modulus transition mechanism as well as the difference between Helmholtz free energy and Gibbs free energy. It deserves mentioning that the eight-chain model is an affine model, involves a very small number of model parameters, and is much simpler in calculation. However, as indicated in previous studies, cross-linked polymer networks are inherently non-affine [6771]. Therefore, the eight-chain model is an ideal simplification, although it works pretty well for a large number of rubbery networks. By considering all of these points and the identification of a lot four-chain features in the molecular dynamics modeling for this type of polymer networks, we believe that the non-affine four-chain model in our study, which includes stiffer Langevin chains and considers the pressure term, the dependence of modulus on deformation, and the difference between the Helmholtz free energy and Gibbs free energy, is appropriate for this type of shape memory polymer network with high steric hindrance.
In addition, in the original four-chain model, the energy density is defined as [44]
W=12μ(I1(C)3)
(51)
where µ is the initial shear modulus. Hence, the corresponding Cauchy stress is
σof=2JFWCFT
(52)

Furthermore, by assuming the incompressibility for the TSMP in the rubbery state and using the segment number of single representative polymer chain n = 1.5 (the same as the parameter nr in Table 2), the compressive stresses for the TSMP at the rubbery state corresponding to the new four-chain model, previous Treloar four-chain model, and Arruda–Boyce eight-chain model can be compared in Fig. 8. It should be mentioned that the constant shear modulus for the Arruda–Boyce model is obtained by measuring the initial average Young’s modulus within the strain range −0.2 < ε < 0, which is about 12 MPa and is also used for shear modulus in the original Treloar four-chain model. According to the new four-chain model in this study, the microscopic constraint changes because of the transition between two microscopic deformation mechanisms, thus the stress shows different slopes under small deformation and under comparatively large deformation. Clearly, our new four-chain model well captures the experimental trend. On the other hand, because the original Treloar four-chain model and the classical Arruda–Boyce eight-chain model only consider constant shear modulus, the corresponding stress slope cannot dramatically change as the strain gradually increases (the segment number n does not affect the Arruda–Boyce eight-chain model very much when it changes from 1.5 to 100, see Fig. S1 in Supplemental Materials), which results in a considerable deviation from the experimental data. In other words, the driving forces for both the original Treloar four-chain model and Arruda–Boyce eight-chain model, which come from entropy, underestimate the real driving forces, thus the new enthalpy-driven four-chain model is an important supplement for predicting the thermomechanical characteristics of the TSMPs such as EPON-IPD.

Fig. 8
Comparison between the stresses for the driving forces of the new four-chain model, the original Treloar four-chain model, and Arruda–Boyce eight-chain model
Fig. 8
Comparison between the stresses for the driving forces of the new four-chain model, the original Treloar four-chain model, and Arruda–Boyce eight-chain model
Close modal

5 Conclusions

A three-dimensional four-chain constitutive model based on the basic thermodynamics framework is developed in this study, which is capable of capturing the thermomechanical behaviors for an enthalpy-driven EPON-IPD shape memory polymer programmed by uniaxial compression loading. Some important conclusions can be drawn as follows:

First, on the basis of the first law and second law of thermodynamics, the total driving force for the TSMP is represented by Gibbs free energy. By using this free energy, the thermomechanical behaviors of TSMPs can be better predicted, which is significantly underestimated by the entropy-driven Arruda–Boyce eight-chain model and the original Treloar four-chain model under finite deformation.

Second, the new four-chain model is formulated in this study by modifying the Langevin chain. Because this unit cell is similar to the microscopic molecular structure of EPON-IPD, it is more suitable for describing the deformation process for this type of high recovery stress TSMPs. Also, this model may be expanded to describe other types of TSMPs with similar unit cells.

Third, the unusual stiffness increase for the TSMP in the rubbery state under large deformation is due to the transition between two different microscopic deformation mechanisms with the increase in macroscopic deformation. In particular, the strain hardening behavior is modeled by a switchable hyperelastic branch based on the understanding of the physical deformation mechanisms of the polymer network.

Fourth, from the parameter sensitivity analysis in Supplemental Materials, the saturation value of the athermal shear deformation resistance and the critical equivalent strain are the most essential parameters to this model, while the evolution of the athermal shear deformation resistance and the initial value of athermal shear deformation resistance do not contribute too much to the model output.

Acknowledgment

The authors gratefully acknowledge the financial support by the National Science Foundation under grant number 1736136 and NASA cooperative agreement NNX16AQ93A under contract number NASA/LEQSF(2016-19)-Phase3-10.

References

1.
Leng
,
J.
,
Lan
,
X.
,
Liu
,
Y.
, and
Du
,
S.
,
2011
, “
Shape-Memory Polymers and Their Composites: Stimulus Methods and Applications
,”
Prog. Mater. Sci.
,
56
(
7
), pp.
1077
1135
. 10.1016/j.pmatsci.2011.03.001
2.
Hager
,
M. D.
,
Bode
,
S.
,
Weber
,
C.
, and
Schubert
,
U. S.
,
2015
, “
Shape Memory Polymers: Past, Present and Future Developments
,”
Prog. Polym. Sci.
,
49–50
, pp.
3
33
. 10.1016/j.progpolymsci.2015.04.002
3.
Hu
,
J.
,
Zhu
,
Y.
,
Huang
,
H.
, and
Lu
,
J.
,
2012
, “
Recent Advances in Shape-Memory Polymers: Structure, Mechanism, Functionality, Modeling and Applications
,”
Prog. Polym. Sci.
,
37
(
12
), pp.
1720
1763
. 10.1016/j.progpolymsci.2012.06.001
4.
Liu
,
C.
,
Qin
,
H.
, and
Mather
,
P. T.
,
2007
, “
Review of Progress in Shape-Memory Polymers
,”
J. Mater. Chem.
,
17
(
16
), pp.
1543
1558
. 10.1039/b615954k
5.
Meng
,
H.
, and
Li
,
G.
,
2013
, “
A Review of Stimuli-Responsive Shape Memory Polymer Composites
,”
Polymer
,
54
(
9
), pp.
2199
2221
. 10.1016/j.polymer.2013.02.023
6.
Lendlein
,
A.
, and
Kelch
,
S.
,
2002
, “
Shape Memory Polymers
,”
Angew. Chem. Int. Ed.
,
41
(
12
), pp.
2034
2057
. 10.1002/1521-3773(20020617)41:12<2034::aid-anie2034>3.0.co;2-m
7.
Li
,
G.
, and
Xu
,
T.
,
2011
, “
Thermomechanical Characterization of Shape Memory Polymer-Based Self-Healing Syntactic Foam Sealant for Expansion Joints
,”
J. Transp. Eng.
,
137
(
11
), pp.
805
814
. 10.1061/(ASCE)TE.1943-5436.0000279
8.
Li
,
G.
, and
Xu
,
W.
,
2011
, “
Thermomechanical Behavior of Thermoset Shape Memory Polymer Programmed by Cold-Compression: Testing and Constitutive Modeling
,”
J. Mech. Phys. Solids
,
59
(
6
), pp.
1231
1250
. 10.1016/j.jmps.2011.03.001
9.
Liu
,
Y.
,
Gall
,
K.
,
Dunn
,
M. L.
,
Greenberg
,
A. R.
, and
Diani
,
J.
,
2006
, “
Thermomechanics of Shape Memory Polymers: Uniaxial Experiments and Constitutive Modeling
,”
Int. J. Plast.
,
22
(
2
), pp.
279
313
. 10.1016/j.ijplas.2005.03.004
10.
Xie
,
T.
,
2010
, “
Tunable Polymer Multi-Shape Memory Effect
,”
Nature
,
464
(
7286
), pp.
267
270
. 10.1038/nature08863
11.
Xu
,
W.
, and
Li
,
G.
,
2011
, “
Thermoviscoplastic Modeling and Testing of Shape Memory Polymer Based Self-Healing Syntactic Foam Programmed at Glassy Temperature
,”
ASME J. Appl. Mech.
,
78
(
6
), p.
061017
. 10.1115/1.4004554
12.
Yan
,
C.
, and
Li
,
G.
,
2019
, “
Design Oriented Constitutive Modeling of Amorphous Shape Memory Polymers and Its Application to Multiple Length Scale Lattice Structures
,”
Smart Mater. Struct.
,
28
(
9
), p.
095030
. 10.1088/1361-665X/ab230c
13.
Yu
,
K.
,
Ge
,
Q.
, and
Qi
,
H. J.
,
2014
, “
Reduced Time as a Unified Parameter Determining Fixity and Free Recovery of Shape Memory Polymers
,”
Nat. Commun.
,
5
(
1
), p.
3066
. 10.1038/ncomms4066
14.
Chung
,
T.
,
Romo-Uribe
,
A.
, and
Mather
,
P. T.
,
2008
, “
Two-Way Reversible Shape Memory in a Semicrystalline Network
,”
Macromolecules
,
41
(
1
), pp.
184
192
. 10.1021/ma071517z
15.
Lu
,
L.
,
Cao
,
J.
, and
Li
,
G.
,
2018
, “
Giant Reversible Elongation Upon Cooling and Contraction Upon Heating for a Crosslinked Cis Poly(1,4-Butadiene) System at Temperatures Below Zero Celsius
,”
Sci. Rep.
,
8
(
1
), p.
14233
. 10.1038/s41598-018-32436-9
16.
Meng
,
H.
, and
Li
,
G.
,
2013
, “
Reversible Switching Transitions of Stimuli-Responsive Shape Changing Polymers
,”
J. Mater. Chem. A
,
1
(
27
), pp.
7838
7865
. 10.1039/c3ta10716g
17.
Zhao
,
Q.
,
Qi
,
H. J.
, and
Xie
,
T.
,
2015
, “
Recent Progress in Shape Memory Polymer: New Behavior, Enabling Materials, and Mechanistic Understanding
,”
Prog. Polym. Sci.
,
49–50
, pp.
79
120
. 10.1016/j.progpolymsci.2015.04.001
18.
Lu
,
L.
, and
Li
,
G.
,
2016
, “
One-Way Multishape-Memory Effect and Tunable Two-Way Shape Memory Effect of Ionomer Poly(Ethylene-Co-Methacrylic Acid)
,”
ACS Appl. Mater. Interfaces
,
8
(
23
), pp.
14812
14823
. 10.1021/acsami.6b04105
19.
Lendlein
,
A.
,
2002
, “
Biodegradable, Elastic Shape-Memory Polymers for Potential Biomedical Applications
,”
Science
,
296
(
5573
), pp.
1673
1677
. 10.1126/science.1066102
20.
Lendlein
,
A.
,
Jiang
,
H.
,
Jünger
,
O.
, and
Langer
,
R.
,
2005
, “
Light-Induced Shape-Memory Polymers
,”
Nature
,
434
(
7035
), pp.
879
882
. 10.1038/nature03496
21.
Ma
,
M.
,
Guo
,
L.
,
Anderson
,
D. G.
, and
Langer
,
R.
,
2013
, “
Bio-Inspired Polymer Composite Actuator and Generator Driven by Water Gradients
,”
Science
,
339
(
6116
), pp.
186
189
. 10.1126/science.1230262
22.
Yang
,
Q.
,
Fan
,
J.
, and
Li
,
G.
,
2016
, “
Artificial Muscles Made of Chiral Two-Way Shape Memory Polymer Fibers
,”
Appl. Phys. Lett.
,
109
(
18
), p.
183701
. 10.1063/1.4966231
23.
Santo
,
L.
,
Quadrini
,
F.
,
Accettura
,
A.
, and
Villadei
,
W.
,
2014
, “
Shape Memory Composites for Self-Deployable Structures in Aerospace Applications
,”
Procedia Eng.
,
88
, pp.
42
47
. 10.1016/j.proeng.2014.11.124
24.
Liu
,
Y.
,
Du
,
H.
,
Liu
,
L.
, and
Leng
,
J.
,
2014
, “
Shape Memory Polymers and Their Composites in Aerospace Applications: A Review
,”
Smart Mater. Struct.
,
23
(
2
), p.
23001
. 10.1088/0964-1726/23/2/023001
25.
Jin
,
B.
,
Song
,
H.
,
Jiang
,
R.
,
Song
,
J.
,
Zhao
,
Q.
, and
Xie
,
T.
,
2018
, “
Programming a Crystalline Shape Memory Polymer Network With Thermo- and Photo-Reversible Bonds Toward a Single-Component Soft Robot
,”
Sci. Adv.
,
4
(
1
), p.
3865
. 10.1126/sciadv.aao3865
26.
Li
,
G.
,
2014
,
Self-Healing Composites: Shape Memory Polymer Based Structures
,
John Wiley & Sons Ltd
,
West Sussex, UK
.
27.
Feng
,
X.
,
Fan
,
J.
,
Li
,
A.
, and
Li
,
G.
,
2019
, “
Multireusable Thermoset With Anomalous Flame-Triggered Shape Memory Effect
,”
ACS Appl. Mater. Interfaces
,
11
(
17
), pp.
16075
16086
. 10.1021/acsami.9b03092
28.
Li
,
A.
,
Challapalli
,
A.
, and
Li
,
G.
,
2019
, “
4D Printing of Recyclable Lightweight Architectures Using High Recovery Stress Shape Memory Polymer
,”
Sci. Rep.
,
9
(
1
), p.
7621
. 10.1038/s41598-019-44110-9
29.
Yang
,
Q.
, and
Li
,
G.
,
2016
, “
Temperature and Rate Dependent Thermomechanical Modeling of Shape Memory Polymers With Physics Based Phase Evolution Law
,”
Int. J. Plast.
,
80
, pp.
168
186
. 10.1016/j.ijplas.2015.09.005
30.
Baghani
,
M.
,
Arghavani
,
J.
, and
Naghdabadi
,
R.
,
2014
, “
A Finite Deformation Constitutive Model for Shape Memory Polymers Based on Hencky Strain
,”
Mech. Mater.
,
73
, pp.
1
10
. 10.1016/j.mechmat.2013.11.011
31.
Yan
,
C.
,
Yang
,
Q.
, and
Li
,
G.
,
2020
, “
A Phenomenological Constitutive Model for Semicrystalline Two-Way Shape Memory Polymers
,”
Int. J. Mech. Sci.
, 10.1016/j.ijmecsci.2020.105552
32.
Tobushi
,
H.
,
Hashimoto
,
T.
,
Hayashi
,
S.
, and
Yamada
,
E.
,
1997
, “
Thermomechanical Constitutive Modeling in Shape Memory Polymer of Polyurethane Series
,”
J. Intell. Mater. Syst. Struct.
,
8
(
8
), pp.
711
718
. 10.1177/1045389X9700800808
33.
Tobushi
,
H.
,
Okumura
,
K.
,
Hayashi
,
S.
, and
Ito
,
N.
,
2001
, “
Thermomechanical Constitutive Model of Shape Memory Polymer
,”
Mech. Mater.
,
33
(
10
), pp.
545
554
. 10.1016/S0167-6636(01)00075-8
34.
Balogun
,
O.
, and
Mo
,
C.
,
2014
, “
Shape Memory Polymers: Three-Dimensional Isotropic Modeling
,”
Smart Mater. Struct.
,
23
(
4
), p.
045008
. 10.1088/0964-1726/23/4/045008
35.
Nguyen
,
T. D.
,
Qi
,
H. J.
,
Castro
,
F.
, and
Long
,
K. N.
,
2008
, “
A Thermoviscoelastic Model for Amorphous Shape Memory Polymers: Incorporating Structural and Stress Relaxation
,”
J. Mech. Phys. Solids
,
56
(
9
), pp.
2792
2814
. 10.1016/j.jmps.2008.04.007
36.
Li
,
G.
, and
Shojaei
,
A.
,
2012
, “
A Viscoplastic Theory of Shape Memory Polymer Fibres With Application to Self-Healing Materials
,”
Proc. R. Soc. A Math. Phys. Eng. Sci.
,
468
(
2144
), pp.
2319
2346
. 10.1098/rspa.2011.0628
37.
Dai
,
L.
,
Tian
,
C.
, and
Xiao
,
R.
,
2020
, “
Modeling the Thermo-Mechanical Behavior and Constrained Recovery Performance of Cold-Programmed Amorphous Shape-Memory Polymers
,”
Int. J. Plast.
,
127
, p.
102654
. 10.1016/j.ijplas.2019.102654
38.
Diani
,
J.
, and
Gall
,
K.
,
2007
, “
Molecular Dynamics Simulations of the Shape-Memory Behaviour of Polyisoprene
,”
Smart Mater. Struct.
,
16
(
5
), p.
1575
1583
. 10.1088/0964-1726/16/5/011
39.
Qi
,
H. J.
,
Nguyen
,
T. D.
,
Castro
,
F.
,
Yakacki
,
C. M.
, and
Shandas
,
R.
,
2008
, “
Finite Deformation Thermo-Mechanical Behavior of Thermally Induced Shape Memory Polymers
,”
J. Mech. Phys. Solids
,
56
(
5
), pp.
1730
1751
. 10.1016/j.jmps.2007.12.002
40.
Aruda
,
E. M.
, and
Boyce
,
M. C.
,
1993
, “
A Three-Dimensional Constitutive Model for the Large Stretch Behavior of Rubber Elastic Materials
,”
J. Mech. Phys. Solids
,
41
(
2
), pp.
389
412
. 10.1016/0022-5096(93)90013-6
41.
Westbrook
,
K. K.
,
Kao
,
P. H.
,
Castro
,
F.
,
Ding
,
Y.
, and
Jerry Qi
,
H.
,
2011
, “
A 3D Finite Deformation Constitutive Model for Amorphous Shape Memory Polymers: A Multi-Branch Modeling Approach for Nonequilibrium Relaxation Processes
,”
Mech. Mater.
,
43
(
12
), pp.
853
869
. 10.1016/j.mechmat.2011.09.004
42.
Fan
,
J.
, and
Li
,
G.
,
2018
, “
High Enthalpy Storage Thermoset Network With Giant Stress and Energy Output in Rubbery State
,”
Nat. Commun.
,
9
(
1
), p.
642
. 10.1038/s41467-018-03094-2
43.
Wu
,
C.
, and
Xu
,
W.
,
2006
, “
Atomistic Molecular Modelling of Crosslinked Epoxy Resin
,”
Polymer
,
47
(
16
), pp.
6004
6009
. 10.1016/j.polymer.2006.06.025
44.
Treloar
,
L. R. G.
,
1946
, “
The Elasticity of a Network of Longchain Molecules—III
,”
Trans. Faraday Soc.
,
42
(0), pp.
83
94
. 10.1039/TF9464200083
45.
Lu
,
H.
, and
Du
,
S.
,
2014
, “
A Phenomenological Thermodynamic Model for the Chemo-Responsive Shape Memory Effect in Polymers Based on Flory-Huggins Solution Theory
,”
Polym. Chem.
,
5
(
4
), pp.
1155
1162
. 10.1039/C3PY01256E
46.
Lu
,
H.
,
Liu
,
Y.
,
Leng
,
J.
, and
Du
,
S.
,
2010
, “
Qualitative Separation of the Physical Swelling Effect on the Recovery Behavior of Shape Memory Polymer
,”
Eur. Polym. J.
,
46
(
9
), pp.
1908
1914
. 10.1016/j.eurpolymj.2010.06.013
47.
Onck
,
P. R.
,
Koeman
,
T.
,
Van Dillen
,
T.
, and
Van Der Giessen
,
E.
,
2005
, “
Alternative Explanation of Stiffening in Cross-Linked Semiflexible Networks
,”
Phys. Rev. Lett.
,
95
(
17
), pp.
19
22
. 10.1103/physrevlett.95.178102
48.
Ogden
,
R. W.
,
1997
,
Non-Linear Eleastic Deformation
,
Dover Publication
,
New York
.
49.
Kim
,
N.-H.
,
2015
,
Introduction to Nonlinear Finite Element Analysis
,
Springer US
,
New York
.
50.
Yeoh
,
O. H.
,
1993
, “
Some Forms of the Strain Energy Function for Rubber
,”
Rubber Chem. Technol.
,
66
(
5
), pp.
751
771
. 10.5254/1.3538343
51.
Pucci
,
E.
, and
Saccomandi
,
G.
,
2011
, “
A Note on the Gent Model for Rubber-Like Materials
,”
Rubber Chem. Technol.
,
75
(
5
), pp.
839
852
. 10.5254/1.3547687
52.
Kroon
,
M.
,
2011
, “
An 8-Chain Model for Rubber-Like Materials Accounting for Non-Affine Chain Deformations and Topological Constraints
,”
J. Elast.
,
102
(
2
), pp.
99
116
. 10.1007/s10659-010-9264-7
53.
Kuhn
,
W.
, and
Grün
,
F.
,
1942
, “
Beziehungen Zwischen Elastischen Konstanten Und Dehnungsdoppelbrechung Hochelastischer Stoffe
,”
Kolloid-Z.
,
101
(
3
), pp.
248
271
. 10.1007/BF01793684
54.
Hossain
,
M.
,
Amin
,
A. F. M. S.
, and
Kabir
,
M. N.
,
2015
, “
Eight-Chain and Full-Network Models and Their Modified Versions for Rubber Hyperelasticity: A Comparative Study
,”
J. Mech. Behav. Mater.
,
24
(
1–2
), pp.
11
24
. 10.1515/jmbm-2015-0002
55.
Bowden
,
P. B.
,
1968
, “
The Elastic Modulus of an Amorphous Glassy Polymer
,”
Polymer
,
9
(
C
), pp.
449
454
. 10.1016/0032-3861(68)90054-2
56.
Lyons
,
W. J.
,
1958
, “
Theoretical Values of the Dynamic Stretch Moduli of Fiber-Forming Polymers
,”
J. Appl. Phys.
,
29
(
10
), pp.
1429
1433
. 10.1063/1.1722962
57.
Shojaei
,
A.
,
Sharafi
,
S.
, and
Li
,
G.
,
2015
, “
A Multiscale Theory of Self-Crack-Healing With Solid Healing Agent Assisted by Shape Memory Effect
,”
Mech. Mater.
,
81
, pp.
25
40
. 10.1016/j.mechmat.2014.10.008
58.
Shojaei
,
A.
, and
Li
,
G.
,
2014
, “
Thermomechanical Constitutive Modelling of Shape Memory Polymer Including Continuum Functional and Mechanical Damage Effects
,”
Proc. R. Soc. A Math. Phys. Eng. Sci.
,
470
(
2170
), p.
20140199
. 10.1098/rspa.2014.0199
59.
Bergström
,
J.
,
2015
,
Mechanics of Solid Polymers: Theory and Computational Modelling
,
Elsevier
,
New York
.
60.
Boyce
,
M. C.
,
Parks
,
D. M.
, and
Argon
,
A. S.
,
1988
, “
Large Inelastic Deformation of Glassy Polymers. Part I: Rate Dependent Constitutive Model
,”
Mech. Mater.
,
7
(
1
), pp.
15
33
. 10.1016/0167-6636(88)90003-8
61.
Boyce
,
M. C.
,
Parks
,
D. M.
, and
Argon
,
A. S.
,
1988
, “
Large Inelastic Deformation of Glassy Polymers. Part II: Numerical Simulation of Hydrostatic Extrusion
,”
Mech. Mater.
,
7
(
1
), pp.
35
47
. 10.1016/0167-6636(88)90004-X
62.
van Melick
,
H. G. H.
,
Govaert
,
L. E.
, and
Meijer
,
H. E. H.
,
2003
, “
On the Origin of Strain Hardening in Glassy Polymers
,”
Polymer
,
44
(
8
), pp.
2493
2502
. 10.1016/S0032-3861(03)00112-5
63.
Boyce
,
M. C.
, and
Arruda
,
E. M.
,
1990
, “
An Experimental and Analytical Investigation of the Large Strain Compressive and Tensile Response of Glassy Polymers
,”
Polym. Eng. Sci.
,
30
(
20
), pp.
1288
1298
. 10.1002/pen.760302005
64.
Wu
,
P. D.
, and
Van Der Giessen
,
E.
,
1993
, “
On Improved Network Models for Rubber Elasticity and Their Applications to Orientation Hardening in Glassy Polymers
,”
J. Mech. Phys. Solids
,
41
(
3
), pp.
427
456
. 10.1016/0022-5096(93)90043-F
65.
Xiao
,
R.
, and
Tian
,
C.
,
2019
, “
A Constitutive Model for Strain Hardening Behavior of Predeformed Amorphous Polymers: Incorporating Dissipative Dynamics of Molecular Orientation
,”
J. Mech. Phys. Solids
,
125
, pp.
472
487
. 10.1016/j.jmps.2019.01.008
66.
Boyce
,
M.
,
Kear
,
K.
,
Socrate
,
S.
, and
Shaw
,
K.
,
2001
, “
Deformation of Thermoplastic Vulcanizates
,”
J. Mech. Phys. Solids
,
49
(
5
), pp.
1073
1098
. 10.1016/S0022-5096(00)00066-1
67.
Miehe
,
C.
,
Göktepe
,
S.
, and
Lulei
,
F.
,
2004
, “
A Micro-Macro Approach to Rubber-Like Materials—Part I: The Non-Affine Micro-Sphere Model of Rubber Elasticity
,”
J. Mech. Phys. Solids
,
52
(
11
), pp.
2617
2660
. 10.1016/j.jmps.2004.03.011
68.
Cioroianu
,
A. R.
,
Spiesz
,
E. M.
, and
Storm
,
C.
,
2016
, “
Disorder, Pre-Stress and Non-Affinity in Polymer 8-Chain Models
,”
J. Mech. Phys. Solids
,
89
, pp.
110
125
. 10.1016/j.jmps.2016.01.014
69.
Ehret
,
A. E.
,
2015
, “
On a Molecular Statistical Basis for Ogden’s Model of Rubber Elasticity
,”
J. Mech. Phys. Solids
,
78
, pp.
249
268
. 10.1016/j.jmps.2015.02.006
70.
Basu
,
A.
,
Wen
,
Q.
,
Mao
,
X.
,
Lubensky
,
T. C.
,
Janmey
,
P. A.
, and
Yodh
,
A. G.
,
2011
, “
Nonaffine Displacements in Flexible Polymer Networks
,”
Macromolecules
,
44
(
6
), pp.
1671
1679
. 10.1021/ma1026803
71.
Drozdov
,
A. D.
, and
Gottlieb
,
M.
,
2005
, “
Constitutive Equations for Non-Affine Polymer Networks With Slippage of Chains
,”
Continuum Mech. Thermodyn.
,
17
(
3
), pp.
217
246
. 10.1007/s00161-004-0200-6

Supplementary data