A Transient Microsphere Model for Nonlinear Viscoelasticity in Dynamic Polymer Networks

The mechanical properties of soft polymeric materials feature an enormous range of both fluid-like and solid-like behaviors that are ultimately governed by interactions at the chain and network levels. Materials with supramolecular architecture such as physically bonded polymers are particularly diverse in their properties, despite similarities in the underlying structure. For instance, physical gels have been reported with different levels of viscoelasticity [1–3], self-healing [4,5], and hyperextensibility [6]. These properties can be attributed to the reversibility of physical bonds, which has also been shown to provide sacrificial energy dissipation and increase the toughness of hybrid (covalent and physical) materials [7]. This makes them desirable for a variety of applications, but there remains a large gap in our knowledge of the connection between chemical design and emergent properties. We thus need to gain a deeper understanding of the mechanisms that govern chain kinetics and, perhaps more importantly, how this drives the emergent network response.

While the chemical makeup of many physically bonded networks may be different, the physics that governs their behavior is fundamentally the same. Physical bonds have an activation energy close in magnitude to the thermal background, causing them to break and reform at relatively high rates. This type of reaction is well-described by kinetic theories such as those proposed by Eyring [8] and Kramers [9]. For instance, the rheological properties of some associating polymer solutions were well-described by relating the relaxation timescale of the material to the average lifetime of a bond [10,11]. This material template has remained common for studying the viscoelastic properties of physical networks, but discrepancies in the response are observed even for networks with similar architecture (see works by Erk et al. [12,13], Berret and Séréro [1,14], and Skrzeszewka et al. [3,15]). Further, while much effort has gone into studying the topology of these networks [16,17], there has been little work dedicated to directly characterizing the role of physical bonds. Indeed, a large variety of unique network architectures have been created [18], but their diversity hinders our ability to finely tune the desired mechanical properties. A primary goal of this study is thus to relate the governing physics of different networks to their resultant material behaviors.

Statistical approaches such as the transient network theory (TNT) [19,20] offer a promising template for bridging the connection between chain-level physics and material properties. The TNT is based on tracking the configuration of polymer chains within a network during deformation. With this approach, the behavior of a material can be fully described by knowing the distribution function $ϕ(r)$ of end-to-end vectors r of polymer chains in the network. While this description was capable of describing the complex mechanics of elastic damage [21] and the rheology of fire ant aggregates [22], it is computationally expensive to track the full distribution $ϕ(r)$ throughout deformation. To combat this, a reduced, tensorial form of the TNT was developed, which instead tracked the covariance μ of the distribution [20,23]. This has been applied to study a variety of transient behaviors such as biological growth [24], fracture [25], and the mechanics of topological materials [26], but the gain in numerical convenience comes at a loss of directional information. These reduced forms are thus insufficient for properly capturing complex chain-level interactions.

The difficulty in adding complexity to chain-level physics lies in the homogenization method used to bridge chain and network length scales. Many physically based formulations, such as Treloar’s three-chains model [27] and the Arruda–Boyce eight-chains model [28], make significant simplifications to the environment felt by a single chain to facilitate convenient constitutive laws. This comes at the cost of a loss in the directional properties of the network, however, which motivated the development of microsphere models that maintain a fully anisotropic description of the network [29–31]. This approach has been used to explore the effects of non-affine deformations [30,32] and provide physical interpretation to the Mullins effect in rubber [31,33]. Current viscoelastic microsphere models are restricted to rubber-like materials such as elastomers [34,35]. These models have done well to capture the time-dependency of topological effects originating from reptation and the release of chain restrictions (entanglements), but the dissipative mechanisms in dynamic polymer systems are quite different [20]. As such, behaviors unique to these materials such as stress overshoot [2,36], self-healing [4,5], and creep [10,14] cannot be described by rubber-like viscoelasticity. Development in this area could thus provide a vital connection between bond characterization and viscoelastic behavior in physically bonded networks.

The goal of this study is to unify experiments and theory in the context of dynamic polymer networks. We propose here a microsphere formulation of the TNT, also known as the transient microsphere model (TMM), which accounts for variable detachment rates. Our contributions are thus two-fold: (i) generalizing the microsphere approach to transient networks and (ii) further supplementing previous works on transient networks by exploring force-dependent bond kinetics. We demonstrate that this formulation is versatile and predicts a large variety of material responses using a minimal parametric space. In Sec. 2, we present the theoretical framework of the model and derive thermodynamically consistent constitutive laws. Section 3 then outlines our solution procedure and is dedicated to exploring the predictions of our model, which begins with defining its limiting cases, and is followed by direct comparisons to published experiments. We finish the paper with concluding remarks and an overview of the limitations and key findings of the model.

This section presents a theoretical basis for the TMM. We begin the detailed analysis of our system by considering flexible polymer chains that detach from the network and reattach at given rates. Expanding on our previous works, which have assumed constant rates of attachment and detachment [25,37], we derive a microsphere theory for dynamic networks which incorporates a force-dependent kinetic rate of detachment.

The full transient network theory approach describes the kinematics of a polymer chain by following the distribution P(r) of its end-to-end vector r under external loading. While this provides a wealth of information regarding the chain-level statistics, it can become computationally expensive for the study of macroscopic problems. Conversely, the fully reduced TNT proposed in Refs. [20,37] track only the statistical moments of the distribution and can be solved at the continuum level. They can thus be incorporated into larger-scale simulations using the finite element method [25,42], but the resulting model is fully isotropic. We therefore propose a semi-reduced model in the context of the microsphere formulations originally introduced by Treloar [43] and formalized by Wu and Van der Geissen [29,44]. In this approach, physical quantities of the network are replaced by a single representative chain in each direction, parameterized by their polar and azimuthal angles, $ϕ$ and $θ$⁠, respectively. To fully characterize the distribution, we require knowledge of two physical states of the network: the directional chain density $ρ(θ,ϕ)$ and the directional chain stretch $λ¯(θ,ϕ)$⁠.

The model presented in the previous section possesses a limited set of intrinsic variables as displayed in Table 1. In this section, we demonstrate that a large number of the reported properties of dynamic networks can be predicted using this limited set of physical parameters. We first briefly outline our discretization and solution procedure to provide a foundation for the illustrations to be presented. Following this, we compare the TMM to earlier hyperelastic and viscoelastic models. Finally, we outline the key findings of our study by identifying the response regimes exhibited by the model.

The finite domain of the microsphere facilitates a reduced numerical procedure without compromising the anisotropic properties of the network. To solve the differential equations presented in the previous section, we first discretize the microsphere on a 200 × 200 grid for $−π/2≤θ≤π/2$ and $0≤ϕ≤π/2$⁠, noting that symmetry allows us to only consider half of the unit sphere. With the initial conditions $ρ=ρ0$ and μ = μ₀, we use a backward-Euler integration scheme to track the chain conformation given a specified velocity gradient L. A typical iteration thus consists of (i) evolving $ρ$ and μ at the current time, (ii) solving for the Cauchy stress and energy dissipation by trapezoidal integration, and (iii) determining $ρ˙$ and $μ˙$ for the next time-step. We note here that, unless stated otherwise, the initial concentration of attached chains c₀ is maintained at 75% of the total population (c₀ = 0.75c_t) and the chain length is kept at N = 10 for each illustration in this section. The Cauchy stress is normalized by the shear modulus c₀k_bT in all plots for generality.

To demonstrate some of the simple predictions of the TMM, we first investigate the response of a single network consisting of bonds that detach according to Eq. (11) and reattach at a constant rate k_a. We take f₀ to be very large (f₀ > 100), so that the detachment rate is approximately constant for small loadings (Fig. 1(b)). The network is subjected to simple shear in the x₁ direction at a rate much greater than its resting kinetics (W = 100). We load to a shear of F₁₂ = 1, hold at for a short period of time before loading in uniaxial tension at the same Weissenberg number to F₁₁ = 1.5, and then hold at this stretch indefinitely (Fig. 3(a)). The decrease in F₁₂ during the uniaxial loading phase is due to the incompressibility of the material (requiring $detF=1$⁠). The corresponding stress versus time curve displays nonlinear stiffening behavior in the loading regimes (I–II and III–IV) and Maxwellian relaxation in the hold regimes (II–III and IV–V).

To illustrate the microstructural evolution predicted by the TMM, we plot in Fig. 3(b) the directional density $ρ$ and average chain stretch $λ¯$ within the plane of loading (⁠$ϕ=π/2$ for uniaxial tension in the x₁ direction). Here, the length of each line in a direction $θ$ corresponds to the average chain stretch $λ¯$⁠, while the color is mapped to the directional chain density $ρ$ in that direction. The top-left panel shows that the network is initially isotropic, as each line is the same length and color. During loading, the model predicts an increase in both the average stretch and density of chains in the principal loading direction (⁠$θ=π/4$⁠). A partial recovery is observed at time III, which is due to chain detachment from their stretched state and reattachment in their initial configuration. The network maintains a small amount of shear when loading in uniaxial tension as observed in panel IV. We predict full stress relaxation in the limit of t → ∞, as indicated by the V in the top-left panel. The illustrations depicted in Fig. 3(b) will be used for the remainder of the paper, but the scale and colorbars will be omitted to avoid crowding the figures. All future illustrations maintain the same scale system established in this figure.

The behavior of physical networks can exhibit a range of responses from that of elastic solids at high strain rates (W ≫ 1), to that of viscous fluids at low rates (W ≪ 1) [20]. In the context of the proposed model, the limit of large f₀ (f₀ → ∞) corresponds to regimes where bond dynamics are independent of force (k_d remains constant). The time dependence of the network is therefore fully decoupled from its deformation history. In these conditions, we compare the predictions of the TMM to earlier formulations under cyclic loading conditions and at various rates.

Elastic response. If the rate of loading is large (W ≫ 1), chains do not have time to detach during loading, and the network behaves elastically. There is no energy dissipation, and the curves in their loading and unloading regimes look identical (Fig. 4(a)). The TMM predicts severe strain stiffening around $λ=1/N$ as chains align and reach their contour length along the loading direction. In fact, at high strain rates, the TMM converges to the elastic microsphere models originally proposed by Treloar [43,46], for which generalized solutions were derived by Wu and Van der Giessen [29,44]. Thus, the TMM may be thought of as an extension of these earlier models with incorporated chain kinetics. The microsphere formulation predicts more stiffening than isotropic stretch models, such as the eight-chains model of Arruda and Boyce [28], which average the chain-stretch over all directions. It is worth noting that isotropic stretch models have actually been shown to outperform Microsphere formulations under certain loading conditions (for instance, when considering both uniaxial and biaxial response simultaneously). This point is further addressed in Sec. 4.

Linear viscoelastic response. Let us now investigate the prediction of the TMM for constant bond dynamics, but moderate to low strain rates (i.e., W ≤ 1). In these regimes, the network exhibits a combination of creep and elasticity that can be measured by hysteresis loops during cyclic loading. As shown in Fig. 4(b), the form of the stress-strain curve depends heavily on the loading rate. Previous models for capturing the effects of nonlinear elasticity and constant bond dynamics include the reduced transient network theory [37]. In this earlier formulation, the average elastic deformation of the chain was described by the conformation tensor μ, which is allowed to relax over time at rate k_d. Using an Arruda-Boyce like approximation for the network elasticity, the model predicted hysteresis loops shown as dotted lines in Fig. 4(b). The predictions of the TMM show notable discrepancies to the reduced transient network theory when the strain rate is large (W = 1 in the figure), but a convergence at very low strain rates (W = 0.5 in the figure). Similar to the elastic comparisons of Fig. 4(a), we observe that the TMM predicts more stiffening of the network at larger strains, and this effect becomes predominant as the strain rate increases. This observation can be explained by the degree of anisotropy developed in elastic networks; at high loading rates, chains tend to align and deform elastically in the stretch directions, thereby inducing a large proportion of highly stretched chains along the principal axis. This phenomenon, responsible for network stiffening, is not well-captured by the “average” tensorial form of the reduced Transient Network Theory. At low strain rates however, chains have time to detach from the network before becoming severely deformed. The network therefore remains more isotropic, with less elastic deformation—a regime that can be captured by an average tensorial formulation. To summarize, the TMM better reflects the anisotropic properties of the network either for elastic networks or at higher loading rates (W ≥ 0.5).

To set a reference, we first show in Fig. 5, the stress response of a network with constant bond dynamics (f₀ → ∞) at varying loading rates. It can be shown that when bond dynamics are decoupled from mechanics, the overall chain concentration remains constant and equal to its base value c₀ (determined in Eq. (13)) regardless of deformation history. In this case, the viscoelastic response—which is due to bond exchange only, but with no change in chain concentration—exhibits a stiffening regime that becomes less and less noticeable as W decreases. This is accompanied with an increase in viscous dissipation (due to bond exchange) that can visually be observed by the area enclosed within the hysteresis loops in Fig. 5. For these networks, the TMM predicts that the evolution of the local chain stretch $λ¯p$ and density $ρp$ closely follow their macroscopic counterparts; the directional stretch linearly increases with the macroscopic deformation, while the directional density (or chain alignment) exhibits similar hysteresis loop as the Cauchy stress. This second trend results from a lag in the chain alignment (compared to the applied strain) due to bond dynamics.

To illustrate the effect of force-sensitive bond detachment, let us now consider the fastest loading rate (W = 10) in Fig. 5 and explore the network response as f₀ varies from 1 to 10 (Fig. 6). For the largest value of f₀, a similar curve is recovered, with very little hysteresis loop. As bond kinetics become more sensitive to force, however, a clear trend of softening is observed, with a more pronounced hysteresis. This is explained by considering the increasing number of highly stretched chains with deformation; as this number increases, detachment events begin to overwhelm reattachment events due to accelerated bond kinetics. Indeed, the induced drop in the total concentration of attached chains, c, is most prevalent for the most sensitive chains (lowest f₀). We note that the change in c is a temporary consequence of the deformation and the TMM predicts full recovery of the chain concentration after unloading. The mechanisms by which chains stretch and detach are readily observed in the evolution of the directional stretch $λ¯p$ and chain density $ρp$⁠. We see here that at low deformations, the chains tend to realign with stretch, but this realignment is accompanied with a fast detachment of the most elongated chains. As the second mechanism prevails, we observe a drop in both the density of chains along the principal axis and their average stretch. These events together are responsible for the noticeable stress softening at the network level. As expected, this softening initiates earlier if the force sensitivity of the bond is large (i.e., f₀ is low). Thus, the chain sensitivity can be an important factor in activating energy dissipation, as observed by the increasing area within the hysteresis loop for lower f₀. This dissipation mechanism is however quite different from the networks with constant bond kinetics presented in Fig. 5, and this difference can be observed in the distinct signature of the stress-strain curves. In real polymers, energy dissipation may arise from a combination of these two processes, yielding more complex responses discussed in the next section.

We here present a detailed study of the TMM by considering the deformation of a network at a constant loading rate. Three dynamic processes govern the response of the network: (i) the elastic distortion of attached chains, (ii) the accelerated detachment of highly stretched chains, and (iii) the realignment of chains in the direction of loading. Each of these processes are coupled with each other: elastic deformation induces a faster detachment rate of highly stretched chains, chain alignment creates a higher concentration of chains detaching at a faster rate, and the reattachment of chains in their stress-free state temporarily resets the process in all directions for each increment of deformation. The evolution of the network thus consists of a transient regime, during which these processes equilibrate with the constant loading rate, and an equilibrium state, at which point the network topology has reached a state of dynamic equilibrium. The duration and difference between these two regimes are determined by the Weissenberg number W and the bond sensitivity f₀. In this section, we demonstrate a large variety of responses that can be captured by different combinations of these two parameters.

As the system equilibrates with loading, we expect the stress response to closely reflect the evolution of the microstructure. Indeed, the stress response predicted by the TMM typically consists of an initially increasing regime, followed by a creeping flow at large deformations. In Fig. 7(a), we illustrate the uniaxial tension response predicted by the TMM for networks with f₀ = 1, f₀ = 2, and f₀ = 5 at a loading rate of W = 1. In each case, a noticeable stress overshoot is observed, which marks the transition from the transient response to the equilibrium state. Increasing f₀ delays the acceleration of bond kinetics, creating a longer transient regime and a larger stress overshoot. We see that the point of overshoot corresponds with a sudden drop in $ρp$ and $λp$⁠, which indicates the existence of a critical stretch at which point the detachment rate of highly aligned chains overwhelms the loading rate. Indeed, the peak value of $λ¯p$ increases with increasing f₀, and the drop in $ρp$ becomes steeper. The secondary rise in $λp$ and $ρp$ is due to the stretching and alignment of fresh chains that have reattached to the network in a stress-free state after a detachment event. Figure 7(b) summarizes the stress overshoot predicted by the TMM as a function of W and f₀. Here, we normalize the overshoot $σo$ by the corresponding peak stress $σmax$ to promote an easier comparison. As expected, the lowest value of f₀ achieves the smallest relative overshoot, as the detachment rate of chains accelerates very quickly for this network. Increasing the Weissenberg number further delays yielding as elastic deformation is more influential. This creates a longer transient regime with a larger stress overshoot.

As we further increase the Weissenberg number, chains do not have time to detach from the network during loading, and there is a buildup of highly stretched chains. From Fig. 7, we expect this to be accompanied by a more delayed and extreme stress overshoot. In Fig. 8, we illustrate the same range of f₀ values as before, but at a higher Weissenberg number of W = 10. For small stretches, loading at this rate is fast enough to maintain a nearly elastic response, and each curve follows the same stiffening behavior. Once again, yielding coincides with a sudden decrease in $ρp$ and $λp$⁠, but this time, we also observe a markable drop in the total chain concentration c/c₀. This is explained by the higher degree of alignment at this loading rate; very few chains have detached during deformation, so the detachment event consists of a larger population of chains that have been stretched since the onset of loading. Once again, increasing f₀ delays the yielding and creates a higher peak stress.

This parametric study of the TMM reveals a large number of characteristic responses, which can potentially be used to physically interpret the behavior of previously studied physical networks. For instance, the fluid-like creep of associative polymer solutions is well-approximated by the TMM in the flow regime (left plots of Fig. 9). The bottom-left plot shows excellent agreement with shear startup experiments performed by Berret and Séréro [1,14] in the shear-thinning regime of their telechelic associative polymer solutions. Each curve reflects shear at a different loading rate, which varied between W = 0.2 and W = 1.8, with the highest stress being achieved at the highest loading rate. For this material, we determined a fitted value of f₀ (see Table 2I) based on the length scale of the backbone polyethylene oxide (PEO) [14]. We particularly note that the TMM accurately predicts both the transient and equilibrium response of the network. In contrast, the TMM underpredicts the stress overshoot observed in the faster experiments (2.5 ≤ W ≤ 5) done by Skrzeszewka et al. [36] for a similar family of telechelic polymers (top-left plot in Fig. 9). A larger overshoot induces a stress softening response that is known to be associated with strain localization. In this context, a higher loading rate induces the onset of shear banding within the material [1], which cannot be predicted by a local continuum model [48]. The polymers studied by Skrzeszewka et al. had a more complex collagen-like structure [3,36], which makes the determination of the Kuhn length b less well-defined. The normalized value of f₀ is, therefore, presented in Table 2 (II).

The failure regime of the TMM consists of two distinct responses characterized by the rate of loading. If the loading is moderate (W ≈ 1), the model predicts a response similar to unsteady creep; the stress initially softens due to the low loading rate, but begins to stiffen at higher stretches as some chains approach their nonlinear regime. This is therefore only observed for bonds with high f₀. The corresponding “S”-shaped curve was observed from the ionically bonded triblock copolymers studied by Henderson et al. [49] (bottom-right plot in Fig. 9). This network was similar in composition to the associative solutions studied previously, but the ionically modified end-groups creates an even more robustly connected network. Based on the length scale of the middle-group, poly(methacrylic acid) (PMAA), we fitted f₀ with a loading rate of W = 1.35 (Table 2(III)). It is interesting to note that the predictions of the TMM in this regime match closely with the rubbery viscoelasticity microsphere model of Miehe and Göktepe [34]. In this light, the dynamics of reptation and the release of entanglements in elastomers could be interpreted as the detachment of physical bonds between neighboring polymer chains. At larger Weissenberg numbers (W ≥ 10), the model predicts a purely stiffening response as discussed in Sec. 3.2. This type of rubber-like response was observed in the supramolecular network of Cordier et al. [6] (top-right plot in Fig. 9). The network consisted of a robust hydrogen-bonding network of multitopic molecules, which can be interpreted as scattered chains of clustered molecules held together by weak interactions. These clusters can be represented as freely jointed chains if there are minimal interactions between neighboring units on the same chain. The stiffening behavior and point of failure is well-captured by the TMM. We use a normalized value of f₀ for fitting (Table 2 (IV)) given the unorthodox chain topology.

In this study, we presented the transient microsphere model (TMM)—a microsphere formulation of the Transient Network Theory—to explore the role of force-sensitive bond kinetics in physical networks. The elastic response at high loading rates was compared to previous hyperelastic models, such as the full network theory of Wu and Van der Geissen, which the TMM converges to in the limit of large Weissenberg number. We also compared the TMM to previous reductions of the transient network theory and demonstrated a higher degree of stiffening, which is not well-approximated by the average stretch at large deformations. Upon incorporating force sensitivity into the bond kinetics, we observed two sources of energy dissipation in physical networks: (i) bond exchange without losses in the total chain concentration and (ii) network softening due to force-accelerated bond exchange. These processes yielded distinct regimes of response for simple loading scenarios such as uniaxial tension and simple shear, which were dependent on the rate of loading (W) and the bond sensitivity (f₀). When chain kinetics dominated the response of the material (low W and low f₀), the TMM predicts a stable stress-strain curve that ends in steady-state creep. If the material deforms more elastically (high W and high f₀), localized directions of sudden detachment events offered a metric for predicting failure. Using published experimental data, we demonstrated that each regime predicted by the model has been observed in physical polymer networks of varying composition. The TMM could therefore be used to quantify these networks based on physical parameters such as their rate of detachment and bond strength.

The presented model rests on several assumptions which could be broadened depending on the system being studied. In particular, it is known that the affine assumption breaks down at large deformations. As we are mainly interested in the viscoelastic properties of the network, we maintain the affine assumption and note that non-affine characteristics could be added later. We also note here that our formulation is only instantaneously affine with the macroscopic velocity gradient (hence, we write $r˙=Lr$⁠). Due to bond kinetics, the stretch felt by an individual chain does not reflect the global deformation gradient F except for the purely elastic case (W → ∞). Other characteristics of non-affine behavior have further been accounted for in a microsphere model by Miehe et al. [30]. It is also interesting to note that the stiffening behavior of rubbers is better approximated by the seemingly less-precise Arruda–Boyce model, which uses an isotropic approximation for the chain stretch [50]. The divergence of materials from their predicted alignment is likely also due to non-affine behaviors and topological constraints such as entanglements [30]. We believe that the added precision in modeling the anisotropic viscoelastic response sufficiently makes up for the loss of precision in this one regime.

A clear advantage to the microsphere approach is a finite numerical domain that still accounts for directional properties. Previous microsphere models [30,32] have further reduced their computational domain by using quadrature schemes such as those presented by Bazant and Oh [51]. This allows for the incorporation of fully anisotropic constitutive laws into more sophisticated studies using the finite element method [30]. This remains possible with the TMM, but the spatial derivatives in Eqs. (9) and (10) would be challenging to compute on a discrete quadrature grid. Future implementations of this model could surpass this by using an embedded finite element method such as Arbitrary Lagrangian–Eulerian, which uses an Eulerian material description in combination with a moving Lagrangian mesh. This model could then be used to predict more complex material processes such as fracture.

Physically bonded networks have been getting increased attention for their ability to strengthen materials via sacrificial bonding [7]. The desire to understand physical networks alone is thus vital to their success, but combining networks into hybrid materials is a primary application. Many methods of modeling these hybrid networks have been proposed, which could incorporate the effects of polydispersity [33] and various network interactions [30]. We note that it is possible to qualitatively predict these behaviors by summing networks in parallel (equal strain assumption) [21]. This could be done with the TMM to model hybrid materials with constituent networks ranging anywhere from purely elastic to linearly viscoelastic. On this note, the TMM, as presented in this work, is restricted to describing one characteristic timescale of a given network (determined by the rate k_d,0). Many real physical networks are however known to be better represented by a spectrum of timescales. These aspects may be easily incorporated into the TMM by considering multiple networks in parallel as discussed in prior works [20], an approach that is akin to a prony series decomposition. Another interesting feature of weak bonding interactions is their tendency to form clusters or a strong junction consisting of many weak interactions. The strength of these junctions varies depending on the number of contributing bonds [52]. This induces complex viscoelastic behaviors such as yielding and stress overshoot, which the TMM captures with the force sensitivity f₀. The added strength of a cluster of weak bonds can, therefore, be explained by a delayed acceleration of bond kinetics. From these observations, we believe that the insight provided by the TMM will be valuable for the rational design of materials for a desired function.

FJV gratefully acknowledges the support of the National Science Foundation under Award No. 1761918. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Science Foundation.

There are no conflicts of interest.

The authors attest that all data for this study are included in the paper.

Rearranging and reordering terms, the Cauchy stress and energy dissipation of Eqs. (16) and (17) are recovered.