## 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 [1–5]. 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) [6–13]; two-way SMPs are capable of switching reversibly between two different shapes upon stimuli, i.e., two-way SME [14–18]. 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 [31–34], viscoelasticity and viscoplasticity model [8,35–37], 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]).

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; *A*_{10}, *A*_{01}, *A*_{30}, *A*_{40}, and *A*_{50} 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.

Model name | Formula |
---|---|

Neo-Hookean model [48] | W = A_{10}(I_{1} − 3) |

Mooney–Rivlin model [49] | W = A_{10}(I_{1} − 3) + A_{01}(I_{2} − 3) |

Yeoh model [50] | $W=A10(I1\u22123)+A01(I2\u22123)2+A30(I1\u22123)3$ |

Gent model [51] | $W=\u2212A402Jmlog(1\u2212I1\u22123Jm)$ |

Ogden model [49] | $W=\u2211i=1NA50\alpha i(\lambda 1\alpha i+\lambda 2\alpha i+\lambda 3\alpha i\u22123)$ |

Model name | Formula |
---|---|

Neo-Hookean model [48] | W = A_{10}(I_{1} − 3) |

Mooney–Rivlin model [49] | W = A_{10}(I_{1} − 3) + A_{01}(I_{2} − 3) |

Yeoh model [50] | $W=A10(I1\u22123)+A01(I2\u22123)2+A30(I1\u22123)3$ |

Gent model [51] | $W=\u2212A402Jmlog(1\u2212I1\u22123Jm)$ |

Ogden model [49] | $W=\u2211i=1NA50\alpha i(\lambda 1\alpha i+\lambda 2\alpha i+\lambda 3\alpha i\u22123)$ |

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.

*A*,

*B*,

*C*, and

*D*) be

*r*

_{0}. 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

*r*

_{0}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

*l*

_{1},

*l*

_{2},

*l*

_{3},

*m*

_{1}

*, m*

_{2}

*, m*

_{3},

*n*

_{1,}

*n*

_{2}, and

*n*

_{3}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

Since the derivative $\u2202\lambda \xafmic/\u2202I1$ can be obtained from Eq. (11), we still need to solve $\u2202\psi /\u2202\lambda \xafmic$, 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.

*dP*for the conformation of a chain falls in the geometrical length interval [$\lambda \xafmic$, $\lambda \xafmic+d\lambda \xafmic$] is [52]

*r*is a normalization constant,

_{v}*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

*µ*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]

*θ*is the angle between the bond direction and the chain axis, which is assumed as a constant [56];

_{i}*l*

_{0}and

*b*

_{0}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);

*k*and

_{i}*d*are the interatomic force constant and angular force constant;

_{i}*n*

_{1}is the number of bond type in a segment.

*µ*as

*f*and

_{s}*f*are the probabilities for the bond stretch and bond angle opening, respectively.

_{a}*f*is assumed to accord with quasi-Gaussian distribution, which has been commonly adopted in polymer science [29,57,58]

_{s}*p*

_{0}is the normalization factor, and Σ is the standard deviation. Therefore, the fraction of possibility for modulus from chain stretch at certain microscopic stretch $\lambda \xafmic$ is

*ψ*for representative polymer chains with respect to the average microscopic stretch reads

**I**is the second-order identity tensor and M reads

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.

*f*and

_{g}*f*are the volume fractions of the glassy phase and rubbery phase, respectively, which are defined as

_{r}*c*and

_{f}*n*

_{2}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

*θ*>>

*θ*and thus the fraction of initial glassy phase is about zero.

_{g}

*σ**and*

_{g}

*σ**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}Under relaxation, the stress is *σ _{rel}* =

*σ*· exp(−

_{r}*t*/

*τ*), where

*τ*is the relaxation time of the TSMP.

**σ**

*is divided into two parts, i.e., the viscoplastic stress $\sigma gp$ and the switchable hyperelastic (rubbery) stress $\sigma gr$ (see Fig. 4)*

_{g}*H*(

*x*) is the Heaviside function as the “switch” and can be written as

**F**

_{c}and

**E**

_{c}are the deformation gradient and strain tensor corresponding to the critical equivalent strain

*E*, respectively.

_{c}*E*is the equivalent strain and reads

_{eq}*E*is dependent on temperature and strain rate.

_{c}*M*in Eq. (32) reads

_{g}*J*is the

^{e}*Jacobin*for the elastic part,

**C**

*is the isotropic elasticity tensor,*

^{e}**I**is the second-order identity tensor, and

*κ*is the bulk modulus for the glassy phase.

_{g}*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*is the saturation value of the athermal shear deformation resistance. The initial condition is

_{s}*s*

_{0}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]

*F**and*

^{n}

*F*

^{n}^{+1}are the overall deformation gradient at the instant

*t*and

_{n}*t*

_{n+}_{1}, and the total deformation gradient evolves as

*F**will be used to replace*

_{g}**in the glassy phase.**

*F**η*

_{e}, i.e.,

*g*=

_{re}*η*

_{e}

*g*

_{r,}which will work as a driving force in the recovery process. Under the boundary condition

**F**=

**F**

_{0}, the Cauchy stress for the rubbery phase at the rubbery state can be written as

## 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 (*f _{g}* = 0) or 100% glassy state (

*f*= 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.

_{g}Description | Parameter | Value |
---|---|---|

Common parameters for both rubbery and glassy phases | ||

Type of bonds for the single segment | n_{1} (−) | 4 |

Average bond length | l_{0} (Å) | 1.5 |

Average bond angle between bond direction and chain axis | θ (deg)_{i} | 35 |

Cross-sectional area assigned to each chain molecule | b_{0} (Å^{2}) | 2.32 |

Exclusive parameters for the rubbery phase | ||

Mean of average microscopic stretch | $\lambda \xafrm$ (−) | 1.04 |

Standard deviation of microscopic stretch | Σ (−) | 0.4 |

Interatomic force constant | k (dyn/cm)_{ir} | 0.160 |

Angular force constant | d (dyn/cm)_{ir} | 0.467 |

Conversion efficiency for Gibbs free energy | η | 0.215 |

Segment number of the single representative polymer chain | n (−)_{r} | 1.5 |

Exclusive parameters for glassy phase | ||

Mean of average microscopic stretch | $\lambda \xafgm$ (−) | 1.01 |

Standard deviation | Σ (−) | 0.3 |

Interatomic force constant for equilibrium time-independent behavior | k (dyn/cm)_{ig} | 0.604 |

Angular force constant for equilibrium time-independent behavior | d (dyn/cm)_{ig} | 0.982 |

Elastic modulus for non-equilibrium time-dependent behavior | E (MPa)_{g} | 1350 |

Zero level activation energy | ΔG (Pa) | 0.4 × 10^{−19} |

Pre-exponential factor | $\gamma \u02d90$ (s^{−1}) | 0.05 |

Initial value of athermal shear strength | s_{0} (MPa) | 69 |

Saturation value of athermal shear strength | s_{s} | 63 |

Segment number of the single representative polymer chain | n_{g} | 1.25 |

Parameter for material softening | h_{0} (MPa) | 2.405 × 10^{3} |

Critical equivalent strain | ε_{c} | 0.27 |

Other parameters | ||

Relaxation time | τ (min) | 11.01 |

Description | Parameter | Value |
---|---|---|

Common parameters for both rubbery and glassy phases | ||

Type of bonds for the single segment | n_{1} (−) | 4 |

Average bond length | l_{0} (Å) | 1.5 |

Average bond angle between bond direction and chain axis | θ (deg)_{i} | 35 |

Cross-sectional area assigned to each chain molecule | b_{0} (Å^{2}) | 2.32 |

Exclusive parameters for the rubbery phase | ||

Mean of average microscopic stretch | $\lambda \xafrm$ (−) | 1.04 |

Standard deviation of microscopic stretch | Σ (−) | 0.4 |

Interatomic force constant | k (dyn/cm)_{ir} | 0.160 |

Angular force constant | d (dyn/cm)_{ir} | 0.467 |

Conversion efficiency for Gibbs free energy | η | 0.215 |

Segment number of the single representative polymer chain | n (−)_{r} | 1.5 |

Exclusive parameters for glassy phase | ||

Mean of average microscopic stretch | $\lambda \xafgm$ (−) | 1.01 |

Standard deviation | Σ (−) | 0.3 |

Interatomic force constant for equilibrium time-independent behavior | k (dyn/cm)_{ig} | 0.604 |

Angular force constant for equilibrium time-independent behavior | d (dyn/cm)_{ig} | 0.982 |

Elastic modulus for non-equilibrium time-dependent behavior | E (MPa)_{g} | 1350 |

Zero level activation energy | ΔG (Pa) | 0.4 × 10^{−19} |

Pre-exponential factor | $\gamma \u02d90$ (s^{−1}) | 0.05 |

Initial value of athermal shear strength | s_{0} (MPa) | 69 |

Saturation value of athermal shear strength | s_{s} | 63 |

Segment number of the single representative polymer chain | n_{g} | 1.25 |

Parameter for material softening | h_{0} (MPa) | 2.405 × 10^{3} |

Critical equivalent strain | ε_{c} | 0.27 |

Other parameters | ||

Relaxation time | τ (min) | 11.01 |

It can be clearly seen from Figs. 5–7 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.

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\epsilon c>Xng>X\theta i>X\lambda \xafm>Xdig>Xl0$$>Xn1>Xb0>X\Sigma >X\gamma \u02d90>X\Delta G>XEg>Xkig>Xh0>Xs0$ and $Xss>X\epsilon c>Xng>X\theta i>X\lambda \xafm>Xdig>Xl0>Xn1$$>Xb0>X\Sigma >X\gamma \u02d90>X\Delta 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

**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 [67–71]. 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.

*µ*is the initial shear modulus. Hence, the corresponding Cauchy stress is

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 *n _{r}* 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.

## 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.