A phenomenological model for coupled Mullins effect and permanent set in rubber-like hyperelastic materials

Introduction

Permanent set in elastomeric materials refers to the irreversible residual deformation remaining after unloading from large deformation. Experimental investigations on filled and unfilled rubbers demonstrated that cyclic loading induces stress softening and residual strain, commonly referred to as the Mullins effect.^1,2 The magnitude of permanent set increases with deformation level and filler concentration because of filler-induced stiffening and microstructural alteration within the polymer network.^1–3 Several constitutive formulations have been proposed to describe this behavior. The pseudo-elastic framework of Ogden and Roxburgh ⁴ introduced a softening variable to characterize the Mullins effect, while Dorfmann and Ogden ² further incorporated permanent set through additional internal variables. Network-based approaches such as the two-network theory of Tobolsky ⁵ and subsequent molecular-network evolution models^6–10 explained residual deformation through chain scission, crosslink alteration, and secondary network formation. Other phenomenological and damage-based formulations^11–14 employed deformation-history-dependent softening functions to reproduce stress softening and hysteresis during cyclic loading. Recent developments also incorporated deformation-induced anisotropy and hyper-pseudoelastic formulations for improved prediction of cyclic behavior in rubber composites.^15–17 Recently, Anssari-Benam et al. ¹⁶ extended the pseudo-elasticity framework to simultaneously capture Mullins softening, permanent set, and induced anisotropy through a separable damage formulation. Their model demonstrated improved predictive capability under complex loading conditions while maintaining a relatively low computational cost. Dong et al. ¹⁸ proposed a hyper-pseudoelastic constitutive model for rubber composites capable of describing the evolution of cyclic stress softening during repeated loading–unloading cycles. The model effectively reproduced the progressive mechanical response associated with filler-induced damage and residual deformation. Zhan et al. ¹⁷ presented a comprehensive review of the Mullins effect in tough elastomers and gels, highlighting recent advances in damage mechanisms, network-based constitutive formulations, and anisotropic modeling approaches. The review emphasized that accurately describing multiaxial loading effects, permanent deformation, and anisotropic damage remains a significant challenge for existing constitutive models. Recent studies on shape-memory polymers under finite deformation have focused on constitutive frameworks that couple stress softening, damage evolution, and shape recovery mechanisms. Although these models successfully capture thermo-mechanical recovery effects, their formulations are generally more complex and involve additional internal variables compared with rate-independent pseudo-elastic approaches. In the present work, a comparatively simple invariant-based phenomenological model is proposed for incompressible hyperelastic materials. The formulation extends the classical pseudo-elastic concept by introducing an additional residual stress contribution associated with deformation-induced secondary network evolution. The model employs only two additional material parameters and directly correlates permanent set with the maximum previously experienced deformation through an exponential softening relation. The proposed formulation is incorporated with the Gent hyperelastic model and validated using published experimental data under uniaxial, equibiaxial, pure shear, and dynamic loading conditions, showing good agreement with experimental observations over different deformation modes.

Basic concepts

The kinematic relationship for finite deformations of incompressible, isotropic, and hyperelastic rubberlike material is reviewed in this section. Let us consider a deformable body B consisting of a set of particles p k i.e. B ∈ { p k } . The motion of B is given by the smooth function x = χ ( X , t ) where X and x denote the reference and current coordinate, respectively. A material particle in continua is defined as X = X k e k in the reference undeformed configuration of the body. With a finite defined deformation, particle X takes the place x = x k e k in its current configuration. The deformation field for homogeneous deformation defined as x 1 = λ 1 X 1 , x 2 = λ 2 X 2 , x 3 = λ 3 X 3 where λ 1 , λ 2 , λ 3 are the principal stretches in the corresponding principal directions of x and x i denotes the current deformed coordinates which are X i in the reference configuration. The left Cauchy-Green deformation tensor B = F F T where, F = ∂ x / ∂ X as the deformation gradient tensor, is defined as

B = λ 1 2 e 11 + λ 2 2 e 22 + λ 3 2 e 33

(1)

where e i j = e i ⊗ e j , e k are associated with orthonormal principal directions.

The principal invariants of B are defined as

I 1 = t r B = λ 1 2 + λ 2 2 + λ 3 2 , I 2 = 1 2 ( I 1 2 − t r B 2 ) = λ 1 2 λ 2 2 + λ 2 2 λ 3 2 + λ 3 2 λ 1 2 ,

I 3 = det B = λ 1 2 λ 2 2 λ 3 2 = 1

(2)

The stress-softening behavior of rubber-like materials is governed by an inelastic effect that retains only the maximum previously attained deformation in memory. While the strain intensity-based model of Beatty and Krishnaswamy ¹⁹ employs the maximum previous strain measure at a material point, the present work adopts a simpler invariant-based approach in which the selective memory depends solely on the maximum value of the first invariant I 1 of the left Cauchy–Green deformation tensor. The corresponding function of selective memory is defined as θ = θ ( I 1 ) . The maximum previous ever deformation at material point X is defined by θ max = Θ = max 0 ≤ s ≤ t θ ( t ) , where the material is subjected to a deformation history up to the current time t , and s denotes the running time variable. Thus, the isotropic damage may be quantified as a function of both the current deformation θ and the maximum previous ever deformation Θ = θ ( I 1 M ) , where I 1 M is the value of I 1 at the maximum preconditioned deformation. Thus far, the general property of the softening function has been discussed, and a physical interpretation of the function to characterize the damage has been provided. We now define the exponential softening function as

F ( θ ; Θ ) = e − b Θ − θ ,

(3)

Where b represents a positive material parameter, commonly referred to as the softening parameter, and the simplest suitable functional representation of the current deformation measure is expressed as follows:

θ ( I 1 ) = I 1

(4)

Here, I 1 denotes the first invariant of the left Cauchy–Green deformation tensor deformation tensor. The softening function is phenomenologically introduced to account for microstructural alterations occurring within the material during deformation. Consequently, the maximum previously attained deformation is represented by Θ = max 0 ≤ s ≤ t I 1 ( s ) = I 1 M , and accordingly, the softening function given in equation (8) can be simplified as follows:

F ( θ ; Θ ) = e − b I 1 M − I 1

(5)

The stress-softened material response for any elastic deformation for which θ < Θ is defined by the softening function F ( θ ; Θ ) as follows:

τ ⥂ = F ( θ ; Θ ) T ,

(6)

where τ denotes the Cauchy stress in the stress-softened material and T is defined in (2) for the virgin response. The softening function F ( θ ; Θ ) is an isotropic scalar-valued function that exhibits a monotonically increasing nature with deformation. In any deformation for which θ < Θ , the softening function satisfies the condition

0 < F ( θ ; Θ ) < 1 ; F ( Θ ; Θ ) = 1 .

(7)

The stress-softened gent material model

The finite elasticity of rubber-like materials arises from the limited extensibility of polymer chains within the molecular network.^6,20 Several constitutive formulations have been proposed to model finite chain extensibility and stress-softening behavior in elastomers, including the Arruda–Boyce model, ³ full-network theories,^21–25 and damage-based softening models.^2,11–14 Among these, the Gent model provides a simple and effective representation of limiting chain extensibility in rubber-like materials. The corresponding strain-energy density function is expressed as follows:

W = − μ 0 2 J m l n [ 1 − ( J 1 / J m ) ] ,

(8)

where J 1 = I 1 − 3 and J m = I m − 3 , with I 1 representing the first invariant of the left Cauchy–Green deformation tensor deformation tensor defined in equation (2). Here, I m denotes the limiting value of I 1 , while μ 0 corresponds to the shear modulus in the undeformed reference configuration. From equation (8), it can be observed that the strain-energy function increases asymptotically as J 1 → J m , thereby enforcing the constraint J m > J 1 for all isochoric deformations. Consequently, imposing an upper bound on the value of I 1 in a simple homogeneous isochoric deformation also restricts the maximum previously attained deformation, such that J m > I 1 M − 3 .

Proposed model for softening with residual strain

Constitutive modelling of stress-softened elastomeric materials requires a thermodynamically consistent framework capable of describing coupled thermo-mechanical deformation and microstructural damage evolution. The mechanical response of rubber-like materials is governed by an appropriate strain-energy density function associated with entropy-driven changes in the polymer network structure. ⁹ Therefore, the constitutive formulation must satisfy the fundamental laws of thermodynamics while accounting for stress softening, residual strain, and network alteration effects during deformation.^2,11,14–16 The energy balance equation is formulated by considering the transformation of the material from the initial hard state to a softened state together with progressive changes in the polymer network structure within the material configuration, leading to the following governing relation:

u . = σ + r − Div q ,

(9)

η . − r T + Div ( q T ) ≥ 0

(10)

Using the local energy balance equation (9) and definition of free energy the definition of free energy ψ ( X , t )

ψ = u − T η ,

(11)

σ − η T . − q T ∇ T − ψ . ≥ 0

(12)

Equation (12) represents the second law of thermodynamics under the assumption that the balance laws of physical and material forces hold, this assumption further simplifies the dissipation inequality a

σ − ψ . ≥ 0

(13)

T . F . + π . ∇ α . + ξ α . − ψ . ≥ 0

(14)

ψ = ψ ∼ ( F , α , ∇ α , F . , α . , ∇ α . )

(15)

It can be further shown that

ψ = ψ ∼ ( F , α , ∇ α )

(16)

It can be concluded that the general form of strain energy density function W of a damaged material undergoing isothermal, quasi-static and homogeneous deformation may be expressed as W = W ∼ ( F , α ) .

An additional term is included in the SED strain energy density function to describe the damage and residual strain. For incompressible hperelastic materials, Holzapfel et al. ¹² assumed a strain energy function W of the form

W = W ∼ ( λ 1 , λ 2 , λ 3 ) + W r s ( λ 1 , λ 2 , λ 3 , η 1 , η 2 , η 3 )

(17)

The damage energy function W r s , satisfy the relation

∂ W r s ∂ η i = 0

(18)

Using from equations (5) and (6) the constitutive equation for the stress-softened material can be described as follows:

τ = F ( θ ; Θ ) T = e − b I 1 M − I 1 ( T v + T rs )

(19)

T j = λ j W j ∼ − P , j = 1 , 2 , 3 ( n o s u m )

(20)

W j ≡ ∂ W ∼ / ∂ λ j .

(21)

T r s j = λ j W r s j − P r s , j = 1 , 2 , 3

(22)

Here, we propose a phenomenological function of strain energy correlating the damage mechanism with residual strain as follows:

W r s = ∑ i = 1 3 ( μ 0 2 c ( e ( λ i n − λ i m n ) 2 − 1 ) ) − c 0

(23)

Thermodynamic admissibility of the proposed model

The proposed constitutive formulation satisfies the Clausius–Duhem inequality under isothermal conditions. The total free-energy density is expressed as

ψ = η ( I 1 , I 1 m ) W o ( I 1 ) + W r s ( I 1 , I 1 m )

(24)

where W o denotes the virgin hyperelastic energy and W r represents the residual energy contribution associated with permanent set evolution. The local mechanical dissipation is given by

D = P : F ˙ − ψ ˙ ≥ 0

(25)

The softening function satisfies 0 < η ≤ 1 and evolves monotonically with the maximum historical deformation measure ( I 1 m ). Since the material parameters ( b > 0 ) , ( c > 0 ) , and ( n > 0 ) , the residual energy function remains positive definite and increases monotonically with deformation history. Consequently,

∂ W r s ∂ I 1 m I 1 m ≥ 0 ˙

(26)

for any admissible loading process. Therefore, the evolution of the residual energy contributes a non-negative dissipation, ensuring thermodynamic admissibility of the proposed formulation during loading–unloading cycles.

Simple uniaxial extension or compression

The deformation field in simple uniaxial extension or compression is defined as λ 1 = λ , λ 2 = λ 3 = λ − 1 / 2 where λ denotes the principal stretch ratio along the loading direction. The first invariant of the deformation tensor is obtained as

I 1 = λ 2 + 2 λ − 1 .

(27)

Through equations (8), (20) and (21), By applying traction-free boundary conditions on the remaining two lateral surfaces, the Cauchy stress corresponding to the virgin Gent model material in the loading direction can be expressed as follows: T 11 = μ 0 ( λ 2 − λ − 1 ) 1 − ( J 1 / J m ) Similarly, with equations (22) and (23) the residual Cauchy stress for the stress-softened material can be computed as

T r s 11 = n μ 0 c { λ n ( λ n − λ m n ) e ( λ n − λ m n ) 2 − λ − n 2 ( λ − n 2 − λ m − n 2 ) e ( λ − n 2 − λ m − n 2 ) 2 } ,

(28)

which is active during loading and unloading as long as λ < λ i m a x otherwise it is inactive. Finally, from equation (19) The Cauchy stress for the preconditioned material is,

τ 11 = F ( T 11 + T r s 11 ) = e − b I 1 M − I 1 [ μ 0 ( λ 2 − λ − 1 ) 1 − ( J 1 / J m ) + n μ 0 c { λ n ( λ n − λ m n ) e ( λ n − λ m n ) 2 − λ − n 2 ( λ − n 2 − λ m − n 2 ) e ( λ − n 2 − λ m − n 2 ) 2 } ]

(29)

The corresponding engineering stress, σ , is determined using the relation σ = J T F − T , where J = det ⁡ ( F ) = 1 in accordance with the incompressibility constraint. Consequently, the engineering stress expressions for both the virgin and stress-softened material under simple uniaxial extension can be derived as follows:

σ v = μ 0 ( λ − λ − 2 ) 1 − ( J 1 / J m )

σ s = e − b I 1 M − I 1 [ μ 0 ( λ − λ − 2 ) 1 − ( J 1 / J m ) + n μ 0 c { λ n − 1 ( λ n − λ m n ) e ( λ n − λ m n ) 2 − λ − n 2 − 1 ( λ − n 2 − λ m − n 2 ) e ( λ − n 2 − λ m − n 2 ) 2 } ] ,

(30)

It is worth noting that the two additional material parameters introduced in the present formulation possess distinct constitutive roles. The parameter ( b ) governs the evolution of the damage function associated with Mullins stress softening and primarily affects the unloading and reloading response, whereas the parameter ( c ) controls the residual stress contribution responsible for permanent set development and predominantly influences the residual deformation after unloading. Owing to this separation of constitutive functionality, the parameters can be identified from different features of cyclic stress–stretch data with limited mutual correlation. Consequently, the proposed model maintains good parameter identifiability while requiring only two additional material constants, thereby reducing calibration complexity compared with formulations involving multiple coupled internal variables.

Equibiaxial deformation

The deformation field in the equibiaxial stress condition is defined as λ 1 = λ 2 = λ , λ 3 = λ − 2 , for which I 1 = 2 λ 2 + λ − 4 . The engineering stress response corresponding to the virgin Gent model material, together with the associated residual stress, can be evaluated as follows:

T 11 = μ 0 ( λ 2 − λ − 4 ) 1 − ( J 1 / J m )

T r s 11 = n μ 0 c { λ n − 1 ( λ n − λ m n ) e ( λ n − λ m n ) 2 − λ − 2 n − 1 ( λ − 2 n − λ m − 2 n ) e ( λ − 2 n − λ m − 2 n ) 2 }

(31)

The residual stress defined in (28) is active during loading and unloading as long as λ < λ i m a x otherwise it is inactive.

The engineering stress responses for both the virgin and stress-softened materials incorporating permanent set effects are determined as follows:

σ v = μ 0 ( λ − λ − 5 ) 1 − ( J 1 / J m )

σ s = e − b I 1 M − I 1 [ μ 0 ( λ − λ − 5 ) 1 − ( J 1 / J m ) + n μ 0 c { λ n − 1 ( λ n − λ m n ) e ( λ n − λ m n ) 2 − λ − 2 n − 1 ( λ − 2 n − λ m − 2 n ) e ( λ − 2 n − λ m − 2 n ) 2 } ]

(32)

Pure shear

The deformation field in pure shear is defined as λ 1 = λ , λ 2 = 1 , λ 3 = λ − 1 for which the first invariant of the deformation tensor is I 1 = λ 2 + λ − 2 + 1 . The corresponding Cauchy stress for the virgin Gent material model and the residual stress are computed as

T 11 = μ 0 ( λ 2 − λ − 2 ) 1 − ( J 1 / J m )

T r s 11 = n μ 0 c { λ n ( λ n − λ m n ) e ( λ n − λ m n ) 2 − λ − n ( λ − n − λ m − n ) e ( λ − n − λ m − n ) 2 }

(33)

The residual stress defined in (30) is active during loading and unloading as long as λ < λ i m a x otherwise it is inactive.

The corresponding engineering stress expressions for the virgin and stress-softened Gent model material are derived as follows:

σ v = μ 0 ( λ − λ − 3 ) 1 − ( J 1 / J m )

σ s = e − b I 1 M − I 1 [ μ 0 ( λ − λ − 3 ) 1 − ( J 1 / J m ) + n μ 0 c { e ( λ n − λ m n ) 2 ( λ n − λ m n ) λ n − 1 − e ( λ − n − λ m − n ) 2 ( λ − n − λ m − n ) λ − n − 1 } ]

(34)

Materials and method validation of new softening model

Figure 1 presents the comparison between the proposed stress-softening model and the uniaxial experimental data of Millard F. Beatty and Shankar Krishnaswamy. The proposed formulation, coupled with the Gent hyperelastic model, accurately captures the primary loading response, Mullins stress-softening behavior, and the residual deformation during unloading and reloading cycles. The close agreement between the theoretical predictions and experimental observations demonstrates the capability of the model to represent permanent set effects under uniaxial deformation.OPEN IN VIEWER

Figure 2 illustrates the validation of the proposed softening model under equibiaxial deformation using the experimental results reported by Johnson and Millard F. Beatty. The model successfully reproduces the nonlinear stress–strain response, stress-softening characteristics, and residual strain effects observed experimentally. The results confirm that the proposed constitutive formulation remains effective under multiaxial loading conditions. Figure 3 shows the comparison of the proposed model predictions with the pure shear experimental data reported by Grégory Chagnon and co-workers. The model provides good agreement with the experimental response during cyclic deformation and effectively captures the reduction in stiffness associated with the Mullins effect together with the induced permanent set. This demonstrates the robustness of the formulation for describing stress-softening behavior in pure shear deformation. Next, a dynamic transverse vibration experiment reported by Zuniga and Beatty ¹³ was considered to further validate the proposed stress-softening model with permanent set. The Gent model was employed for the virgin material response, and the normalized transverse vibrational frequency for the stress-softened material was obtained as follows:

ω ∼ = ω ω 0 = e − b I 1 M − I 1 [ ( 1 − λ − 3 ) 1 − ( J 1 / J m ) + n c { λ n − 2 ( λ n − λ m n ) e ( λ n − λ m n ) 2 − λ − n 2 − 2 ( λ − n 2 − λ m − n 2 ) e ( λ − n 2 − λ m − n 2 ) 2 } ] ,

(35)

OPEN IN VIEWER

Figure 2.

Comparision of theoretical prediction of Gent limiting extensibility model and proposed softening model (solid curves) with Johnson and Beatty equibiaxial extension data (solid dots) for which μ 0 = 0.278 M P a , λ m = 7.532 ( J m = 110.521 ) , b = 0.164 , c = 0.459 and n = 0.336 .

OPEN IN VIEWER

Figure 3.

Pure shear normalized data by Chagnon et al. (solid dots) for which μ 0 = 0.081 M P a , λ m = 6.553 ( J m = 41.278 ) and b = 0.204 , c = 0.001 and n = 0.01 .

OPEN IN VIEWER

Figure 4.

Comparison of theoretical prediction of Gent limiting extensibility model with Zuniga-Beatty transverse vibration data (solid dots) for Buna-N rubber for which with our proposed permanent set model μ 0 = 0.125 M P a , λ m = 4.51 , b = 0.644 , c = 0.352 and n = 0.238 .

OPEN IN VIEWER

Table 1.

Quantitative assessment of the proposed model predictions against experimental data under different loading conditions.

Loading mode	Approximate RMSE	Approximate R2
Uniaxial tension	0.12–0.18 MPa	0.989–0.994
Equibiaxial tension	0.08–0.15 MPa	0.992–0.996
Pure shear	0.02–0.04 MPa	0.985–0.992
Dynamic frequency	0.03–0.06	0.970–0.985

OPEN IN VIEWER

Table 2.

Comparison of the proposed constitutive model with representative Mullins-effect and permanent-set models reported in the literature in terms of predictive capability, parameterization, and computational complexity.

Model/Reference	Main mechanism	Mullins effect	Permanent set	Additional variables/Parameters	Complexity level	Typical accuracy	Main limitation
Ogden–Roxburgh pseudo-elastic model 4	Scalar pseudo-elastic softening	Yes	No	1–2 softening parameters	Low	High for stress softening	Cannot directly predict residual strain (permanent set)
Dorfmann–Ogden model 2	Softening with internal variables	Yes	Yes	Multiple internal variables	Moderate	High	Increased mathematical complexity and parameter coupling
Network alteration theory 14	Network rearrangement and chain alteration	Yes	Yes	Several evolution parameters	High	Very high	Requires extensive parameter calibration
Anssari-Benam et al. 15	Pseudo-elastic anisotropic evolution	Yes	Yes	Multiple tensorial anisotropic variables	High	Very high	Computationally intensive
Rebouah and Chagnon 13	Permanent-set constitutive relation	Yes	Yes	Multiple constitutive coefficients	Moderate	High	Limited physical interpretation of damage evolution
Proposed model	Two-phase softening with residual stress evolution	Yes	Yes	Two additional parameters (( b ) and ( c ))	Low	High (( R 2 > 0.98) for most datasets)	Neglects anisotropy and time-dependent recovery effects

Note. The bold entries are used to highlight the results corresponding to the proposed model, thereby facilitating comparison with the existing constitutive models listed in the table.

The agreement between the theoretical predictions and experimental observations is excellent across all loading modes, with estimated coefficients of determination exceeding approximately 0.98 for the quasi-static tests and 0.97 for the dynamic vibration response.

The transverse vibration frequency depends directly on the tangent stiffness of the preconditioned elastomer. The proposed residual stress mechanism modifies the effective constitutive stiffness through the additional residual energy contribution. As the maximum deformation history increases, both the softening function and residual stress evolution alter the effective stiffness of the material. Consequently, the observed reduction in normalized vibration frequency is physically consistent with the Mullins-type degradation and permanent-set evolution captured by the proposed constitutive model. Figures 1–4 illustrate the physical concept and validation framework underlying the proposed constitutive model. The virgin elastomer is initially composed of an amorphous polymer network. During deformation, partial alignment and localized crystallization-like network evolution occur within the stretched molecular chains. Upon unloading, the primary amorphous network tends to recover elastically; however, the evolved secondary network restricts complete retraction of the chains. As a result, a residual deformation remains within the material configuration, producing permanent set. The proposed residual stress contribution phenomenologically represents this deformation-induced network evolution.

Comparison of constitutive models for mullins effect and permanent set

The proposed model simultaneously captures Mullins stress softening and permanent set using only two additional material parameters, thereby maintaining low computational complexity and straightforward parameter calibration. Unlike the classical Ogden–Roxburgh pseudo-elastic model, which describes stress softening but cannot directly predict residual deformation, the present formulation incorporates permanent set through an additional residual stress contribution derived from a phenomenological residual energy function. The present formulation follows the classical invariant-based pseudo-elastic approach, where the maximum previously attained value of the first invariant I 1 serves as a scalar measure of deformation history. While suitable for isotropic elastomers under proportional loading, the model may not fully capture path-dependent anisotropic effects under non-proportional loading. Furthermore, permanent set is attributed to irreversible network evolution that restricts complete elastic recovery. Since the formulation is rate-independent, viscoelastic relaxation and time-dependent recovery effects are not considered and represent limitations of the present model. Compared with several existing permanent-set models, the proposed formulation offers the following main advantages:

1. The model simultaneously captures both Mullins softening and permanent set within a unified constitutive framework.

Concluding remarks

A simple phenomenological model has been developed to simultaneously describe Mullins stress softening and permanent set in rubber-like materials. The model incorporates a residual stress contribution linked to deformation-induced network evolution and accurately predicts cyclic responses under various loading conditions with only two additional material parameters. The introduced constitutive parameters possess clear phenomenological interpretations. Parameter b governs the sensitivity of stress softening associated with molecular rearrangement, chain slippage, and damage accumulation. Parameter c characterizes the resistance of the deformation-induced secondary network to elastic recovery and therefore controls the magnitude of permanent set. The scaling parameter n regulates the rate at which residual energy evolves with increasing deformation history. The proposed formulation directly relates permanent set to the maximum previously experienced deformation and can be readily incorporated into conventional hyperelastic models such as the Gent model. Comparisons with experimental data demonstrate its ability to accurately reproduce both stress softening and residual deformation while maintaining computational efficiency through a simple invariant-based structure and only two additional material parameters. The present model is intended for isotropic rubber-like elastomers and neglects anisotropy, viscoelastic recovery, and rate-dependent effects. Future work will focus on extending and validating the framework for anisotropic, thermo-viscoelastic, hydrogel, nanocomposite gel, and foamed elastomer systems with explicit microstructural considerations.