Multi-cycle deformation of supramolecular elastomers: Constitutive modeling and structure-property relations

International Journal of Engineering Science 133 (2018) 311–335

Contents lists available at ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

Multi-cycle deformation of supramolecular elastomers: Constitutive modeling and structure-property relations A.D. Drozdov∗, J. deClaville Christiansen Department of Materials and Production, Aalborg University Fibigerstraede 16, Aalborg 9220, Denmark

a r t i c l e

i n f o

Article history: Received 27 June 2018 Revised 24 September 2018 Accepted 8 October 2018

Keywords: Supramolecular elastomer Viscoelasticity Cyclic viscoplasticity Mullins effect

a b s t r a c t Supramolecular elastomers are rubber-like polymers with double-network structure formed by chains bridged by covalent and non-covalent bonds. Due to high mobility of temporary bonds, whose rate of recombination is comparable with the strain rate under loading, these materials demonstrate such properties as self-healing and self-recovery at ambient temperature. A constitutive model is developed for the viscoelastic and viscoplastic behavior of supramolecular elastomers. Stress–strain relations and the kinetic equations for plastic deformation are derived from the free energy imbalance inequality for an isothermal threedimensional deformation. The viscoelastic response reflects breakage and reformation of temporary junctions in an inhomogeneous transient network (transition of chains from their active to dangling state and vice versa). The viscoplastic response reflects slippage of permanent junctions with respect to their reference positions. A junction starts to slide when it becomes unbalanced due to transformation of one of the chains connected by this junction from its active state into the dangling state. The sliding process (plastic flow) proceeds until the junction reaches its new equilibrium state. The model is applied to fit experimental data in tensile relaxation tests, loadingunloading tests, and multi-cycle tests on supramolecular elastomers and triblock copolymers. Numerical simulation shows that the governing equations describe adequately the experimental stress–strain diagrams, the material parameters evolve consistently with chemical composition and experimental conditions, and predictions of the model are in qualitative agreement with observations. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Covalently cross-linked elastomers (vulcanized natural and synthetic rubbers) are widely used as seals, tires, shock absorbers and biomedical materials owing to their unique high elasticity. A shortcoming of these materials is that they are incapable of being healed after damage or injuries due to accidental cuttings or scratching during the service lifetime. Development of elastomers demonstrating fast repair at room temperature (without heating of the damaged area) has recently became a focus of attention (Chen, Kushner, Williams, & Guan, 2012; Cordier, Tournilhac, Soulie-Ziakovic, & Leibler, 2008; Kim et al., 2018; Li et al., 2016; Liu et al., 2017). The intrinsic self-healing ability of these materials allows not only to enhance their reliability, but opens new areas of their application as smart anticorrosion coatings (Xu, Ye, Ding, Tan, & Fu,

∗

Corresponding author. E-mail address: [email protected] (A.D. Drozdov).

https://doi.org/10.1016/j.ijengsci.2018.10.002 0020-7225/© 2018 Elsevier Ltd. All rights reserved.

312

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

2018; Zhang et al., 2017), biomimetic prostheses (Kang et al., 2018; Lei & Wu, 2018), energy harvesters (Parida et al., 2017) and substrates for stretchable electronic devices (Luo, Wan, Yang, Shah, & Chen, 2017; Yan et al., 2018). Autonomous self-healing of elastomers is driven by formation of a supramolecular architecture (Herbst, Dohler, Michael, & Binder, 2013): introduction of non-covalent junctions with finite lifetimes and reversible (dynamic) covalent bonds into their polymer networks. The healing process is associated with reconstruction (breakage and reformation) of non-covalent bonds at room temperature (Neal, Mozhdehi, & Guan, 2015) and recombination of reversible covalent junctions at elevated temperatures (Wu, Cai, & Weitz, 2017). The physical mechanisms and the kinetics of self-healing and selfrecovery were recently analyzed in Sheridan and Bowman (2013), Stukalin, Cai, Kumar, Leibler, and Rubinstein (2013) and Zhang, Chen, and Colby (2018). A constitutive framework for modeling the self-healing process was developed in Voyiadjis, Shojaei, and Li (2011), Hui and Long (2012), Wang, Gao, and Yu (2017), Kumar, Francfort, and Lopez-Pamies (2018) and Vernerey (2018), to mention a few. Design of novel supramolecular elastomers (polymer networks with chains bridged by permanent chemical cross-links and temporary non-covalent and dynamic junctions) was discussed in De Espinosa, Fiore, Weder, Johan Foster, and Simon (2015), Kajita, Noro, and Matsushita (2017) and Zhang, Rong, and Zhang (2018). Formation of non-covalent bonds in these materials is grounded on (i) electrostatic interaction, (ii) metal-ligand coordination, (iii) hydrophobic association, (iv) hydrogen bonding, and (v) host-guest recognition. The following necessary conditions for self-repair of elastomers at room temperature are formulated (Stukalin et al., 2013): (i) high concentration of temporary junctions, (ii) low activation energy for their rearrangement, and (iii) extraordinary mobility of polymer chains. Commercial diblock and triblock thermoplastic elastomers (TPEs) provide examples of physically cross-linked networks (where non-covalent bonds between polymer chains are formed by hard micro-domains arising due to phase separation) that demonstrate no self-healing and a limited self-recovery (Drozdov, 2007). This may be explained by the fact that soft domains in their phase-separated micro-structure (whose characteristic size does not exceed a fee nanometers; Li et al., 2016; Tomita et al., 2017) are severely constrained by hard domains. The latter results in a strong reduction of mobility of chains in the rubbery domains (evaluated by means of the characteristic relaxation time that exceeds 103 s; Cho et al., 2017; Tomita et al., 2017). Weakening of these constrains leads to a pronounced acceleration of the relaxation process (its characteristic time is reduced to 102 s; Tang, Huang, Guo, Zhang, & Liu, 2016; Zhao et al., 2012) accompanied by substantial enhancement of the self-healing ability (Chen & Guan, 2014; Kim et al., 2018; Yoshida, Ejima, & Yoshie, 2017). This study deals with constitutive modeling of the response of supramolecular elastomers in uniaxial tensile relaxation tests and multi-cycle loading–unloading tests with various deformation programs. A characteristic feature of the timedependent behavior of these materials is a strong and rapid relaxation of stresses in tests with large (several hundred percent) elongation ratios, where tensile stress decays by 30–70 per cent within a couple of min (Liu, Liu, et al., 2017; Tang et al., 2016; Wu, Qiu, Tang, Liu, & Guo, 2017; Xu et al., 2018). This implies that the relaxation rate becomes comparable with the strain rate, and the entire relaxation spectrum should be taken into account to describe the stress–strain diagrams (Grindy et al., 2015). Another interesting property of the response of supramolecular elastomers is observed in multi-cycle tests. Both supramolecular elastomers and conventional TPEs demonstrate the Mullins effect (Diani, Fayolle, & Gilormini, 2009): • Under deformation with a fixed maximum elongation ratio kmax , the loading and unloading paths differ substantially. • In tests with a fixed kmax , the hysteresis energy (estimated as the area between the loading and unloading curves) along the first cycle exceeds strongly that along the second cycle and decreases with number of cycles n. • Under deformation with monotonically increasing elongation ratios kmax n , the stress–strain diagram under reloading at the (n + 1 )th cycle coincides with that under uniaxial tension, provided that elongation ratio k exceeds kmax n . The qualitative difference between conventional TPEs and supramolecular elastomers consists in the shape of stress– strain diagrams under reloading. For TPEs, these curves are concave (Buckley, Prisacariu, & Martin, 2010; Cho, Rinaldi, & Boyce, 2013; Wang & Chester, 2018), which is also typical of vulcanized rubbers (Clough, Creton, Craig, & Sijbesma, 2016; Plagge & Kluppel, 2017). For supramoleculer elastomers, the reloading paths are S-shaped: these diagrams are convex at small stresses and concave at larger stresses (Beniah, Fortman, Heath, Dichtel, & Torkelson, 2017; Higaki et al., 2017; Jiang, Fang, Zhang, Wang, & Wang, 2017; Zhao et al., 2012). The objective of this work is threefold: (i) to develop a constitutive model for the mechanical behavior of supramolecular elastomers that describes the above features of their response, (ii) to determine adjustable parameters in the governing equations by fitting experimental data on several elastomers and to establish correlation between their mechanical properties and molecular architecture, and (iii) to examine the ability of the model to predict observations in multi-cycle tests with various deformation programs and arbitrary elongation ratios at which the strain rate changes its sign. Cyclic deformation of vulcanized (filled and unfilled) rubbers and conventional TPEs has recently been modeled in Cho et al. (2013), Dargazany and Itskov (2013), Guo, Guo, Amirkhizi, Zou, and Yuan (2016), Li, Tang, Kroger, and Liu (2016), Shen et al. (2016), Bacca, Creton, and McMeeking (2017), Cho et al. (2017), Plagge and Kluppel (2017), Lu, Wang, Yang, and Wang (2017), Pan and Zhong (2017), Mossi Idrissa, Wang, Ahzi, Patlazhan, and Remond (2018), Wang, Shan, Liao, and Xia (2018) and Zhou, Jiang, and Khayat (2018), to mention a few. Unlike the above studies that deal mainly with the description of the first cycle of loading–retraction, we focus on modeling the mechanical response in multi-cycle tests with stress- and strain-controlled deformation programs.

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

313

An equivalent polymer network in a supramolecular elastomer is modeled as a combination of permanent and transient networks. Dissipation of energy in the equivalent network (estimated as the hysteresis area between the loading and unloading paths under cyclic deformation) is induced by two processes: (i) breakage and reformation of physical bonds in the transient network driven by thermal fluctuations (viscoelasticity), and (ii) slippage of permanent cross-links with respect to their reference positions (viscoplasticity). The viscoelastic mechanism of energy dissipation is conventional for transient networks, where a dangling chain adopts the actual state at the instant of its reattachment as the reference (stress-free) state (Green & Tobolsky, 1946; Tanaka & Edwards, 1992). To describe observations in relaxation tests (which require the entire relaxation spectrum to be accounted for), we adopt the approach suggested in Drozdov and Gupta (2003) and Drozdov and Yuan (2003) and treat the equivalent network as heterogeneous and composed of meso-domains with various activation energies for separation of chains from temporary junctions. The viscoplastic mechanism of energy dissipation in non-covalent networks was recently proposed in Filippidi et al. (2017), see also Yanagisawa, Nan, Okuro, and Aida (2018). According to it, a junction between chains becomes unbalanced when one of the chains connected by this junction is transformed from the active state into the dangling state (which means that stress in this chain vanishes suddenly). As a result, the junction begins to slide with respect to the network (plastic flow) until it reaches a new equilibrium state. The novelty of our aproach consists in the description of both phenomena (viscoelasticity and viscoplasticity) within a unified constitutive model grounded on the free energy imbalance inequality. The exposition is organized as follows. Governing equations for the mechanical response of a supramolecular elastomer under uniaxial deformation are discussed in Section 2. Section 3 focuses on the analysis of observations on several elastomers and triblock copolymers in tensile relaxation tests, loading–retraction tests, and multi-cycle tests. Concluding remarks are formulated in Section 4. Derivation of the constitutive relations for an arbitrary three-dimensional deformation is provided in Appendix. 2. Model A supramolecular elastomer is modeled as an equivalent network of polymer chains. The network is treated as a superposition of two sub-networks: permanent and transient. Chains in the permanent network are connected by covalent (permanent) junctions, while chains in the transient network are bridged by non-covalent (temporary) bonds whose rearrangement (breakage and reformation) is driven by thermal fluctuations. According to the affinity hypothesis, deformations of the permanent and transient networks coincide with macro-deformation of the elastomer. The reference (stress-free) state of the permanent network before application of external loads coincides with its initial (undeformed) state. According to the multiplicative decomposition formula, the deformation gradient F for transition from the initial state into the actual state reads

F = Fe · Fp ,

(1)

where Fe and Fp are the deformation gradients for elastic and plastic deformations, and the dot stands for inner product. To describe the response of the transient network, we denote by τ an instant when an active chain (both ends are attached to different junctions) is bridged with the network, and distinguish chains that joined the network under preparation (τ = 0) and those connected with the network under deformation (τ > 0). The reference (stress-free) state of an active chain with τ = 0 coincides with the initial (undeformed) state. The reference state of an active chain with τ > 0 coincides with the actual state of the elastomer at instant τ (stresses in dangling chains relax entirely). The transient network is presumed to be composed of meso-domains with various activation energies U for separation of chains from temporary junctions. The inhomogeneity is characterized by the distribution function f(u) of meso-domains with dimensionless activation energy u = U/(kB T ), where kB is the Boltzmann constant, and T stands for the absolute temperature. The quasi-Gaussian expression is adopted for this function



f (u ) = f0 exp −

u2 2 2



( u ≥ 0 ),

(2)

where  > 0 is a material constant, and the coefficient f0 is determined from the normalization condition



∞ 0

f ( u )d u = 1 .

(3)

Rearrangement of the transient network is determined by the rate of detachment of active chains from temporary junctions

 = γ exp(−u ),

(4)

where γ > 0 stands for the attempt rate. Constitutive equations for the mechanical response of a supramolecular elastomer under three-dimensional deformation are developed in Appendix. These relations are derived by means of the free energy imbalance inequality for an arbitrary specific mechanical energy of the permanent network

We (Ie1 , Ie2 ) + Wp (Ip1 , Ip2 )

314

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

and an arbitrary strain energy of an active chain in the transient network

w(Iτ 1 , Iτ 2 ). The specific energy We stored in chains of the permanent network depends on the principal invariants Iem (m = 1, 2) of the Cauchy–Green tensor for elastic deformation Be = Fe · F e , where  denotes transpose. The specific energy of inter-chain interaction Wp is treated as a function of the principal invariants Ipm (m = 1, 2) of the Cauchy–Green tensor for plastic deformation bp = Fp · F p . The energy w stored in a chain of the transient network is presumed to depend on the principal −1 invariants Iτ m (m = 1, 2) of the Cauchy–Green tensor bτ = fτ · f τ , where fτ (t ) = F (t ) · F (τ ) is the deformation gradient for transition from the actual state at time τ to the actual state at time t. In the analysis of observations, we apply the Gent equation for the strain energy density of the permanent network and the neo–Hookean formulas for the specific energy of inter-chain interaction and the strain energy of active chains in the transient network





1 Ie1 − 3 We = − Ge K ln 1 − , 2 K

Wp =

1 Gp (Ip1 − 3 ), 2

w=

1 g(Iτ 1 − 3 ) 2

( τ ≥ 0 ),

(5)

where Ge , Gp , g and K are constants. The Gent equation (Gent, 1996) is adopted for the function We to account for strainhardening under tension with large elongation ratios. The neoHookean expression for the energy of inter-chain interaction serves as the first term in the formal expansion of the function Wp into the Taylor series with respect to its arguments. Governing equations for a supramolecular elastomer subjected to uniaxial tensile cyclic deformation with a constant strain rate ˙ involve (i) The formula for engineering tensile stress

   S2 k3 − 1 σ = G (1 − κ )V e + κ S1 k − 2 kke

(6)

k

with V given by Eq. (A-62), (ii) The kinematic equation for elongation ratio k under macro-deformation

k˙ = ±˙ ,

k ( 0 ) = 1,

(7)

where ˙ denotes the strain rate, the superscript dot stands for the derivative with respect to time t, and the signs “+” and “−” correspond to loading and unloading, respectively, (iii) The kinematic equation for elongation ratio ke for elastic deformation

k˙ p k˙ e k˙ = − , ke k kp

ke ( 0 ) = 1,

(8)

(iv) The kinetic equation for elongation ratio kp under plastic deformation

 k3 − 1 k3p − 1  k˙ p =P V e −R , kp ke kp

kp ( 0 ) = 1,

(9)

(v) The integral relations for the functions S1 and S2 (these functions characterize tensile stress in the transient network)



S1 =

∞ 0

f (u )s1 (t, u )du,



S2 =

0

∞

f (u )s2 (t, u )du,

(10)

(vi) The kinetic equations for the functions s1 and s2

1  ∂ s1 =  2 − s1 , ∂t k

∂ s2 = (k − s2 ), ∂t

s1 ( 0, u ) = s2 ( 0, u ) = 1.

(11)

Eqs. (2), (4) and (6)–(11) involve four material constants: G is the elastic modulus, K −1 is the Gent constant, κ is the ratio of elastic moduli of transient and permanent networks, and  is a measure of inhomogeneity of the transient network. The governing equations involve also three adjustable functions: γ , P and R. Parameter γ stands for the attempt rate for separation of active chains from temporary junctions. To account for slowing down of stress relaxation under retraction and its acceleration under reloading (Lu et al., 2017), we presume γ to be affected by plastic deformation and set

γ = γ∗ exp[−α (Ip2 − 3 )],

(12)

where γ ∗ is the attempt rate at infinitesimal strains, α is a constant, and Ip2 = 2kp + k−2 p denotes the second principal invariant of the Cauchy–Green tensor bp . Another approach to describe nonlinearity of the viscoelastic response was suggested in Drozdov and Kalamkarov (1996). In the analysis of multi-cycle deformation, three regimes are distinguished: (i) stretching of a virgin sample (k˙ > 0), (ii) retraction (k˙ < 0) down to the zero stress, and (iii) reloading (k˙ > 0) up to some maximum stress per cycle σ ∗ . We suppose

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

315

that the coefficient P vanishes under stretching and stress relaxation (k˙ = 0) and accepts different values, P1 and P2 , under retraction and reloading:

P=0

(stretching ),

P = P1

(retraction ),

P = P2

(reloading ).

(13)

The rate of plastic flow under retraction is given by

P1 = P1∗ exp[β1 (σ∗ − σ )],

(14)

where P1∗ is the rate of sliding of junctions at the instant when unloading starts and tensile stress reaches its local maximum σ ∗ , and the constant β 1 describes slowing down of plastic deformation at unloading driven by a decrease in tensile stress. To describe plastic deformation in elastomers with S-shaped reloading paths (convex at stresses σ < σ¯ and concave at stresses σ > σ¯ , where σ¯ is a material constant), the following equations are introduced:

P2 = P2∗ exp(−β2 σ )

(σ ≤ σ¯ ),

P2 = P¯2 [1 + β (σ − σ¯ )]

(σ > σ¯ ).

(15)

Here P2∗ is the rate of sliding of junctions at the instant when retraction starts and tensile stress vanishes, the constant β 2 describes slowing down of the sliding process at σ ≤ σ¯ , the coefficient P¯2 is the rate of sliding at σ = σ¯ where the shape of the stress–strain diagram changes, and the constant β characterizes stress-induced acceleration of plastic deformation at σ > σ¯ . By analogy with Eq. (13), the coefficient R is presumed to accept different values under retraction and reloading,

R = R1

(retraction ),

R = R2

(reloading ).

(16)

The coefficients P1∗ , P2∗ , R1 , R2 , remain constant along each retraction and reloading path, but they are affected by elongation ratios kmax and kmin at which the strain rate changes its sign. Evolution of these parameters with kmax and kmin reflects damage accumulation (changes in the micro-structure of an elastomer) under multi-cycle loading. To reduce the number of phenomenological relations, we treat R2 as a material constant (independent of kmin ). The influence of kmax on parameters P1∗ , R1 is described by the equations

log P1∗ = P10∗ − P11∗



Ie2 max − 3,

log R1 = R01 − R11



Ie2 max − 3,

(17)

where log = log10 , the coefficients are determined by the least-squares method, and the second principal invariant of the Cauchy–Green tensor for elastic deformation Ie2 = 2ke + k−2 is calculated at the point where the elongation ratio reaches e its local maximum, k = kmax . Eqs. (17) imply that the quantities P1∗ and R1 (parameters characterizing plastic deformation under retraction) are governed by the elastic deformation stored in the permanent network at the instant when the strain rate ˙ becomes negative. Two approaches can be mentioned to the description of the effect of kmin on the rate of plastic deformation under reloading P2∗ . According to the first, the coefficientP2∗ is affected by plastic deformation accumulated under retraction. Adopting an analog of Eq. (17) to describe this effect, we write

log P2∗ = P20∗ − P21∗



Ip2 min − 3,

(18)

where the coefficients are calculated by the least-squares technique, and Ip2 min stands for the second principal invariant of the Cauchy–Green tensor for plastic deformation Ip2 = 2kp + k−2 p at the point with k = kmin . According to the other approach, P2∗ is governed by the elastic deformation stored in the permanent network at the instant when reloading starts (this deformation differs from zero at the point with σmin = 0 due to viscoelastic effects). Within this concept, the dependence of P2∗ on kmin is describe by the relation similar to Eq. (17),

log P2∗ = P20∗ − P21∗



Ie2 min − 3,

(19)

where Ie2 min stands for the second principal invariant of the Cauchy–Green tensor for elastic deformation Ie2 at the point with k = kmin . Applicability of Eqs. (18) and (19) to fitting observations in multi-cycle tests with various programs will be discussed in Section 3. An important advantage of the model, Eqs. (2)–(19), is its ability to describe experimental stress–strain diagrams on supramolecular elastomers under multi-cycle deformations with arbitrary programs (elongation ratios kmax and kmin at which the strain rate changes its sign). To reach this goal, mutual dependencies are introduced between stresses and plastic flow (Eqs. (14) and (15)), between the relaxation rate and plastic deformation (Eq. (12)), and between the rates of plastic flow and elastic and plastic deformations (Eqs. (17)–(19)). As a result, the number of material constants in the governing equations increases. This number remains, however, comparable with the number of adjustable parameters in other models (mentioned in the introductory section) for the viscoelastic and viscoplastic response of elastomers under cyclic deformation. 3. Fitting of observations The ability of the model is examined to describe experimental stress–strain diagrams on supramolecular elastomers under loadings with various deformation programs.

316

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

Fig. 1. A: Tensile stress σ versus elongation ratio k. Symbols: experimental data in tensile tests with strain rate ˙ = 0.083 s−1 on IR-CND composites with various concentrations r phr (Wu, Qiu, et al., 2017). Solid lines: results of simulation. B: Ratio of stresses S versus relaxation time trel . Symbols: experimental data in relaxation tests with strain rate under stretching ˙ = 0.083 s−1 and elongation ratio k = 2. Solid lines: results of simulation. C: Parameters G and κ versus concentration r of CND. Circles: treatment of observations. Solid lines: approximation of the data by Eq. (20). D: Tensile stress σ versus elongation ratio k. Circles: experimental data in cyclic test with strain rate ˙ = 0.083 s−1 and kmax = 2.5 on IR-CND composite with r = 3 phr. Solid line: results of simulation.

3.1. Comparison of the response in relaxation and cyclic tests We demonstrate that the governing equations ensure good agreement with observations on three series of elastomers (with varied concentrations of physical bonds) in relaxation tests at finite strains and cyclic (loading–unloading) tests.

3.1.1. IR-CND nanocomposite elastomers We begin with matching observations on maleic anhydride (MAH) grafted polyisoprene rubber (IR) reinforced with amine-passivated carbon nanodots (CND) (Wu, Qiu, et al., 2017). Preparation of the composites involves three stages. First, isoprene rubber was melt-blended with MAH to produce MAH-grafted IR (4 mol.% of MAH). Then, amine-passivated CNDs were manufactured by microwave heating of carbon black (CB) particles in a phosphate solution. Finally, MAH-grafted IR was mixed with various amounts of CNDs (whose concentration r varied from 0 to 5 phr), sulfur-cured (1.5 phr of sulfur), and hot-pressed (143 °C) to prepare samples for mechanical tests. The permanent network in IR-CND nanocomposite elastomers is formed by IR chains (i) bridged by covalent (sulfur) cross-links and (ii) grafted on the surfaces of CNDs under thermal treatment. Physical junctions in the transient network are formed by hydrogen bonds between amino groups on the surfaces of CNDs and anhydrid moeties grafted on IR chains. Mechanical tests were conducted at room temperature on IR-CND nanocomposites with r = 0, 1, 2, 3 and 5 phr. The experimental program involves: (i) tensile tests with strain rate ˙ = 0.083 s−1 , (ii) relaxation tests with elongation ratio k = 2 and the same strain rate under stretching, and (iii) cyclic (loading–unloading) test with maximum elongation ratio kmax = 2.5 and the strain rate ˙ on IR-CND nanocomposite with r = 3. Observations in these tests are depicted in Fig. 1 together with results of simulation with the material constants listed in Table 1. In Fig. 1A and D, tensile stress σ is plotted versus elongation ratio k. In Fig. 1B, the dimensionless ratio of stresses S = σ (trel )/σ (trel = 0 ) is depicted versus relaxation time trel = t − t0 , where t0 is the instant when the required elongation ratio is reached.

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

317

Table 1 Material parameters for IR-CND composites (Fig. 1). Parameter

Value

K

14.7 0.5 4.0 0.0 0.0

γ ∗ (s−1 )  α β1

Table 2 Material parameters for IR-ATA-Zn elastomers (Fig. 2). Parameter

Value

K

110.0 0.03 2.6 0.0 0.0

γ ∗ (s−1 )  α β1

As relaxation tests are conducted at elongation ratio k = 2, while the characteristic times for stretching and stress relaxation are comparable, the standard procedure of fitting relaxation curves is inapplicable. Adjustable parameters are determined by means of the following iterative algorithm. We start with matching experimental data on the nanocomposite with r = 5 phr, fix γ ∗ , κ and  , and determine G and K from the best-fit condition for the stress–strain diagram under stretching in Fig. 1A. Then, for the fixed pair (G, K), we find γ ∗ , κ and  that ensure the best fit of the corresponding relaxation curve in Fig. 1B. Taking these parameters as a new approximation for the triplet (γ ∗ , κ ,  ), we return to approximation of the stress–strain curve in Fig. 1A. Simulation shows that three iterations are sufficient to ensure good agreement between the data in Fig. 1A and B and the results of numerical analysis. Observations on nanocomposites with r = 0, 1, 2 and 3 are approximated by using the same technique with fixed K, γ ∗ and  (these quantities are presumed to be independent of r). The effect of concentration of CNDs r on parameters G and κ is illustrated in Fig. 1C, where the data are approximated by the linear equations

G = G0 + G1 r,

κ = κ 0 + κ 1r

(20)

with the coefficients calculated by the least-squares method. Simulation of the loading path of the stress–strain diagram in Fig. 1D is performed with the constants collected in Table 1 and Fig. 1C (to reach good agreement with observations, G = 1.73 MPa is replaced with G = 1.67 MPa, the difference is less than 4%). Afterwards, we set α = 0, β1 = 0, and find P1∗ = 7.3 · 10−3 s−1 and R1 = 2.6 from the best-fit condition for the retraction path of the stress–strain curve in Fig. 1D. 3.1.2. IR-ATA-Zn elastomers We proceed with fitting observations on a series of supramolecular elastomers prepared by means of a three-stage procedure (Liu, Liu, et al., 2017). At the first stage, polyisoprene rubber (IR) is grafted with maleic anhydride (MAH) and 3-amino1,2,4-triazole (ATA) (the ATA concentration is 4.8 mol.%, and ATA:MAH proportion equals 2:1). At the other stage, ATA-grafted IR chains are dissolved in a solution of zinc chloride ZnCl2 in tetrahydrofuran (THF) to form metal-coordination bonds between zinc ions and triazole groups. At the final stage, THF is evaporated by drying the samples at 40 °C under vacuum. The permanent network in IR-ATA-Zn elastomers is formed due to entanglement between IR chains. Physical junctions in the transient network are formed by (i) hydrogen bonds between amide-triazole groups (ATA) and anhydrid moeties (MAH) and (ii) metal-ligand bonds between Zn ions and triazole groups (ATA). Mechanical tests were conducted at room temperature on elastomers with molar ratios r of ZnCl2 and ATA ranging from 0 to 0.5. The experimental program involves: (i) cyclic (loading–unloading) tests with strain rate ˙ = 0.0056 s−1 and maximum elongation ratio kmax = 5, and (ii) relaxation tests with strain rate under stretching ˙ = 0.083 s−1 and elongation ratio k = 2. Observations in these tests are reported in Fig. 2 together with results of simulation with the material constants collected in Table 2. In Fig. 2A, tensile stress σ is plotted versus elongation ratio k. In Fig. 2B, the dimensionless ratio of stresses S = σ (trel )/σ (trel = 0 ) is depicted versus relaxation time trel = t − t0 , where t0 is the instant when the required elongation ratio is reached. The loading paths of the stress–strain curves in Fig. 2A and the relaxation curves in Fig. 2B are characterized by 5 adjustable parameters: G, K, γ ∗ , κ and  . These quantities are found by using the same iterative algorithm as for the observations in Fig. 1. Three parameters (K, γ ∗ ,  ) that are independent of concentration of physical bonds are listed in Table 2. The effect of r on the remaining two parameters (G and κ ) is illustrated in Fig. 2C, where the data are approximated by Eq. (20).

318

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

Fig. 2. A: Tensile stress σ versus elongation ratio k. Symbols: experimental data in cyclic tests with strain rate ˙ = 0.0056 s−1 and maximum elongation ratio kmax = 5 on IR-ATA-Zn elastomers with various molar ratios r of ZnCl2 and ATA (Liu et al., 2017). Solid lines: results of simulation. B: Ratio of tensile stresses S versus relaxation time trel . Symbols: experimental data in relaxation tests with strain rate under stretching ˙ = 0.083 s−1 and elongation ratio k = 2. Solid lines: results of simulation. C: Parameters G and κ versus ratio r of molar fractions of ZnCl2 and ATA. Circles: treatment of observations. Solid lines: approximation of the data by Eq. (20). D: Parameters P1∗ and R1 versus ratio r. Circles: treatment of observations. Solid lines: approximation of the data by Eq. (21).

Setting α = 0 and β1 = 0, we conclude that each retraction path in Fig. 2A is determined by two parameters: P1∗ and R1 . Their values are found from the best-fit condition for each diagram in Fig. 1A separately. The effect of r on coefficients P1∗ and R1 is shown in Fig. 2D. The data are approximated by the linear relations

P1∗ = P10∗ + P11∗ r,

R1 = R01 + R11 r,

(21)

where the coefficients are calculated by the least-squares technique. 3.1.3. ENR-ZDA-Fe elastomers Finally, experimental data are approximated on supramolecular elastomers prepared on the basis of epoxidized natural rubber (ENR) (Zhang, Tang, Guo, & Zhang, 2016). To develop a transient network where physical junctions are formed by metal-oxigen coordination bonds, ENR is dissolved in a mixture of chloroform and iron chloride hexahydrate FeCl3 · 6H2 O. After evaporation of solvent, the master-batch was mixed with neat ENR and zinc diacrylate (ZDA) and hot-pressed at 160 °C. Covalent cross-links in the permanent network are formed under vulcanization of rubber with ZDA. Mechanical tests were performed at room temperature on ENR-ZDA-Fe elastomers with a fixed concentration of ZDA (2 phr with respect to rubber gum) and concentrations r of iron chloride ranging from 0 to 3 phr. The experimental program involves: (i) tensile tests with strain rate ˙ = 0.21 s−1 , and (ii) cyclic (loading–unloading) tests with strain rate ˙ = 0.042 s−1 and maximum elongation ratio kmax = 2. Observations in these tests are depicted in Fig. 3 (in Fig. 3A and B, tensile stress σ is plotted versus elongation ratio k under tension and cyclic deformation, respectively) together with results of numerical analysis with the material constants collected in Table 3. Adjustable parameters are determined by means of the following algorithm. First, the stress–strain diagram under tension of elastomer with r = 3 phr (Fig. 3A) is fitted with the help of 5 parameters (G, K, γ ∗ , κ ,  ). Then, the triplet (K, γ ∗ ,  ) is fixed, and the other diagrams in Fig. 1A and the loading paths of the stress–strain curves in Fig. 3B are approximated. Each diagram is matched separately with the help of two parameters (G and κ ). The effect of r on these quantities is reported in Fig. 3C, where the data are approximated by Eq. (20).

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

319

Fig. 3. A: Tensile stress σ versus elongation ratio k. Symbols: experimental data in tensile tests with strain rate ˙ = 0.21 s−1 on ENR-ZDA-Fe elastomers with various concentrations r phr of iron chloride (Zhang et al., 2016). Solid lines: results of simulation. B: Tensile stress σ versus elongation ratio k. Symbols: experimental data in cyclic tests with strain rate ˙ = 0.042 s−1 and maximum elongation ratio kmax = 2. Solid lines: results of simulation. C: Parameters G and κ versus concentration r. Symbols: treatment of observations in tensile and cyclic tests. Solid lines: approximation of the data by Eq. (20). D: Parameter P1∗ versus concentration r. Circles: treatment of observations in cyclic tests. Solid line: approximation of the data by Eq. (21). Table 3 Material parameters for ENR-ZDA-Fe elastomers (Fig. 3). Parameter

Value

K

34.0 3.0 0.8 0.0 −0.2

γ ∗ (s−1 )  α β1

Afterwards, the retraction paths of the curves in Fig. 3B are fitted. Setting α = 0, we calculate β 1 and R1 from the best-fit condition for the unloading curve of the elastomer with r = 1. These values are fixed, and each retraction path in Fig. 3B is matched separately with the help of the only parameter P1∗ . Evolution of this quantity with r is illustrated in Fig. 3D, where the data are approximated by Eq. (21). 3.1.4. Discussion Figs. 1–3 demonstrate good agreement between results of simulation and experimental data on IR-CND, IR-ATA-Zn and ENR-ZDA-Fe supramolecular elastomers in tensile tests, relaxation tests, and cyclic tests. Small deviations are observed in Fig. 3A only for the elastomer with r = 3 at relatively small strains. Adjustable parameters in the governing equations evolve consistently with concentration of physical junctions. Some scatter, however, should be mentioned in Figs. 2D and 3C between the data and theoretical predictions. Inaccuracy of measurements appears to be one of the reasons for these discrepancies. For example, tensile stress in the elastomer with r = 1 along the loading path of the cyclic test (Fig. 3B) exceeds noticeably that under tension (Fig. 3A) despite the fact that the strain rate under tensile deformation exceeds that under cyclic loading by a factor of five.

320

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335 Table 4 Material constants for sPS-IR-sPS triblock copolymer (Fig. 4). Parameter

Value

K

75.0 0.2 1.6 0.0 0.3 0.24 4.0 23.0 105.0

γ ∗ (s−1 )  α β1 σ¯ (MPa) R2 · 10−3

β2 β

Results of numerical analysis show that the mechanical response of supramolecular elastomers is described adequately by the model under the assumption that (i) the distribution of temporary bonds with various activation energies (characterized by  ), (ii) the attempt rate for separation of chains from temporary bonds γ ∗ , and (iii) the Gent constant K are affected insignificantly by concentration of extra physical junctions r. An increase in concentration of these junctions results in the growth of the elastic modulus G (Figs. 1C, 2C and 3C). It follows from Eq. (A-67) that G equals the sum of the elastic moduli of permanent, Ge , and transient, Ga , networks. An increase in G with r is explained by two reasons: (i) Ga grows with r being proportional to this parameter, and (ii) Ge increases with r as some non-covalent bonds whose lifetime exceeds the duration of test should be treated as permanent junctions. The growth of concentration of non-covalent junctions r can induce both an increase and a decrease in coefficient κ . This seems to be natural as κ is defined by Eq. (A-67) as the ratio of Ga to the sum Ga + Ge . When Ga increases noticeably, while Ge remains practically unaffected by r, the parameter κ grows (Figs. 1C and 3C). On the contrary, when Ge increases strongly with r (for example, due to replacement of “weak” hydrogen bonds having short lifetimes with “strong” metal-coordination bonds whose lifetime exceeds the duration of test), κ is reduced with r (Fig. 2C). Figs. 2D and 3D demonstrate that the rate of plastic flow under retraction P1∗ increases with r. This conclusion is in accord with the mechanism of plastic deformation (Filippidi et al., 2017), according to which sliding of permanent junctions with respect to their reference positions occurs when these junctions become unbalanced due to sudden relaxation of stresses in chains separated from physical bonds (which means that the higher is the concentration of non-covalent junctions, the more pronounced plastic flow is observed). Fig. 2D reveals that the growth of concentration of metal-ligand bonds r induces an increase in energy of inter-chain interaction R1 . This result appears to be natural: the energy of interaction between chains grows when “weak” hydrogen bonds are replaced with “strong” metal-coordination bonds. On the contrary, good agreement with observations in Fig. 3B is reached under the assumption that R1 is unaffected by r. This difference can be explained by the structure of permanent networks in these materials: permanent junctions in IR-ATA-Zn elastomers are formed by entanglements whose lifetime exceeds the duration of tests, while those in ENR-ZDA-Fe elastomers are formed by covalent cross-links developed under vulcanization. 3.2. Mechanical response in multi-cycle tests For two series of triblock copolymers (with varied lengths of blocks), we show that the governing equations ensure good agreement with observations in multi-cycle tests with monotonically increasing maximum elongation ratios kmax n and the zero stress under retraction, and examine the ability of the model to predict the mechanical response in cyclic tests with different deformation programs (oscillations between a fixed maximum elongation ratio under stretching kmax and the zero minimum stress under retraction, and ratcheting between the maximum σ max and minimum σ min tensile stresses). 3.2.1. sPS-IR-sPS triblock copolymer We begin with fitting observations on sPS-IR-sPs triblock copolymer prepared by cationic copolymerization of syndiotactic polystyrene (sPS) and hydrogenated polyisoprene (IR) with composition 1:8:1, number average molecular weight 105 kg/mol, polydispersity index 1.18, and volume fraction of sPs 0.27 (Higaki et al., 2017). To develop a double netwotk in the triblock copolymer, it was dissolved in dichlorobenzene at 140 °C, dried under vacuum at the same temperature (at this stage, the sPS chains were aggregated to form semicrystalline domains serving as permanent and transient junctions between IR chains), and annealed at 120 °C for 12 h. Multi-cycle tensile tests were performed at room temperature with a fixed strain rate ˙ = 0.0167 s−1 under tension and retraction. A sample is stretched up to the maximum elongation ratio kmax 1 , unloaded down to the zero stress σmin = 0, stretched up to the maximum elongation ratio kmax 2 , unloaded down to the zero stress σmin = 0, etc. Observations in the test with kmax 1 = 2, kmax 2 = 2, kmax 3 = 3, kmax 4 = 3, kmax 5 = 4, kmax 6 = 4 are depicted in Fig. 4A together with results of simulation with the material constants listed in Table 4.

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

321

Fig. 4. A: Tensile stress σ versus elongation ratio k. Symbols: experimental data in multi-cycle test on sPS-IR-sPS triblock copolymer with strain rate ˙ = 0.0167 s−1 and various kmax n (Higaki et al., 2017). Solid line: results of simulation. B: Parameters P1∗ and R1 versus the principal invariant Ie2 at maximum elongation ratio kmax . Circles: treatment of observations. Solid lines: approximation of the data by Eq. (17). C: Parameter P2∗ versus the principal invariant Ip2 at minimum elongation ratio kmin . Circles: treatment of observations. Solid line: approximation of the data by Eq. (18). D: Tensile stress σ versus elongation ratio k. Solid lines: predictions of the model. Upper panel: program (22); lower panel: program (23).

Adjustable parameters are determined by means of the following algorithm. The parameters G, K, γ ∗ , κ ,  are found from the best-fit condition for the envelope of all loading paths in Fig. 4A. Then, we set α = 0 and determine P1∗ , R1 and β 1 from the best-fit condition for the retraction path of the curve with kmax 6 . Fixing β 1 , we match unloading paths of the stress–strain curves with other kmax n with the help of two parameters, P1∗ and R1 . Each path is fitted separately. The stress σ¯ is determined from the condition that the curvature of the reloading paths vanishes. After approximation of the stress– strain diagram for the cycle with kmax 5 , the reloading path along the last cycle is fitted by means of parameters P2∗ , R2 , β 2 , β . Then, R2 , β 2 and β are fixed, and all other reloading paths are matched separately with the help of the only parameter P2∗ . For each cycle, the value ke max is calculated corresponding to k = kmax n (by integration of the governing equations), and the quantities P1∗ and R1 are plotted versus Ie2 max in Fig. 4B together with their approximation by Eq. (17). In a similar way, we calculate for each cycle kp min at the instant when tensile stress σ vanishes, determine Ip2 min at this point, and plot P2∗ versus Ip2 min in Fig. 4C. The data in this figure are matched by Eq. (18). Coefficients in Eqs. (17) and (18) are determined by the least-squares method. Eq. (18) is chosen to fit the experimental data (instead of Eq. (19)) because the quantity Ie2 min − 3 remains rather small for all cycles of loading–retraction. The latter can be explained by the fact that the relaxation rate γ∗ = 0.2 s−1 for the copolymer under consideration exceeds strongly the strain rate ˙ = 0.0167 s−1 . To show the ability of the model to predict the Mullins effect, numerical analysis is conducted of the governing equations for deformation programs with 3 and 6 cycles:

kmax n = 1 + n

kmax n = 1 +

1 n 2

( n = 1, . . . , 3 ),

( n = 1, . . . , 6 ).

(22)

(23)

322

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

Fig. 5. Results of simulation for sPS-IR-sPS triblock copolymer in multi-cycle tests with strain rate ˙ = 0.0167 s−1 . A: Tensile stress σ versus elongation ratio k for program P-1 and various kmax . B: Maximum stress per cycle σ max versus number of cycles n for program P-1 and various kmax . C: Tensile stress σ versus elongation ratio k for program P-2 and various σ max MPa. D: Maximum elongation ratio per cycle kmax versus number of cycles n for program P-2 with various σ max MPa.

Results of simulation are reported in Fig. 4D. This figure reveals all characteristic features of the Mullins phenomenon: (i) at each cycle of deformation, the loading and unloading paths differ noticeably (stress-softening), (ii) when elongation ratio k exceeds the maximum elongation ratio at the previous cycle, the reloading curve coincides with the stress–strain curve under tension (the Mullins effect), (iii) two retraction curves corresponding to different deformation programs coincide when they start from the same kmax (fading memory of deformation history). To demonstrate that the mechanical behavior in multi-cycle tests with other deformation programs can be predicted adequately by the governing equations, two programs are chosen: • P-1: oscillations between a fixed maximum elongation ratio kmax and the zero minimum stress σmin = 0, • P-2: ratcheting between a fixed maximum tensile stress σ max and the zero minimum stress. Results of numerical analysis for multi-cycle (20 cycles) tests are presented in Fig. 5, where stress–strain diagrams with various kmax and σ max are depicted in Fig. 5A and C, respectively. Fig. 5B reveals a decay in the maximum stress per cycle σ max with number of cycles n for program P-1. Fig. 5D demonstrates an increase in the maximum elongation ratio kmax with number of cycles n for program P-2. The data in these figures are in qualitative agreement with observations reported in Bartolome, Aurrekoetxea, Urchegui, and Tato (2013) and Wang et al. (2018, 2017), to mention a few. 3.2.2. IBA-(BA-co-VI)-IBA triblock copolymers We proceed with fitting experimental data on a series of ABA triblock copolymers with A-block formed by poly(isobornyl acrylate) (IBA) and B-block formed by copolymer of poly(n-butyl acrylate) (BA) and poly(1-vinylimidazole) (VI). The copolymers were synthesized by means of RAFT (reversible addition-fragmentation chain transfer) technique (Jiang et al., 2017). Elastomers designated here as E−1, E−2 and E−3 have the same content of IBA (32 wt.%) and VI (2 wt.%), but differ by chain lengths. Their number-average molecular weights read 136, 66 and 34 kg/mol with polydispercity indices 1.6, 1.4 and 1.3, respectively. Due to phase separation under preparation conditions, aggregation of IBA blocks induces development of glassy micro-domains that serve as physical bonds between BA-co-VI chains in the rubbery state. Entanglements between

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

323

Fig. 6. A, B, C: Tensile stress σ versus elongation ratio k. Symbols: experimental data in multi-cycle tests with various kmax n on elastomers E−1, E−2 and E−3 (Liang et al., 2017). Solid lines: results of simulation. D: Parameter P1∗ versus the principal invariant Ie2 at maximum elongation ratio kmax . Symbols: treatment of observations. Solid lines: approximation of the data by Eq. (17).

Table 5 Material parameters for triblock copolymers E−1, E−2 and E−3 (Fig. 6). Parameter

E−1

E−2

E−3

G (MPa) K γ ∗ (s−1 )

0.43 46.0 0.06 0.53 0.8 2.3 0.1 0.5 0.31 8.0 8.0 70.0

0.205 70.0 0.06 0.55 0.8 2.3 0.4 0.5 0.132 24.0 23.0 195.0

0.092 51.0 0.06 0.53 0.8 2.3 0.3 0.5 0.044 200.0 67.0 550.0

κ  α

R1

β1 σ¯ (MPa)

R2 · 10−2

β2 β

chains and glassy micro-domains whose lifetime exceeds the characteristic time of experiments are treated as permanent junctions, while those whose lifetime is lower than the duration of tests are modeled as temporary bonds. Multi-cycle tensile tests were conducted at room temperature with a fixed strain rate ˙ = 0.0083 s−1 under tension and retraction. In each experiment, a sample is stretched up to elongation ratio kmax 1 , retracted down to the zero stress σmin = 0, reloaded up to elongation ratio kmax 2 > kmax 1 , unloaded down to σmin = 0, etc. Observations of elastomers E−1, E−2 and E−3 in tests with kmax 1 = 1.5, kmax 2 = 2.0, kmax 3 = 2.5, kmax 4 = 3.0, kmax 5 = 3.5, kmax 6 = 4.0 are plotted in Fig. 6A–C, respectively, together with results of simulation with the material constants listed in Table 5. Adjustable parameters are determined by the same algorithm that was used to match the data in Fig. 4A. The only difference is that the fitting procedure is conducted under the assumption that R1 remains independent of kmax n . This implies that each retraction path is determined by the only parameter P1∗ , and each reloading path is characterized by the only

324

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

Fig. 7. A: Parameter P2∗ versus the principal invariant Ie2 at minimum elongation ratio kmin . Symbols: treatment of observations. Solid lines: approximation of the data by Eq. (19). B: Parameter P2∗ versus the principal invariant Ip2 at the point with elongation ratio kmin . Symbols: treatment of observations. Solid lines: approximation of the data by Eq. (18). C: Tensile stress σ versus elongation ratio k. Solid lines: results of simulation with Eq. (19) for elastomer E−1. Upper panel: program (22); lower panel: program (23). D: Tensile stress σ versus elongation ratio k. Solid lines: results of simulation with Eq. (18) for elastomer E−1. Upper panel: program (22); lower panel: program (23). E: Tensile stress σ versus elongation ratio k. Solid lines: results of simulation with Eq. (19) for elastomer E−1 in tests with various maximum elongation ratios kmax . F: Tensile stress σ versus elongation ratio k. Solid lines: results of simulation with Eq. (18) for elastomer E−1 in tests with various maximum elongation ratios kmax .

coefficient P2∗ . Evolution of P1∗ with maximum elongation ratio per cycle is shown in Fig. 6D. The data are approximated by Eq. (17) with the coefficients calculated by the least-squares technique. The data collected in Table 5 demonstrate that an increase in molecular weight of triblock chains (i) induces a pronounced growth of the elastic modulus G and the stress σ¯ at which the curvature of reloading curves vanishes, (ii) results in a strong decrease in parameter R2 (characterizing the energy of inter-chain interaction) and coefficients β 2 and β (that account for the effect of stress on the rate of plastic flow under reloading), and (iii) does not affect substantially the other parameters in the governing equations.

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

325

Fig. 8. A, C, E: Maximum stress per cycle σ max versus number of cycles n. Symbols: results of simulation for multi-cycle tests with various maximum elongation ratios kmax on elastomers E−1, E−2, E−3. B,D,F: Maximum elongation ratio per cycle kmax versus number of cycles n. Symbols: results of simulation for multi-cycle tests with various maximum stresses σ max MPa on elastomers E−1, E−2, E−3.

To assess the difference between Eqs. (18) and (19), the parameter P2∗ is plotted versus ke min in Fig. 7A and versus kp min in Fig. 7B. The data are approximated by Eqs. (18) and (19), where the coefficients are found by the least-squares method. The same accuracy of fitting observations is revealed in these figures. The response of triblock copolymers in cyclic tests with monotonically increasing maximum elongation ratios kmax n is analyzed numerically. Results of simulation for elastomer E−1 in tests with programs (22) and (23) are reported in Fig. 7C and D. These figures show that assumptions (18) and (19) lead to similar stress–strain diagrams demonstrating the characteristic features of the Mullins effect. We proceed with predictions of the response in multi-cycle (20 cycles) tests with deformation program P-1 and kmax = 1.5, 2.0, 2.5 and 3.0. Results of numerical analysis for elastomer E−1 are reported in Fig. 7E and F. Similar dependencies for elastomers E−2 and E−3 are not presented for the sake of brevity. Fig. 7E and F demonstrate a pronounced difference between the stress–strain curves based on Eqs. (18) and (19). Simulation grounded on Eq. (19) ensures that the

326

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

maximum stress per cycle σ max increases monotonically with kmax and decreases weakly with n. Numerical analysis based on Eq. (18) shows that the dependence of σ max on kmax becomes non-monotonic, and maximum stress per cycle is strongly reduced with n at kmax > 2. As such a response is not observed in experiments, we conclude that Eq. (18) cannot be applied to predict the mechanical behavior of triblock copolymers in multi-cycle tests with arbitrary deformation programs. Finally, we apply Eq. (19) to predict the response of elastomers E−1, E−2 and E−3 in multi-cycle (20 cycles) tests with program P-1 (with kmax = 1.5, 2.0, 2.5, 3.0) and P-2 (with σ max ranging from 0.2 to 0.8 MPa for E−1, from 0.1 to 0.25 MPa for E−2, and from 0.04 to 0.1 MPa for E−3; this choice ensures that kmax does not exceed 4.0 in ratcheting tests). Results of numerical analysis are reported in Fig. 8. In Fig. 8A, C, E, the maximum stress per cycle σ max is plotted versus number of cycles n for program P-1. In Fig. 8B, D, F, the maximum elongation ratio per cycle kmax is depicted versus n for program P-2. The following conclusions are drawn from these figures: (I) Under deformation program P-1, σ max decreases noticeably with n along the first five cycles and remains constant afterwards. The rate of decay in σ max during the transition period increases slightly with kmax . The ultimate maximum stress per cycle after equilibration grows strongly with kmax . These conclusions are in agreement with experimental data on thermoplastic elastomers (Beniah et al., 2017; Drozdov & Christiansen, 2008) and double-network gels (Bai, Yang, Morelle, Yang, & Suo, 2018; Bai et al., 2017; Yang, Wang, Yang, Wang, & Wu, 2018; Zhang et al., 2018). (II) Under deformation program P-2, kmax increases noticeably with n along the first five cycles (primary fatigue). After this transition period, kmax remains practically independent of n at small σ max and grows linearly with number of cycles at moderate σ max (secondary fatigue). The rate of growth of kmax with n at the stage of secondary fatigue increases substantially with maximum stress per cycle σ max . These results are in accord with observations on thermoplastic elastomers (El Fray & Altstadt, 2003; Gotz, Lim, Puskas, & Altstadt, 2014).

4. Concluding remarks A constitutive model is developed for the viscoelastic and viscoplastic behavior of supramolecular elastomers with covalent cross-links and non-covalent bonds formed due to electrostatic interaction, metal-ligand coordination and hydrophobic association. An elastomer is treated as a combination of two networks: chains in the first network are linked by permanent junctions, while chains in the other network are bridged by temporary bonds. The viscoelastic response reflects breakage and reformation of temporary junctions in the transient network (transition of chains connected by physical bonds from their active to dangling state). The transient network is presumed to be inhomogeneous and composed of meso-regions with different rates for rearrangement of junctions. This allows the entire relaxation spectrum to be accounted for in terms of a distribution function for meso-regions with various activation energies for separation of chains. The viscoplastic response reflects slippage of junctions in the permanent network with respect to their reference positions. A junction starts to slide when it becomes unbalanced due to sudden relaxation of stresses in one of the chains connected by this junction caused by transformation of this chain from its active state into the dangling state. The sliding process (plastic flow) proceeds until the junction reaches its new equilibrium state. Kinetic equations for plastic flow under isothermal conditions are derived from the free energy imbalance inequality with account for the energy of inter-chain interaction. Under multi-cycle deformation, coefficients in these relations adopt different values for the first loading, retraction and reloading. The model is applied to fit experimental data in tensile relaxation tests, loading-unloading tests, and multi-cycle tests on a series of supramolecular elastomers and triblock copolymers with covalent and non-covalent bonds. Two characteristic features of these observations are to be mentioned: (i) the rate of relaxation is comparable with the strain rate under stretching and retraction, and (ii) the reloading paths of the stress–strain diagrams are S-shaped: these curves are convex at small stresses and concave at large stresses. Numerical simulation demonstrates that (i) the governing equations describe adequately the experimental stress–strain diagrams, (ii) the material parameters evolve consistently with chemical composition of elastomers (Figs. 1C, 2C and 3C, D) and experimental conditions (Figs. 4B, C, 6D and 7A, B), and (iii) predictions of the model are in quantitative (Figs. 4D and 7C, D) and qualitative (Fig. 8) agreement with observations. Fitting observations on supramolecular elastomers reveals that an increase in concentration of temporary bonds r induces pronounced changes in two parameters: G and κ . The elastic modulus G grows linearly with r. The coefficient κ that characterizes the ratio of the elastic modulus of the transient network to that of the elastomer increases (Figs. 1C and 3C) or decreases (Fig. 2C) depending on strength of non-covalent bonds. Matching experimental data of triblock copolymers with fixed weight fractions of blocks shows that the elastic modulus G increases strongly with molecular weight of chains. Although the constitutive equations can be applied to the analysis of the mechanical response of supramolecular elastomers under an arbitrary three-dimensional loading, this study focuses on tensile cyclic tests with a fixed strain rate, but various (stress- and strain-controlled) deformation programs. An important advantage of the model is its ability to predict observations in multi-cycle tests with one program when material constants are determined by fitting data in a test with another program (Fig. 8). Numerical simulation demonstrates that this property is ensured when the rates of plastic flow P1∗ and P2∗ are governed by the elastic energy stored in the permanent network (Eqs. (17) and (19)).

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

327

Acknowledgments Financial support by the Innovationsfonden (IFD project 5152-0 0 0 02B) is gratefully acknowledged. ADD is thankful to Profs. O. Okay and C.S. Patrickios for fruitful discussions. Appendix A A supramolecular elastomer is modeled as an equivalent network of polymer chains. The network is treated as a superposition of two networks: permanent and transient. The permanent network is formed by chains bridged by covalent cross-links. Chains in the transient network are connected by physical junctions whose rearrangement is driven by thermal fluctuations. A1. Kinematic relations The initial configuration of an elastomer coincides with that of an undeformed specimen. We adopt the affinity hypothesis and presume deformations of permanent and transient networks to coincide with macro-deformation of the elastomer. Transformation of the initial configuration into the actual configuration at time t is described by the deformation gradient F(t, X), where X stand for radius-vector of an arbitrary point. To simplify the notation, some arguments of functions will be omitted. A1.1. Permanent network The reference (stress-free) state of the permanent network before loading coincides with the initial state of the elastomer. Denote by Fe and Fp the deformation gradients for elastic and plastic deformations of the permanent network. These tensors are connected with the deformation gradient for macro-deformation F by the multiplicative decomposition formula (1). The velocity gradient for macro-deformation is given by

L = F˙ · F−1 .

(A-1)

Combination of Eqs. (1) and (A-1) implies that

L = Le + Lp ,

(A-2)

where

Le = F˙ e · F−1 e ,

Lp = Fe · lp · F−1 e ,

lp = F˙ p · F−1 p .

(A-3)

Following the conventional approach, we disregard plastic spin and presume the velocity gradient lp to be symmetric,

lp = dp ,

(A-4)

where

1 (lp + lp ) 2

dp =

(A-5)

is the rate-of-strain tensor for plastic deformation. It follows from Eqs. (A-3) and (A-4) that the tensor Fp obeys the differential equation

F˙ p = dp · Fp

(A-6)

with the initial condition

Fp ( 0 ) = I,

(A-7)

where I is the unit tensor. Combination of Eqs. (A-2)–(A-4) results in

D = De + Dp ,

(A-8)

where

D=

1 ( L + L ), 2

with

Dp =

De =

1 (Le + Le ), 2



Dp =

1 (Lp + Lp ) 2

(A-9)



1 −  Fe · dp · F−1 e + Fe · dp · Fe . 2

(A-10)

The left Cauchy–Green tensors for elastic and plastic deformation read

Be = Fe · F e,

bp = Fp · F p.

(A-11)

328

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

Differentiating Eqs. (A-11) with respect to time and using Eqs. (A-3) and (A-6), we find that

B˙ e = Le · Be + Be · L e,

b˙ p = dp · bp + bp · dp .

(A-12)

Denote by Ie1 , Ie2 , Ie3 the principal invariants of the Cauchy–Green tensor for elastic deformation, and by Ip1 , Ip2 , Ip3 the principal invariants of the Cauchy–Green tensor for plastic deformation. It follows from Eq. (A-12) that the derivatives of these functions with respect to time obey the differential equations





I˙e1 = 2Be : De ,

I˙e2 = 2 Ie2 I − Ie3 B−1 : De , e

I˙p1 = 2bp : dp ,

I˙p2 = 2 Ip2 I − Ip3 b−1 : dp , p





I˙e3 = 2Ie3 I : De ,

(A-13)

I˙p3 = 2Ip3 I : dp ,

(A-14)

where the colon stands for convolution. Replacing De in Eq. (A-13) by means of Eq. (A-8) and using Eq. (A-10), we find that





I˙e1 = 2 Be : D − Ce : dp , I˙e2 = 2











Ie2 I − Ie3 B−1 : D − Ie2 I − Ie3 C−1 : dp , e e





I˙e3 = 2Ie3 I : D − I : dp ,

(A-15)

Ce = F e · Fe

(A-16)

where

stands for the right Cauchy–Green tensor for elastic deformation. Eqs. (A-14) and (A-15) are valid for arbitrary elastic and plastic deformations. We suppose now that these deformations are volume-preserving (isochoric), which means that

Ie3 = 1,

Ip3 = 1.

(A-17)

Combination of Eq. (A-17) with Eqs. (A-14) and (A-15) results in

I : D = 0,

I : dp = 0,

(A-18)

which implies that the tensors D and dp are traceless. It follows from Eqs. (A-14) and (A-18) that the functions Ip1 , Ip2 are governed by the equations

I˙p2 = −2b−1 p : dp ,

I˙p1 = 2bp : dp ,

and the functions Ie1 , Ie2 obey the equations





I˙e1 = 2 Be : D − Ce : dp ,

(A-19)





−1 I˙e2 = −2 B−1 e : D − Ce : dp .

(A-20)

Differentiating Eq. (A-16) with respect to time and using Eqs. (A-3) and (A-9), we find that

C˙ e = 2F e · De · Fe .

(A-21)

Combination of Eqs. (A-21) with Eqs. (A-8) and (A-10) yields

C˙ e = 2F e · D · Fe − ( Ce · dp + dp · Ce ).

(A-22)

A1.2. Transient network To describe the time-dependent behavior of the transient network, we denote by τ an instant when a dangling chain bridges with the network, and distinguish chains bonded to the network before loading (τ = 0) and those attaching the network under deformation (τ > 0). For a chain with τ = 0, the reference state coincides with the initial state, and deformation gradient for elastic deformation f0 (t) is given by

f0 (t ) = F(t ).

(A-23)

For a chain with τ > 0, the reference state is presumed to coincide with the actual state of the network at instant τ . This means that stresses in a dangling chain have relaxed entirely before it joins the network, and the deformation gradient for elastic deformation reads

fτ (t ) = F(t ) · F−1 (τ ).

(A-24)

The corresponding Cauchy–Green tensors are determined by the formulas

b0 = f0 · f 0,

bτ = fτ · f τ.

(A-25)

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

329

Differentiating Eqs. (A-25) with respect to time and using Eqs. (A-3), (A-23) and (A-24), we find that

b˙ 0 = L · b0 + b0 · L ,

b˙ τ = L · bτ + bτ · L .

(A-26)

Denote by I01 , I02 , I03 and Iτ 1 , Iτ 2 , Iτ 3 the principal invariants of the tensors b0 and bτ , respectively. These functions are governed by the differential equations









I˙01 = 2b0 : D,

I˙02 = 2 I02 I − I03 b−1 : D, 0

I˙03 = 2I03 I : D,

I˙τ 1 = 2bτ : D,

I˙τ 2 = 2 Iτ 2 I − Iτ 3 b−1 : D, τ

I˙τ 3 = 2Iτ 3 I : D.

Keeping in mind Eq. (A-18), we conclude from these relations that

I03 = 1,

Iτ 3 = 1,

(A-27)

and the functions I01 , I02 , Iτ 1 , Iτ 2 obey the equations

I˙01 = 2b0 : D,

I˙02 = −2b−1 0 : D,

I˙τ 1 = 2bτ : D,

I˙τ 2 = −2b−1 τ : D.

(A-28)

To describe rearrangement of chains in the transient network (separation of active chains from their junctions and merging of dangling chains with the network), we presume the network to be composed of meso-domains with various activation energies U. The rate of separation of an active chain from a junction  in a meso-domain with activation energy U is governed by the Eyring equation

 U   = γ exp − , kB T

where γ denotes the attempt rate. Introducing the dimensionless activation energy u = U/(kB T ), we present this relation in the form

(t, u ) = γ (t ) exp(−u ),

(A-29)

where the dependence of γ on time reflects the effect of mechanical factors on the rearrangement process. The number of active chains in the transient network (per unit volume in the initial state) is denoted by Na . The number of active chains in meso-domains with dimensionless activation energy u reads

N¯ a (u ) = Na f (u ),

(A-30)

where the distribution function for meso-domains with various activation energies f(u) obeys normalization condition (3). The current state of the transient network is uniquely determined by a function of three variables n(t, τ , u) that equals the number (per unit volume) of active chains at time t ≥ 0 that have returned into the active state before instant τ ≤ t and belong to a meso-domain with activation energy u. In particular, n(t, 0, u) denotes concentration of chains in meso-domains with activation energy u that have been attached to the network before loading and remain active at time t, while n(t, t, u) stands for concentration of active chains at time t ≥ 0 in meso-domains with activation energy u. We suppose that the number of active chains in meso-domains with various activation energies remains independent of time,

n(t , t , u ) = n(0, 0, u ). It follows from this equality and Eq. (A-30) that

n(t , t , u ) = Na f (u ).

(A-31)

The number of chains that were active at the initial instant and detach from their junctions within the interval [t , t + dt ] reads

−

∂n (t, 0, u ) dt, ∂t

the number of dangling chains that return into the active state within the interval [τ , τ + dτ ] is given by ϕ (τ , u)dτ with

ϕ (τ , u ) =

 ∂n  (t, τ , u ) , ∂τ t=τ

(A-32)

and the number of chains that merged (for the last time) with the network within the interval [τ , τ + dτ ] and detach from their junctions within the interval [t , t + dt ] is determined as

−

∂ 2n (t, τ , u ) dtdτ . ∂t∂ τ

Separation of active chains from their junctions is governed by the first-order kinetic equations

∂n (t, 0, u ) = −(t , u )n(t , 0, u ), ∂t

∂ 2n ∂n (t, τ , u ) = −(t , u ) (t , τ , u ), ∂t∂ τ ∂τ

(A-33)

330

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

which mean that the number of active chains detaching from temporary junctions per unit time is proportional to the number of active chains in an appropriate meso-region. Integrating Eq. (A-33) with initial conditions (A-31) and (A-32), we find that

 

n(t, 0, u ) = Na f (u ) exp −

t 0

 (s, u )ds ,

 t  ∂n (t, τ , u ) = ϕ (τ , u ) exp − (s, u )ds . ∂τ τ

(A-34)

Inserting expressions (A-34) into the equality

n(t , t , u ) = n(t, 0, u ) +



∂n (t, τ , u )dτ , ∂τ

t

0

and using Eq. (A-31), we find that

ϕ (t, u ) = γ (t )Na f (u ).

(A-35)

A2. Strain energy density Denote by W the strain energy density (per unit volume) of the equivalent polymer network. For the network consisting of two parts (permanent and transient), the function W is given by

W = W1 + W2 ,

(A-36)

where W1 , W2 are the specific mechanical energies of chains that form the permanent and transient networks. The strain energy density of the permanent network reads

W1 = We + Wp ,

(A-37)

where the mechanical energy stored in chains We depends on the principal invariants Ie1 , Ie2 of the Cauchy–Green tensor for elastic deformation Be , and the energy of inter-chain interaction Wp is treated as a function of the principal invariants Ip1 , Ip2 of the Cauchy–Green tensor bp . To calculate the strain energy density of the transient network W2 , we introduce the mechanical energy per chain w and suppose that it depends on the principal invariants Iτ 1 , Iτ 2 of the Cauchy–Green tensor bτ (τ ≥ 0). At an arbitrary instant t ≥ 0, the quantity W2 is determined by



W2 (t ) =

∞ 0







n(t, 0, u )w I01 (t ), I02 (t ) +



   ∂n (t , τ , u )w Iτ 1 (t ), Iτ 2 (t ) dτ du, ∂τ

t

0

(A-38)

where the first term stands for the strain energy (per unit volume) of chains that have merged with the network before deformation and remain active at time t, and the other term expresses the energy stored in chains that have attached the network at various instants τ ≤ t and have not separated from the network within the intervals [τ , t]. A3. Constitutive equations Constitutive equations for a supramolecular elastomer under isothermal deformation are developed by means of the free energy imbalance inequality

W˙ − mec ≤ 0,

(A-39)

mec = J : D

(A-40)

where

denotes the work (per unit volume and unit time) produced by external loads, and  stands for the Cauchy stress tensor. Eq. (A-39) is fulfilled when the deformation gradient F obeys the incompressibility condition (A-18). To account for the latter condition, we multiply Eq. (A-18) by an arbitrary function  (pressure treated as a Lagrange multiplier) and add the result to Eq. (A-39). Using Eqs. (A-19), (A-20), (A-28), (A-33), (A-36)–(A-38) and (A-40), and arrive at inequality

[2Ke − J ( + I )] : D + 2Kp : dp −  ≤ 0 with



(A-41)



Ke = We,1 Be − We,2 B−1 e  ∞    + n(t, 0, u ) w,1 b0 − w,2 b−1 + 0



0





0

   ∂n (t, τ , u ) w,1 bτ − w,2 b−1 dτ du, τ ∂τ 

− We,1 Ce − We,2 C−1 e  ∞    t = (t, u ) n(t, 0, u )w I01 (t ), I02 (t ) +

Kp =

Wp,1 bp − Wp,2 b−1 p

t

0

0

,

   ∂n (t, τ , u )w Iτ 1 (t ), Iτ 2 (t ) dτ du, ∂τ

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

331

where

∂ We , ∂ Iem

We,m =

Wp,m =

∂ Wp , ∂ Ipm

w,m =

∂w , ∂ Iτ m

and the prime stands for the deviatoric component of a tensor (A = A − 13 A : I for any A). Keeping in mind that D is an arbitrary tensor function, while dp is an arbitrary tensor function that obeys Eq. (A-18), and using the non-negativity of the function , we conclude that inequality (A-41) is satisfied, provided that the Cauchy stress tensor  is given by

 = − I + 2  +





We,1 Be − We,2 B−1 + e



∞



0



n(t, 0, u ) w,1 b0 − w,2 b−1 0



  

∂n (t, τ , u ) w,1 bτ − w,2 b−1 dτ du , τ ∂τ

t 0

(A-42)

and plastic deformation of the permanent network is governed by the equation

dp = P







We,1 Ce − We,2 C−1 − Wp,1 bp − Wp,2 b−1 e p



,

where P is an arbitrary non-negative function. This relation is equivalent to the equation

dp = P −







We,1 Ce − We,2 C−1 − Wp,1 bp − Wp,2 b−1 e p

1 3









We,1 Ie1 − We,2 Ie2 − Wp,1 Ip1 − Wp,2 Ip2

 

I .

(A-43)

Constitutive Eqs. (A-42) and (A-43) are accompanied by the equilibrium equation for the stress tensor

∇ ·  = 0,

(A-44)

where ∇ stands for the gradient operator in the actual state, and the corresponding initial and boundary conditions. For definiteness, we adopt the Gent equation for the strain energy density of the permanent network and the neo– Hookean formulas for the specific energy of inter-chain interaction and the strain energy of active chains in the transient network, see Eq. (5). Insertion of Eq. (5) into Eq. (A-43) implies that

dp = where

 





1 1 1 PGe V Ce − Ie1 I − R bp − Ip1 I 2 3 3



V = 1−

Ie1 − 3 K

−1

,

R=



,

(A-45)

Gp . Ge

(A-46)

Combination of Eqs. (A-12) and (A-45) results in the differential equation for the tensor bp ,









1 1 b˙ p = PGe V Ce · bp + bp · Ce − (Ie1 − RIp1 )bp − Rb2p . 2 3

(A-47)

Inserting expression (A-45) into Eq. (A-22), we arrive at the differential equation for the tensor Ce ,

C˙ e = 2F e · D · Fe + PGe

1  2



R Ce · bp + bp · Ce +



1 (V Ie1 − RIp1 )Ce − V C2e . 3

(A-48)

Combination of Eqs. (A-34), (A-35), (A-42) and (5) results in the constitutive equation

 ∞   t   = − I + Ge (V Be − I ) + Ga f (u ) exp − (s, u )ds b0 0 0  t  t   + (τ , u ) exp − (s, u )ds bτ dτ − I du

(A-49)

Ga = gNa .

(A-50)

0

τ

with

Applying Eqs. (A-23)–(A-25), (3), we find from Eq. (A-49) that

 = − I + Ge (V Be − I ) + Ga (F · S · F − I ), where

 S=

∞ 0

f (u )s(t, u )du

(A-51)

(A-52)

332

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

with

 

s(t, u ) = exp −

t 0

 t  t  (s, u )ds I + (τ , u ) exp − (s, u )ds C−1 (τ )dτ , τ

0

(A-53)

and

C = F · F

(A-54)

is the right Cauchy–Green tensor for macro-deformation. Differentiation of Eq. (A-53) with respect to time implies that the tensor function s(t, u) obeys the equation

  ∂s (t, u ) = (t , u ) C−1 (t ) − s(t, u ) , ∂t

s ( 0, u ) = I.

(A-55)

A4. Uniaxial tension Under uniaxial tension, the deformation gradient for macro-deformation reads

1 F = ki1 i1 + √ (i2 i2 + i3 i3 ), k

(A-56)

where k stands for elongation ratio, k(0 ) = 1, and im (m = 1, 2, 3) are unit vectors of a Cartesian frame. Eqs. (A-1), (A-9) and (A-56) imply that

  1 k˙ i1 i1 − (i2 i2 + i3 i3 ) . k 2

D=

(A-57)

We search the deformation gradient for plastic deformation Fp is the form similar to Eq. (A-56),

F p = k p i1 i1 +

1

 (i2 i2 + i3 i3 ),

(A-58)

kp

where the unknown function kp (t) obeys the initial condition kp (0 ) = 1 (this equality ensures that Eq. (A-7) is satisfied). Combination of Eqs. (A-3), (A-4) and (A-58) yields

dp =

  k˙ p 1 i1 i1 − (i2 i2 + i3 i3 ) . kp 2

(A-59)

We search the deformation gradient for elastic deformation in the form

Fe = ke i1 i1 +

1

 (i2 i2 + i3 i3 ),

(A-60)

ke

where the function ke (t) obeys the initial condition ke (0 ) = 1. Substitution of expressions (A-58), (A-60) into Eqs. (A-11), (A-16), respectively, results in









bp −

1 1 2 2 Ip1 I = k − 3 3 p kp

Ce −

1 1 2 2 Ie1 I = k − 3 3 e ke



i1 i1 −

1 (i2 i2 + i3 i3 ) , 2

i1 i1 −

1 (i2 i2 + i3 i3 ) . 2



Inserting these expressions and Eq. (A-59) into Eq. (A-45), we arrive at the differential equation





k3p − 1 k˙ p k3 − 1 =P V e −R , kp ke kp

(A-61)

where

P= and

1 PGe , 3



V = 1−



1 2 2 k + −3 K e ke

−1

is given by Eq. (A-46). Combination of Eqs. (A-16) and (A-60) yields

C˙ e = 2

  k˙ e 2 1 ke i1 i1 − (i2 i2 + i3 i3 ) . ke 2ke

(A-62)

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

333

It follows from Eqs. (A-57) and (A-60) that

F e · D · Fe =

  1 k˙ 2 ke i1 i1 − (i2 i2 + i3 i3 ) . k 2ke

Eqs. (A-11), (A-16), (A-58) and (A-60) imply that





1 1 R Ce · bp + bp · Ce + (V Ie1 − RIp1 )Ce − V C2e 2 3   3 3 k − 1  2 2 1 k −1 p =− V e −R k e i1 i1 − (i2 i2 + i3 i3 ) . 3 ke kp 2ke Inserting these expressions into Eq. (A-48), we arrive at the equation

 k3 − 1 k3p − 1  k˙ e k˙ = −P V e −R . ke k ke kp

(A-63)

Eqs. (8) and (9) follow from Eqs. (A-61) and (A-63). We search the tensor function s(t, u) in the form

s = s1 i1 i1 + s2 (i2 i2 + i3 i3 ),

(A-64)

where the functions sm (t, u) (m = 1, 2) satisfy the condition s1 (0, u ) = s2 (0, u ) = 1. Keeping in mind that

C−1 = k−2 i1 i1 + k(i2 i2 + i3 i3 ), and inserting Eq. (A-64) into Eq. (A-55), we arrive at Eqs. (11). It follows from Eqs. (A-52) and (A-64) that

S = S1 i1 i1 + S2 (i2 i2 + i3 i3 ), where

 S1 =

∞ 0

f (u )s1 (t, u )du,

(A-65)  S2 =

∞ 0

f (u )s2 (t, u )du.

These relations coincide with Eqs. (10). Substitution of Eqs. (A-56), (A-60) and (A-65) into Eq. (A-51) implies that the Cauchy stress tensor  reads

 = 1 i1 i1 + 2 (i2 i2 + i3 i3 ),

(A-66)

where

1 = − + Ge (V k2e − 1 ) + Ga (S1 k2 − 1 ),

2 = − + Ge

V ke



− 1 + Ga

S

2

k



−1 .

It follows from equilibrium equation (A-44) and the boundary condition in stresses at the lateral surface of a sample that

2 = 0 . Combination of these equalities implies that

    1 S2 1 = GeV k2e − + Ga S1 k2 − . ke

k

Introducing the engineering tensile stress

σ=

1 k

,

we arrive at Eq. (6) with

G = Ge + Ga ,

κ=

Ga . Ge + Ga

(A-67)

References Bacca, M., Creton, C., & McMeeking, R. M. (2017). A model for the Mullins effect in multinetwork elastomers. Journal of Applied Mechanics, 84, 121009. Bai, R., Yang, J., Morelle, X. P., Yang, C., & Suo, Z. (2018). Fatigue fracture of self-recovery hydrogels. ACS Macro Letters, 7, 312–317. Bai, R., Yang, Q., Tang, J., Morelle, X. P., Vlassak, J., & Suo, Z. (2017). Fatigue fracture of tough hydrogels. Extreme Mechanics Letters, 15, 91–96. Bartolome, L., Aurrekoetxea, J., Urchegui, M. A., & Tato, W. (2013). The influences of deformation state and experimental conditions on inelastic behaviour of an extruded thermoplastic polyurethane elastomer. Materials and Design, 49, 974–980. Beniah, G., Fortman, D. J., Heath, W. H., Dichtel, W. R., & Torkelson, J. M. (2017). Non-isocyanate polyurethane thermoplastic elastomer: Amide-based chain extender yields enhanced nanophase separation and properties in polyhydroxyurethane. Macromolecules, 50, 4425–4434. Buckley, C. P., Prisacariu, C., & Martin, C. (2010). Elasticity and inelasticity of thermoplastic polyurethane elastomers: Sensitivity to chemical and physical structure. Polymer, 51, 3213–3224. Chen, Y., & Guan, Z. (2014). Multivalent hydrogen bonding block copolymers self-assemble into strong and tough self-healing materials. Chemical Communications, 50, 10868–10870.

334

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

Chen, Y., Kushner, A. M., Williams, G. A., & Guan, Z. (2012). Multiphase design of autonomic self-healing thermoplastic elastomers. Nature Chemistry, 4, 467–472. Cho, H., Mayer, S., Poselt, E., Susoff, M., in’t Veld, P. J., Rutledge, G. C., & Boyce, M. C. (2017). Deformation mechanisms of thermoplastic elastomers: Stress-strain behavior and constitutive modeling. Polymer, 128, 87–99. Cho, H., Rinaldi, R. G., & Boyce, M. C. (2013). Constitutive modeling of the rate-dependent resilient and dissipative large deformation behavior of a segmented copolymer polyurea. Soft Matter, 9, 6319–6330. Clough, J. M., Creton, C., Craig, S. L., & Sijbesma, R. P. (2016). Covalent bond scission in the Mullins effect of a filled elastomer: Real-time visualization with mechanoluminescence. Advanced Functional Materials, 26, 9063–9074. Cordier, P., Tournilhac, F., Soulie-Ziakovic, C., & Leibler, L. (2008). Self-healing and thermoreversible rubber from supramolecular assembly. Nature, 451, 977–980. Dargazany, R., & Itskov, M. (2013). Constitutive modeling of the Mullins effect and cyclic stress softening in filled elastomers. Physical Review E, 88, 012602. De Espinosa, L. M., Fiore, G. L., Weder, C., Johan Foster, E., & Simon, Y. C. (2015). Healable supramolecular polymer solids. Progress in Polymer Science, 49–50, 60–78. Diani, J., Fayolle, B., & Gilormini, P. (2009). A review on the Mullins effect. European Polymer Journal, 45, 601–612. Drozdov, A. D. (2007). An unusual elastoplastic response of thermoplastic elastomers at cyclic deformation. International Journal of Engineering Science, 45, 660–678. Drozdov, A. D., & Christiansen, J. (2008). Cyclic thermo-viscoplasticity of carbon black-reinforced thermoplastic elastomers. Composites Science and Technology, 68, 3114–3122. Drozdov, A. D., & Gupta, R. K. (2003). Non-linear viscoelasticity and viscoplasticity of isotactic polypropylene. International Journal of Engineering Science, 41, 2335–2361. Drozdov, A. D., & Kalamkarov, A. L. (1996). A constitutive model for nonlinear viscoelastic behavior of polymers. Polymer Engineering and Science, 36, 1907–1919. Drozdov, A. D., & Yuan, Q. (2003). The viscoelastic and viscoplastic behavior of low-density polyethylene. International Journal of Solids and Structures, 40, 2321–2342. El Fray, M., & Altstadt, V. (2003). Fatigue behaviour of multiblock thermoplastic elastomers. 2. Dynamic creep of poly(aliphatic/aromatic-ester) copolymers. Polymer, 44, 4643–4650. Filippidi, E., Cristiani, T. R., Eisenbach, C. D., Waite, J. H., Israelachvili, J. N., Ahn, B. K., & Valentine, M. T. (2017). Toughening elastomers using mussel-inspired iron-catechol complexes. Science, 358, 502–505. Gent, A. N. (1996). A new constitutive relation for rubber. Rubber Chemistry and Technology, 69, 59–61. Gotz, C., Lim, G. T., Puskas, J. E., & Altstadt, V. (2014). The effect of carbon black reinforcement on the dynamic fatigue and creep of polyisobutylene-based biomaterials. Journal of the Mechanical Behavior of Biomedical Materials, 39, 355–365. Green, M. S., & Tobolsky, A. V. (1946). A new approach to the theory of relaxing polymeric media. Journal of Chemical Physics, 14, 80–92. Grindy, S. C., Learsch, R., Mozhdehi, D., Cheng, J., Barrett, D. G., Guan, Z., et al. (2015). Control of hierarchical polymer mechanics with bioinspired metal– coordination dynamics. Nature Materials, 14, 1210–1216. Guo, H., Guo, W., Amirkhizi, A. V., Zou, R., & Yuan, K. (2016). Experimental investigation and modeling of mechanical behaviors of polyurea over wide ranges of strain rates and temperatures. Polymer Testing, 53, 234–244. Herbst, F., Dohler, D., Michael, P., & Binder, W. H. (2013). Self-healing polymers via supramolecular forces. Macromolecular Rapid Communications, 34, 203–220. Higaki, Y., Suzuki, K., Kiyoshima, Y., Toda, T., Nishiura, M., Ohta, N., et al. (2017). Molecular aggregation states and physical properties of syndiotactic polystyrene/hydrogenated polyisoprene multiblock copolymers with crystalline hard domain. Macromolecules, 50, 6184–6191. Hui, C. Y., & Long, R. (2012). A constitutive model for the large deformation of a self-healing gel. Soft Matter, 8, 8209–8216. Jiang, F., Fang, C., Zhang, J., Wang, W., & Wang, Z. (2017). Triblock copolymer elastomers with enhanced mechanical properties synthesized by RAFT polymerization and subsequent quaternization through incorporation of a comonomer with imidazole groups of about 2.0 mass percentage. Macromolecules, 50, 6218–6226. Kajita, T., Noro, A., & Matsushita, Y. (2017). Design and properties of supramolecular elastomers. Polymer, 128, 297–310. Kang, J., Son, D., Wang, G. J. N., Liu, Y., Lopez, J., Kim, Y., et al. (2018). Tough and water-insensitive self-healing elastomer for robust electronic skin. Advanced Materials, 30, 1706846. Kim, S. M., Jeon, H., Shin, S. H., Park, S. A., Jegal, J., Hwang, S. Y., et al. (2018). Superior toughness and fast self-healing at room temperature engineered by transparent elastomers. Advanced Materials, 30, 1705145. Kumar, A., Francfort, G. A., & Lopez-Pamies, O. (2018). Fracture and healing of elastomers: A phase-transition theory and numerical implementation. Journal of the Mechanics and Physics of Solids, 112, 523–551. Lei, Z., & Wu, P. (2018). A supramolecular biomimetic skin combining a wide spectrum of mechanical properties and multiple sensory capabilities. Nature Communications, 9, 1134. Li, C. H., Wang, C., Keplinger, C., Zuo, J. L., Jin, L., Sun, Y., et al. (2016a). A highly stretchable autonomous self-healing elastomer. Nature Chemistry, 8, 618–624. Li, X., Stribeck, A., Schulz, I., Poselt, E., Eling, B., & Hoell, A. (2016b). Nanostructure of thermally aged thermoplastic polyurethane and its evolution under strain. European Polymer Journal, 81, 569–581. Li, Y., Tang, S., Kroger, M., & Liu, W. K. (2016c). Molecular simulation guided constitutive modeling on finite strain viscoelasticity of elastomers. Journal of the Mechanics and Physics of Solids, 88, 204–226. Liu, J., Liu, J., Wang, S., Huang, J., Wu, S., Tang, Z., et al. (2017a). An advanced elastomer with an unprecedented combination of excellent mechanical properties and high self-healing capability. Journal of Materials Chemistry A, 5, 25660–25671. Liu, J., Tan, C. S. Y., Yu, Z., Li, N., Abell, C., & Scherman, O. A. (2017b). Tough supramolecular polymer networks with extreme stretchability and fast room-temperature self-healing. Advanced Materials, 29, 1605325. Lu, T., Wang, J., Yang, R., & Wang, T. J. (2017). A constitutive model for soft materials incorporating viscoelasticity and Mullins effect. Journal of Applied Mechanics, 84, 021010. Luo, C. S., Wan, P., Yang, H., Shah, S. A. A., & Chen, X. (2017). Healable transparent electronic devices. Advanced Functional Materials, 27, 1606339. Mossi Idrissa, A. K., Wang, K., Ahzi, S., Patlazhan, S., & Remond, Y. (2018). A composite approach for modeling deformation behaviors of thermoplastic polyurethane considering soft-hard domains transformation. International Journal of Material Forming, 11, 381–388. Neal, J. A., Mozhdehi, D., & Guan, Z. (2015). Enhancing mechanical performance of a covalent self-healing material by sacrificial noncovalent bonds. Journal of the American Chemical Society, 137, 4846–4850. Pan, Y., & Zhong, Z. (2017). Modeling the Mullins effect of rubber-like materials. International Journal of Damage Mechanics, 26, 933–948. Parida, K., Kumar, V., Wang, J., Bhavanasi, V., Bendi, R., & Lee, P. S. (2017). Highly transparent, stretchable, and self-healing ionic-skin triboelectric nanogenerators for energy harvesting and touch applications. Advanced Materials, 29, 1702181. Plagge, J., & Kluppel, M. (2017). A physically based model of stress softening and hysteresis of filled rubber including rate- and temperature dependency. International Journal of Plasticity, 89, 173–196. Shen, F., Yuan, S., Guo, Y., Zhao, B., Bai, J., Qwamizadeh, M., et al. (2016). Energy absorption of thermoplastic polyurethane lattice structures via 3D printing: Modeling and prediction. International Journal of Applied Mechanics, 8, 1640 0 06. Sheridan, R. J., & Bowman, C. N. (2013). Understanding the process of healing of thermoreversible covalent adaptable networks. Polymer Chemistry, 4, 4974–4979.

A.D. Drozdov, J.d. Christiansen / International Journal of Engineering Science 133 (2018) 311–335

335

Stukalin, E. B., Cai, L. H., Kumar, N. A., Leibler, L., & Rubinstein, M. (2013). Self-healing of unentangled polymer networks with reversible bonds. Macromolecules, 46, 7525–7541. Tanaka, F., & Edwards, S. F. (1992). Viscoelastic properties of physically cross-linked networks. Transient network theory. Macromolecules, 25, 1516–1523. Tang, Z., Huang, J., Guo, B., Zhang, L., & Liu, F. (2016). Bioinspired engineering of sacrificial metal-ligand bonds into elastomers with supramechanical performance and adaptive recovery. Macromolecules, 49, 1781–1789. Tomita, S., Wataoka, I., Igarashi, N., Shimizu, N., Takagi, H., Sasaki, S., & Sakurai, S. (2017). Strain-induced deformation of glassy spherical microdomains in elastomeric triblock copolymer films: Time-resolved 2d-SAXS measurements under stretched state. Macromolecules, 50, 3404–3410. Vernerey, F. J. (2018). Transient response of nonlinear polymer networks: A kinetic theory. Journal of the Mechanics and Physics of Solids, 115, 230–247. Voyiadjis, G. Z., Shojaei, A., & Li, G. (2011). A thermodynamic consistent damage and healing model for self healing materials. International Journal of Plasticity, 27, 1025–1044. Wang, M., Shan, D., Liao, Y., & Xia, L. (2018). Investigation of inelastic behavior of elastomeric composites during loading-unloading cycles. Polymer Bulletin, 75, 561–568. Wang, Q., Gao, Z., & Yu, K. (2017a). Interfacial self-healing of nanocomposite hydrogels: Theory and experiment. Journal of the Mechanics and Physics of Solids, 109, 288–306. Wang, S., & Chester, S. A. (2018). Experimental characterization and continuum modeling of inelasticity in filled rubber-like materials. International Journal of Solids and Structures, 136–137, 125–136. Wang, W., Chen, M., Niu, Y., Tao, Q., Bai, L., Hou Chen, H., & Cheng, Z. (2017b). Facile one-pot synthesis and self-healing properties of tetrazole-based metallopolymers in the presence of iron salts. RSC Advances, 7, 47316–47323. Wu, J., Cai, L. H., & Weitz, D. A. (2017a). Tough self-healing elastomers by molecular enforced integration of covalent and reversible networks. Advanced Materials, 29, 1702616. Wu, S., Qiu, M., Tang, Z., Liu, J., & Guo, B. (2017b). Carbon nanodots as high-functionality cross-linkers for bioinspired engineering of multiple sacrificial units toward strong yet tough elastomers. Macromolecules, 50, 3244–3253. Xu, J. H., Ye, S., Ding, C. D., Tan, L. H., & Fu, J. J. (2018). Autonomous self-healing supramolecular elastomer reinforced and toughened by graphitic carbon nitride nanosheets tailored for smart anticorrosion coating applications. Journal of Materials Chemistry A, 6, 5887–5898. Yan, X., Liu, Z., Zhang, Q., Lopez, J., Wang, H., Wu, H. C., et al. (2018). Quadruple h-bonding cross-linked supramolecular polymeric materials as substrates for stretchable, antitearing, and self-healable thin film electrodes. Journal of the American Chemical Society, 140, 5280–5289. Yanagisawa, Y., Nan, Y., Okuro, K., & Aida, T. (2018). Mechanically robust, readily repairable polymers via tailored noncovalent cross-linking. Science, 359, 72–76. Yang, Y., Wang, X., Yang, F., Wang, L., & Wu, D. (2018). Highly elastic and ultratough hybrid ionic-covalent hydrogels with tunable structures and mechanics. Advanced Materials, 30, 1707071. Yoshida, S., Ejima, H., & Yoshie, N. (2017). Tough elastomers with superior self-recoverability induced by bioinspired multiphase design. Advanced Functional Materials, 27, 1701670. Zhang, B., Zhang, P., Zhang, H., Yan, C., Zheng, Z., Wu, B., & Yu, Y. (2017). A transparent, highly stretchable, autonomous self-healing poly(dimethyl siloxane) elastomer. Macromolecular Rapid Communications, 38, 1700110. Zhang, W., Liu, X., Wang, J., Tang, J., Hua, J., Lu, T., & Suo, Z. (2018). Fatigue of double-network hydrogels. Engineering Fracture Mechanics, 187, 74–93. Zhang, X., Tang, Z., Guo, B., & Zhang, L. (2018a). Enabling design of advanced elastomer with bioinspired metal-oxygen coordination. ACS Applied Materials and Interfaces, 8, 32520–32527. Zhang, Z., Chen, Q., & Colby, R. H. (2018b). Dynamics of associative polymers. Soft Matter, 14, 2961–2977. Zhang, Z. P., Rong, M. Z., & Zhang, M. Q. (2018c). Polymer engineering based on reversible covalent chemistry: A promising innovative pathway towards new materials and new functionalities. Progress in Polymer Science, 80, 39–93. Zhao, Y., Ning, N., Hu, X., Li, Y., Chen, F., & Fu, Q. (2012). Processing temperature dependent mechanical response of a thermoplastic elastomer with low hard segment. Polymer, 53, 4310–4317. Zhou, J., Jiang, L., & Khayat, R. E. (2018). A micro–macro constitutive model for finite-deformation viscoelasticity of elastomers with nonlinear viscosity. Journal of the Mechanics and Physics of Solids, 110, 137–154.