# Modeling of nonlinear viscoelasticity and stress
softening in soft tissues

Anne-Sophie Caro, André Chrysochoos, Sarah Iaquinta, Grégory Chagnon

<figure><img src="image_1.png" /></figure>

## To cite this version:

Anne-Sophie Caro, André Chrysochoos, Sarah Iaquinta, Grégory Chagnon. Modeling of nonlinear viscoelasticity and stress softening in soft tissues. European Journal of Mechanics - A/Solids, 2026, 115, pp.105818. 10.1016/j.euromechsol.2025.105818. hal-05229634

HAL Id: hal-05229634

https://imt-mines-ales.hal.science/hal-05229634v1

Submitted on 1 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Full length article

# Modeling of nonlinear viscoelasticity and stress softening in soft tissues

Anne-Sophie Caro<sup>a ID,*</sup>, André Chrysochoos<sup>b</sup>, Sarah Iaquinta<sup>a</sup>, Grégory Chagnon<sup>c</sup>

## ARTICLE INFO

### Keywords:

- Soft tissue

- Finite strain

- Visco-hyperelasticity

- Mullins effect

- Incompressibility

## ABSTRACT

This paper deals with the mechanical behavior of soft living tissues under load-unload and relaxation cyclic strains. It proposes a thermodynamic model formulated within the Generalized Standard Materials framework that incorporates both Mullins' effect and viscoelasticity, integrating the history dependent behavior of the material under finite strain. A key innovation lies in the use of a common softening function to modulate both the hyperelastic and viscous components, capturing history-dependent behavior more accurately. Viscous parameters are adjusted based on the loading history via the maximum strain invariant. Numerical implementation is validated against uniaxial tensile tests on porcine perineal tissues and a global Sobol sensitivity analysis confirms that elastic, viscous, and Mullins-related parameters are identifiable from different phases of the loading protocol. This model provides a unified, thermodynamical consistent tool for simulating soft tissue mechanics.

## 1. Introduction

The Mullins effect is a well-established phenomenon observed in filled elastomers and soft living tissues, characterized by stress softening upon cyclic loading and unloading. This phenomenon has been widely studied in the context of elastomers such as rubber-like materials (Mullins and Tobin, 1965; Mullins, 1969; Diani et al., 2009) and more recently, in collagenous tissues (Peña et al., 2011; Holzapfel and Ogden, 2020; Ortún-Terrazas et al., 2019; Rickaby and Scott, 2014; Wenguang, 2016). The Mullins effect is often attributed to irreversible microstructural changes such as the breakage of molecular bonds or the interaction between filler particles and polymer chains or even sliding and friction between polymer chains (Diani et al., 2009; Cantournet et al., 2009). However, despite decades of research, the precise mechanisms involved remain controversial, with several competing theories regarding its origin. These uncertainties complicate the development of accurate models, particularly for materials that undergo large deformations, where the microstructural rearrangements play a significant role. To account for the Mullins softening effect, Simo (1987) proposed a model that modifies the classic elastic strain energy by introducing a damage parameter, which evolves based on the deformation history of the material. Over the last few decades, several models have been developed with different softening models, even when the stress-strain response during reloading is the same as during unloading. Indeed, there are two main classes of models: those where the response of the material to reloading is identical to

unloading (ideal Mullins effect), and those where reloading shows a distinct response. The first class is simpler and more widely used, but often fails to capture experimental data, especially in cyclic loading scenarios. More flexible models, such as those using damage-dependent variables or two-phase material models, provide better fits but present their own theoretical and experimental challenges. In addition to these challenges, the viscoelasticity of soft materials adds further complexity to the modeling process. Elastomers and biological tissues are known for exhibiting both hyperelastic and viscous behaviors that are coupled with the rate of loading (strain rate) and the applied stress. This strain-rate dependence introduces additional complexities, especially in the case of large deformations. Several models have been developed to capture the complex mechanical response of soft tissues under large deformations. Classical viscoelastic models often rely on linear formulations or small-strain assumptions, limiting their applicability in finite deformation regimes. More advanced approaches, such as the fully nonlinear viscoelastic formulation of Simo (1987), introduce internal variables governed by time-dependent evolution laws and enabling the simulation of true stress relaxation.

In parallel, pseudo-elastic models have gained attention for their ability to capture strain-rate effects by allowing the parameters of hyperelastic potentials to vary with the applied strain rate (Anssari-Benam and Saccomandi, 2024). In such models, viscous effects are captured implicitly without introducing internal variables. These formulations

<footnote><sup>a</sup> LMGC, IMT Mines Alès, Univ Montpellier, CNRS, Alès, France</footnote>

<footnote><sup>b</sup> LMGC, Univ Montpellier, CNRS, Montpellier, France</footnote>

<footnote><sup>c</sup> Univ. Grenoble Alpes, CNRS, UMR 5525, VetAgro Sup, Grenoble INP, TIMC, 38000 Grenoble, France</footnote>

<footnote>* Corresponding author.
E-mail address: Anne-Sophie.Caro@mines-ales.fr (A.-S. Caro).</footnote>
provide efficient tools to describe monotonic and certain cyclic behaviors, especially under varying loading rates. More recently, some works have proposed unified pseudo-elastic formulations that capture both continuous softening during loading and discontinuous softening during unloading, without the need for explicit damage variables or history functions (Anssari-Benam and Hossain, 2023), Anssari-Benam and Hossain (2024), Anssari-Benam (2023). These approaches rely on hyperelastic models with rate-sensitive parameters that implicitly incorporate the effects of strain history and rate, offering a compact yet effective description of softening behaviors in elastomeric and biological materials under cyclic loading.

In contrast, the present study proposes a visco-hyperelastic model which explicitly separates elastic and viscous contributions: the elastic response is described both by a strain energy function and by the elastic elements within de viscoelastic branches, while the viscous behavior is described through internal variables evolving with time, enabling the simulation of stress relaxation. In addition, a history-dependent function acting on the elastic potential accounts for Mullins-type stress softening. This structure allows the model to capture the coupled effects of stress softening, and dissipation within a thermodynamically consistent framework. The formulation follows the Generalized Standard Materials (GSM) approach (Halphen and Nguyen, 1975), where both state and complementary laws are derived from energy and dissipation potentials.

The theoretical framework and the proposed thermomechanical model combining visco-hyperelasticity and the Mullins effect are presented in Section 2. Section 3 describes the experimental procedure on sow perineal tissues, the numerical implementation, and the global sensitivity analysis. Finally, Section 4 provides the discussion, conclusion, and perspectives for future work.

## 2. Stress softening and visco-hyperelastic model

### 2.1. Kinematics

Let $\mathbf{x} = \chi(\mathbf{x}_0, t) : \Omega_0 \times \mathbb{R}^+ \to \Omega$ the motion mapping and $\mathbf{F}$ the associated transformation gradient tensor. $\mathbf{x}_0$ and $\mathbf{x}$ are the positions of a particle in the (Lagrangian) reference $\Omega_0$ and (Eulerian) current $\Omega$ configurations such that $\mathbf{F} = \partial_{\mathbf{x}_0} \mathbf{x}$. Throughout this work, the notation $\partial_{\mathit{var}}$ denotes the partial derivative with respect to the variable $\mathit{var}$.

The Cauchy-Green C and Green-Lagrange strain tensors $\mathbf{E}$ corresponding to the above transformation gradient are written as:

$$ \mathbf{C} = \mathbf{F}^T \mathbf{F}, \quad \mathbf{E} = \frac{1}{2}(\mathbf{C} - \mathbf{1}), \qquad (1) $$

1 being the second order identity tensor. The stretch ratios $\lambda_i$, $i = 1..3$ correspond to the square roots of the eigenvalues of $\mathbf{C}$.

The invariants of $\mathbf{C}$, commonly used in large strain, can be expressed as:

$$ I_1 = \mathrm{tr}(\mathbf{C}), \quad I_2 = \frac{1}{2}[(\mathrm{tr}(\mathbf{C}))^2 - \mathrm{tr}(\mathbf{C}^2)], \quad I_3 = \det(\mathbf{C}) = \mathbf{J}^2, \qquad (2) $$

where $\mathbf{J} = \det \mathbf{F}$ is the Jacobian determinant of $\mathbf{F}$, representing the local volume change.

Some partial derivatives of these quantities with respect to $\mathbf{C}$ will be useful when presenting the model. Classically, it can be shown that:

$$ \partial_{\mathbf{C}} \det \mathbf{C} = \det \mathbf{C} \mathbf{C}^{-1}, \quad \partial_{\mathbf{C}} \mathbf{J} = \frac{1}{2} \mathbf{J} \mathbf{C}^{-1}. \qquad (3) $$

### 2.2. Experimental observation

The Mullins effect can be observed on filled elastomers and soft living tissues when they are tested through specific loading and unloading tests. These materials are thus subjected to repeated deformation cycles, with the load being increased and decreased at controlled strain rate. In the case of ideal Mullins effect, following the initial loading

<figure><img src="image_3.png" /><figcaption>Fig. 1. Mullins effect and viscoelasticity at finite strain.</figcaption></figure>

cycle (“active loading”), the material shows a decrease in stiffness and deviates from the original stress-strain path in subsequent cycles (“neutral loading”) while the maximum deformation is not reached. In the absence of viscosity, the stress-strain response during the second loading is identical to that during unloading (Mullins, 1969). In some rubber materials, the most commonly cited explanation for the Mullins effect involves the filler particles interacting with polymer chains or between polymer chains themselves via weak bonds. When first stressed, some of these bonds break irreversibly, leading to a loss of rigidity in the material. The weakest bonds break first, while those requiring greater energy break at higher elongation levels, needing a higher stress intensity (Diani et al., 2009). To account for this history effect, some authors (Chagnon et al., 2004; Rebouah and Chagnon, 2014) classically introduce a state variable represented by the maximum value of $\mathbf{C}$:

$$ I_1^*(t) = \max_{\tau \le t} (I_1(\tau)). \qquad (4) $$

This variable acts directly on the elastic energy stored by the material so that it induces a softening of the stress during unloading and then reloading, but also allows the material to recover its elastic energy level as soon as it returns to an elongation $I_1$ equivalent to $I_1^*$ which it had reached before unloading. It should be noted that during the active loadings $I_1^* = I_1$ and during the neutral loadings and unloadings $I_1^* = 0$.

In the case of visco-hyperelastic materials, due to viscosity, the neutral loading and unloading paths differ (see Fig. 1).

### 2.3. Thermomechanical framework

We consider an homogeneous and incompressible material.

The following model is consistent with those proposed by Simo (1987) and Simo and Hughes (1998) and all those who have followed his approach, such as Holzapfel and Gasser (2001) or Peña et al. (2011). Heuristically, viscoelasticity is approached by parallelizing a hyperelastic branch and one or more viscoelastic branches (generalized Maxwell-type model). Geometric and behavioral nonlinearities are carried by the hyperelastic branch, while viscous effects remain linear with respect to strain rate even at large strains (and reasonably low strain rates).

In the isovolumic model proposed below (see Fig. 2), the behavior is schematically represented by $N + 1$ parallel branches:

- a hyperelastic branch, which manages the nonlinear elastic behavior and the Mullins effect through a state variable that recapitulates the effects of the load history on the mechanical response,

- $N$ linear viscoelastic branches, whose viscosities $\mu^{(i)}$, $i = 1..N$ are also a function of the load history. This schematic formulation leads to a total stress response denoted by $\mathbf{S}_{\mathrm{iso}}$.

Although the schematic representation focuses on the isochoric part of the response (related to distortional deformation), the model also
<figure><img src="image_4.png" /><figcaption>Fig. 2. Rheological model of the isovolumic behavior: hyperelastic spring with Mullins effect in parallel with i, i = 1..N Maxwell viscoelastic branches.</figcaption></figure>

includes a volumetric contribution to enforce incompressibility, which is essential for soft biological tissues. This volumetric term is not illustrated in the graphical representation but is fully implemented in the numerical formulation.

### 2.3.1. States variables

In an isothermal context, the set of state variables retained are the overall Green-Lagrange strain tensor E, the N Lagrangian viscous stress tensor S<sub>v</sub><sup>(i)</sup> developed by the viscoelastic branch t, and the history scalar variable I<sub>1</sub><sup>*</sup> (defined in Eq. (4)), which allows the Mullins effect to be managed. The viscous stress tensor S<sub>v</sub><sup>(i)</sup> was chosen instead of a viscous deformation, as conventionally done for small strains, to avoid using any hypothesis regarding the decomposition of the transformation gradient into viscous and elastic parts, an attitude that is still often controversial.

### 2.3.2. Thermodynamic potential and state laws

In isothermal conditions, the Helmholtz free energy ψ of a conventional viscoelastic material corresponds to the portion of the deformation that is stored elastically. In the present context, this stored energy is fully described by the hyperelastic strain energy density W<sub>e</sub> as well as by the elastic components of the viscoelastic branches. The remaining part of the input energy is dissipated through viscous mechanisms.

Consequently, the thermodynamic potential chosen and denoted by G is the Legendre transform of this free energy ψ extended to large deformations. A relation between ψ and G can be given:

$$ \psi = G - \sum_{i=1}^{N} \partial_{S_{v}^{(i)}} G : S_{v}^{(i)}. \qquad (5) $$

Decomposing the potential G into isochoric and volumetric parts, we obtain

$$ G = G(E, S_v, I_1^*) = W_e(C, I_1^*) + S_v : E - \frac{1}{2} \sum_{i=1}^{N} \mathbb{R}^{(i)} : S_v^{(i)} : S_v^{(i)} + G_{vol}(J), \qquad \underbrace{G_{iso}(E, S_v, I_1^*)}_{G_{iso}(E, S_v, I_1^*)} \qquad (6) $$

where S<sub>v</sub> = Σ<sub>i=1</sub><sup>N</sup> S<sub>v</sub><sup>(i)</sup>. The last term in G<sub>iso</sub> involves a double contraction of a fourth-order tensor with two second-order tensors. In index notation, this reads as: R<sup>(i)</sup> : S<sub>v</sub><sup>(i)</sup> : S<sub>v</sub><sup>(i)</sup> = Σ<sub>i=1</sub><sup>3</sup> Σ<sub>j=1</sub><sup>3</sup> Σ<sub>k=1</sub><sup>3</sup> Σ<sub>ijkl</sub> R<sup>(i)</sup><sub>ijkl</sub> S<sub>v</sub><sup>(i)</sup> S<sub>ijkl</sub><sup>(i)</sup>. In (6) each term has the following meaning:

$$ \cdot W_e(C, I_1^*) $$

This term is the modified hyperelastic strain energy density, weighted by a softening function F(I<sub>1</sub>, I<sub>1</sub><sup>*</sup>), which is defined,

continuous and differentiable, to model the Mullins effect. Specifically,

$$ \partial_C W_e(E, I_1^*) = F(I_1, I_1^*) \partial_C W_e(C) \qquad (7) $$

where F(I<sub>1</sub>, I<sub>1</sub><sup>*</sup>) is a positive real smaller than 1 and satisfying:

$$ F(I_1, I_1^*) = 1, \partial_{I_1} F(I_1, I_1^*) = 0 \text{ and } \partial_{I_1}^* F(I_1, I_1^*) = 0 \text{ when } I_1 = I_1^*. \qquad (8) $$

The variable I<sub>1</sub><sup>*</sup> stores the maximum experienced value of I<sub>1</sub> with evolution defined by:

$$ I_1^* = \begin{cases} \dot{I}_1 & \text{if } I_1 = I_1^* \text{ and } \dot{I}_1 \ge 0 \\ 0 & \text{otherwise.} \end{cases} \qquad (9) $$

A possible choice for the softening function is:

$$ F(I_1, I_1^*) = 1 - \eta \left( \frac{I_1^* - I_1}{I_1^* - 3} \right)^\alpha, \qquad (10) $$

where η and α are adimensional material parameters.

A specific form of W<sub>e</sub> can be selected depending on the material behavior. The present framework is compatible with any hyperelastic model defined in terms of strain invariants. For a comprehensive overview of available constitutive models, the reader is referred to standard references on continuum mechanics and hyperelasticity (Marckmann and Verron, 2006).

$$ \cdot G_{vol} $$

The volumic part is defined as:

$$ G_{vol}(J) = p(J - 1), \qquad (11) $$

where p is a Lagrange multiplier enforcing the incompressibility condition J = 1.

$$ \cdot R^{(i)} $$

Each tensor R<sup>(i)</sup> represents a symmetric fourth-order compliance-like tensor and respects the incompressibility hypothesis. A typical form can be:

$$ R_{ijkl}^{(i)} = \frac{1}{2G^{(i)}} (\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) - \frac{1}{3G^{(i)}} \delta_{ij}\delta_{kl}, \quad i, j, k, l = 1..3, \quad i = 1..N \qquad (12) $$

where G<sup>(i)</sup> are shear-like moduli,

The state laws are the partial derivatives of the potential G with respect to the different state variables. We then get:

$$ \left\{ \begin{array}{l} \mathbf{A}_E = \partial_E G = \partial_E G_{iso} + \partial_E G_{vol} \\ \mathbf{A}_{S_v^{(i)}} = \partial_{S_v^{(i)}} G, \quad t = 1..N \\ A_{I_1^*} = \partial_{I_1^*} G. \end{array} \right. \qquad (13) $$

Taking Eq. (6) into account, the conjugate variables can be rewritten as:

$$ \left\{ \begin{array}{l} \mathbf{A}_E = 2\partial_C W_e(C, I_1^*) + S_v + pC^{-1} \\ \mathbf{A}_{S_v^{(i)}} = E - \mathbb{R}^{(i)} : S_v^{(i)}, \quad t = 1..N \\ A_{I_1^*} = \partial_{I_1^*} W_e(C, I_1^*) \end{array} \right. \qquad (14) $$

### 2.3.3. Dissipation potential and thermodynamic forces

As mentioned above, the GSM framework proposes to derive the complementary laws (i.e. thermodynamic forces) from a dissipation potential that depends on the fluxes of the selected state variables.

In what follows we suppose that the dissipative effects are exclusively due to viscosity, and that the viscous behavior remains linear. However, following the thermodynamic framework proposed by Halphen and Nguyen (1975), the dissipation potential is not only a function of the rates of internal variables, but may also depend explicitly on the state variables themselves, which are considered as parameters of the model:

$$ \varphi(\dot{\mathbf{E}}, \dot{\mathbf{S}}_v^{(i)}, \dot{I}_1^*; \mathbf{E}, \mathbf{S}_v^{(i)}, \mathbf{I}_1^*) = \frac{1}{2} \sum_{i=1}^{N} \mu^{(i)}(\mathbf{E}, \mathbf{I}_1^*) \dot{\mathbf{A}}_v^{(i)} : \dot{\mathbf{A}}_v^{(i)}, \qquad (15) $$
where $\dot{\mathbf{A}}_{\mathbf{S}_v}^{(i)} = \dot{\mathbf{E}} - \mathbf{R}^{(i)} : \dot{\mathbf{S}}_v^{(i)}$ and $\mu^{(i)}(\mathbf{E}, \mathbf{I}_1^*)$ is a viscosity coefficient that evolves with the material state. This formulation makes it possible to reflect, for instance, a viscosity that decreases during unloading, thus capturing the behavior of a wide range of soft materials under cyclic or complex loading. For instance, it may be modulated by the same function $\mathcal{F}$ used to describe the Mullins effect, although this formulation is not limited to such a dependency and can accommodate other forms of state evolution. As a result of this choice, the dissipation potential no longer explicitly depends on the internal variable flux $\dot{\mathbf{I}}_1^*$ nor on the viscous stress tensor $\mathbf{S}_v^{(i)}$.

The thermodynamic forces of the model play the same role with respect to the dissipation potential as the conjugate variables with respect to the thermodynamic potential. In this case the forces are defined by:

$$ 
\left\{ 
\begin{aligned}
\mathbf{X}_E &= \partial_E \varphi(\dot{\mathbf{E}}, \dot{\mathbf{S}}_v^{(i)}; \mathbf{E}, \mathbf{I}_1^*) = \sum_{i=1}^N \mu^{(i)}(\mathbf{E}, \mathbf{I}_1^*) (\dot{\mathbf{E}} - \mathbf{R}^{(i)} : \dot{\mathbf{S}}_v^{(i)}), \\
\mathbf{X}_{\mathbf{S}_v^{(i)}} &= \partial_{\mathbf{S}_v^{(i)}} \varphi(\dot{\mathbf{E}}, \dot{\mathbf{S}}_v^{(i)}; \mathbf{E}, \mathbf{I}_1^*) \\
&= -\mu^{(i)}(\mathbf{E}, \mathbf{I}_1^*) \mathbf{R}^{(i)} : (\dot{\mathbf{E}} - \mathbf{R}^{(i)} : \dot{\mathbf{S}}_v^{(i)}), \quad t = 1..N.
\end{aligned}
\right. 
\quad (16) 
$$

### 2.3.4. Dissipation analysis

The link between conjugate variables and thermodynamic forces is classically established by comparing two values of the dissipation $d$. The first is derived directly from the Clausius-Duhem inequality (2nd principle) and the second is obtained via the definition of the dissipation potential itself. The Clausius-Duhem inequality can be read as:

$$ d = \mathbf{S} : \dot{\mathbf{E}} - \dot{\psi} \geq 0 \quad (17) $$

where $\mathbf{S}$ is the second Piola-Kirchhoff tensor. On the other hand, the definition of the dissipation potential gives:

$$ d = \partial_E \varphi : \dot{\mathbf{E}} + \sum_{i=1}^N \partial_{\mathbf{S}_v^{(i)}} \varphi : \dot{\mathbf{S}}_v^{(i)} = \mathbf{X}_E : \dot{\mathbf{E}} + \sum_{i=1}^N \mathbf{X}_{\mathbf{S}_v^{(i)}} : \dot{\mathbf{S}}_v^{(i)}. \quad (18) $$

From (5) we can compute the time derivative of the free energy:

$$ \dot{\psi} = \mathbf{A}_E : \dot{\mathbf{E}} + A_{I_1^*} \dot{I}_1^* - \sum_{i=1}^N \dot{\mathbf{A}}_{\mathbf{S}_v^{(i)}} : \mathbf{S}_v^{(i)}, \quad (19) $$

leading to

$$ \dot{\psi} = (\mathbf{A}_E - \mathbf{S}_v) : \dot{\mathbf{E}} + \sum_{i=1}^N \mathbf{R}^{(i)} : \mathbf{S}_v^{(i)} : \dot{\mathbf{S}}_v^{(i)} + A_{I_1^*} \dot{I}_1^*. \quad (20) $$

During the active loadings $A_{I_1^*} = 0$ (from Eq. (8)) and during the neutral loadings and unloadings $\dot{I}_1^* = 0$ so that $A_{I_1^*} \dot{I}_1^* = 0$ for any loading paths.

Introducing expression (20) in (17) and using (14a) lead to express the dissipation as:

$$ d = (\mathbf{S} - 2\mathcal{F}(I_1, I_1^*) \partial_c W_e - p \mathbf{C}^{-1}) : \dot{\mathbf{E}} - \sum_{i=1}^N \mathbf{R}^{(i)} : \mathbf{S}_v^{(i)} : \dot{\mathbf{S}}_v^{(i)}. \quad (21) $$

On the other hand, one can use Eq. (18) with (16) to obtain another expression for the dissipation:

$$ d = \sum_{i=1}^{N} \mu^{(i)}(\mathbf{E}, \mathbf{I}_1^*) (\dot{\mathbf{E}} - \mathbf{R}^{(i)} : \dot{\mathbf{S}}_v^{(i)}) : \dot{\mathbf{E}} - \sum_{i=1}^{N} \mu^{(i)}(\mathbf{E}, \mathbf{I}_1^*) \mathbf{R}^{(i)} : (\dot{\mathbf{E}} - \mathbf{R}^{(i)} : \dot{\mathbf{S}}_v^{(i)}) : \dot{\mathbf{S}}_v^{(i)}. \quad (22) $$

Comparing then the Eqs. (21) and (22), we finally obtain:

$$ \mathbf{S}_v^{(i)} = \mu^{(i)}(\mathbf{E}, \mathbf{I}_1^*) : (\dot{\mathbf{E}} - \mathbf{R}^{(i)} : \dot{\mathbf{S}}_v^{(i)}), \quad t = 1..N \quad (23) $$

and

$$ \mathbf{S} = 2\mathcal{F} \partial_c W_e + \mathbf{S}_v + p \mathbf{C}^{-1}. \quad (24) $$

The stress tensor $\mathbf{S}$ can be therefore decomposed into an isochoric $\mathbf{S}_{iso}$ and volumetric parts $\mathbf{S}_{vol}$ such that:

$$ \mathbf{S}_{iso} = 2\mathcal{F} \partial_c W_e + \mathbf{S}_v, \quad \mathbf{S}_{vol} = p \mathbf{C}^{-1}, \quad (25) $$

with

$$ 
\left\{ 
\begin{array}{l}
\mathbf{R}^{(i)} : \dot{\mathbf{S}}_v^{(i)} + \frac{1}{\mu^{(i)}(\mathbf{E}, \mathbf{I}_1^*)} \mathbf{S}_v^{(i)} = \dot{\mathbf{E}} \\
\mathbf{S}_v^{(i)}|_{t=0} = \mathbf{0}
\end{array}
\right. 
\quad t = 1, \dots, N \quad (26) 
$$

0 being a second order tensor.

The stress is therefore fully known once the parameter set $\theta = (c_1, c_2, c_3, \{\mu_i\}_{i=1}^N, \eta, \alpha)$ is defined. For convenience and to allow direct comparison with experimental data, we introduce the first Piola-Kirchhoff (PK1) stress tensor, which is derived from $\mathbf{S}$ using the relation:

$$ \mathbf{P} = \mathbf{F}\mathbf{S}. $$

## 3. Uniaxial testing on sow skin tissues - Analytical model

In the present study, we focus on uniaxial tests performed on porcine skin tissues. These biological samples were selected because they exhibit key features relevant to biomechanics, such as nonlinear elasticity, viscoelasticity, and stress softening, along with a significant mechanical variability from one specimen to another. This inter-individual variability is representative of the heterogeneity encountered in biological tissues and provides a meaningful benchmark for validating constitutive models.

Although soft tissues are inherently anisotropic, we adopt here an effective isotropic formulation to describe the response along a single loading direction, consistent with the uniaxial nature of the experiments. This modeling choice allows us to investigate, in a controlled setting, the ability of the proposed formulation to reproduce key dissipative mechanisms, including rate-dependent stress relaxation and Mullins-type softening. The model is formulated using strain energy invariants (in particular $I_1$), and is thus naturally extendable to three-dimensional loading conditions. This uniaxial validation should therefore be seen as a first step towards the application of the model in more complex multiaxial biomechanical contexts.

### 3.1. Experimental procedure

Fresh sow perinea were collected from Ales slaughterhouse. The skin layers were dissected and stored in a sealed bags in a cryoconservant solution (10% glycerol 90% phosphate buffered saline) (Caro-Bretelle et al., 2016). The bags were placed in an ethanol bath and both were stored at -80 °C. 4 weeks later, the bags, still in the ethanol solution, were thawed at 5 °C for 6 h. Samples were dissected from the tissues, following the natural orientation of the vaginal tips (see Fig. 3 (left)). The samples were tested along their longest dimension, with an approximate 1:5 aspect ratio (5 mm width × 25 mm height), corresponding to the maximum usable length that could be clamped.

The uniaxial tensile tests were performed using a Zwick TH010 universal testing machine, controlled with the TestXpert II® software. The samples were clamped in jaws with sandpaper to avoid sliding, and digital image correlation was used to estimate displacement field. The samples were stretched until 30% following a cyclic loading-relaxation-unloading protocol (strain rate 0.01 s⁻¹). The loading, relaxation and unloading phases lasted 30, 120 and 30 s. The protocol is provided in Fig. 3 (right) and the experimental axial stress-time and elongation curves are shown for three samples obtained from different sows (Fig. 4).

### 3.2. Numerical implementation

Let us consider a cartesian frame $(\vec{x}, \vec{y}, \vec{z})$ and examine an idealized uniaxial test along the direction $\vec{x}$ for an incompressible, isotropic visco-hyperelastic material. In this case, the kinematics and stress expressions can be simplified:

$$ \mathbf{F} = \lambda \vec{x} \otimes \vec{x} + \lambda^{-1/2} \vec{y} \otimes \vec{y} + \lambda^{-1/2} \vec{z} \otimes \vec{z}, \quad \mathbf{C} = \lambda^2 \vec{x} \otimes \vec{x} + \lambda^{-1} \vec{y} \otimes \vec{y} + \lambda^{-1} \vec{z} \otimes \vec{z}, \quad (27) $$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Sampling location on the sow perineum (left). Mechanical testing protocol with imposed elongation, including loading, relaxation, and unloading phase (right).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Experimental axial stress and axial elongation over time during loading, relaxation, and unloading cycle for three samples from different sows.</figcaption></figure>

$$S = S_x \vec{x} \otimes \vec{x}, \qquad (28)$$

$$I_1 = \lambda^2 + 2\lambda^{-1}, \qquad (29)$$

where $\lambda$ is the axial stretch ratio.

The viscous stress tensor is evaluated from (26). Since both $\dot{E}$ is diagonal and the initial condition corresponds to a stress-free state, the solution of (26) preserves a diagonal structure. Therefore, $S_v^{(i)}$ takes the form:

$$S_v^{(i)} = S_{v_x}^{(i)} \vec{x} \otimes \vec{x} + S_{v_y}^{(i)} \vec{y} \otimes \vec{y} + S_{vz}^{(i)} \vec{z} \otimes \vec{z}, \quad i = 1..N. \qquad (30)$$

In this particular case, the stress components, are:

$$ \left\{ \begin{array}{l} S_x = 2\mathcal{F}\partial_{I_1} W_e + \sum_{i=1}^{N} S_{vx}^{(i)} + p\lambda^{-2}, \\ S_y = 2\mathcal{F}\partial_{I_1} W_e + \sum_{i=1}^{N} S_{vy}^{(i)} + p\lambda = 0, \\ S_z = 2\mathcal{F}\partial_{I_1} W_e + \sum_{i=1}^{N} S_{vz}^{(i)} + p\lambda = 0. \end{array} \right. \qquad (31) $$

For simplicity, the viscous parameter in (15) is assumed to depend only on the invariants of the deformation, yielding:

$$\mu^{(i)}(\mathbf{E}, I_1^*) = \mu^{(i)}(I_1, I_1^*). \qquad (32)$$

This choice reflects the isotropy of the material and the assumption that viscous effects are primarily governed by the deformation magnitude and loading history, rather than the full tensorial state. Such a formulation also facilitates calibration and reduces the number of internal variables, while still capturing essential history-dependent features like softening during unloading.

Eqs. (26) are reduced in this case to:

$$ \left\{ \begin{array}{l} 2\tau^{(i)} \dot{S}_{vx}^{(i)} + S_{vx}^{(i)} = 3G^{(i)} \tau^{(i)} \lambda \dot{\lambda}, \quad S_{vy}^{(i)} |_{t=0} = 0, \\ 2\tau^{(i)} \dot{S}_{vy}^{(i)} + S_{vy}^{(i)} = -3G^{(i)} \tau^{(i)} \lambda^{-2} \dot{\lambda}, \quad S_{vy}^{(i)} |_{t=0} = 0, \\ S_{vz}^{(i)} = S_{vy}^{(i)}, \quad S_{vz}^{(i)} |_{t=0} = 0, \end{array} \right. \qquad (33) $$

with $\tau^{(i)} = \frac{\mu^{(i)}}{3G^{(i)}}$. In these expressions we have introduced a new parameter, the relaxation time, $\tau^{(i)}$ (in seconds), which is commonly used in the literature on viscoelasticity. Moreover, we assume that $\mu^{(i)}$ and $\tau^{(i)}$ are both dependant on $I_1$ and $I_1^*$ via the formula:

$$\mu^{(i)} = \mathcal{F} \mu_0^{(i)}, \quad \tau^{(i)} = \mathcal{F} \tau_0^{(i)}, \quad i = 1...N, \qquad (34)$$

leading to

$$G^{(i)} = \frac{1}{3} \frac{\mu_0^{(i)}}{\tau_0^{(i)}} = \frac{1}{3} G_0^{(i)}, \quad i = 1...N, \qquad (35)$$

where $\mu_0^{(i)}$, $\tau_0^{(i)}$ and $G_0^{(i)}$ are respectively the initial viscosity coefficients, relaxation times and shear moduli in each viscoelastic branch.

This choice reflects the hypothesis that microstructural alterations responsible for the Mullins effect impact not only the elastic response but also the dissipative mechanisms. In biological soft tissues, such structural changes affect both stiffness and internal friction, enabling the modulation of viscoelastic parameters by the same softening function $\mathcal{F}$ (Xiang et al., 2019).

From (31b) and (31c) we can deduce that $S_{vy}^{(i)} = S_{vz}^{(i)}$ and hydrostatic pressure can be deduced from (31b) and replaced in (31a) leading to:

$$S_{xx} = 2(1 - \lambda^{-3}) \mathcal{F} \partial_{I_1} W_e + \sum_{i=1}^{N} S_{vx}^{(i)} - \lambda^{-3} \sum_{i=1}^{N} S_{vy}^{(i)}. \qquad (36)$$

In this part we focus on isotropic hyperelasticity depending on the invariants of the right Cauchy Green tensor (Gurtin, 1982):

$$W_e(C) = W_e(I_1, I_2, I_3). \qquad (37)$$

More, we consider a cubic form, that is frequently used in the literature to model the hyperelastic behavior of the skin Joodaki and Panzer (2018), Alliliche et al. (2023), proposed by Yeoh (Yeoh, 1993), depending only on $I_1$ and on material parameters $c_1, c_2$ and $c_3$, all of them positive:

$$W_e(C) = c_1(I_1 - 3) + c_2(I_1 - 3)^2 + c_3(I_1 - 3)^3. \qquad (38)$$
<figure><img src="image_7.png" /><figcaption>Fig. 5. Comparison between experimental and numerical data after fitting process: axial PK1 stress versus time and elongation for the sample 2.</figcaption></figure>

The derivation of $W_e$ can be simplified in:

$$ \partial_c W_e = \partial_{l_1} W_e. \qquad (39) $$

All the equations presented above were implemented using Python (version 3.11.12). A dedicated script was developed to solve the differential equations, compute the stress response under given loading conditions, and fit the model parameters. The identification was performed by minimizing a cost function defined as the sum of squared differences between experimental and simulated axial stress values:

$$ \sum_{k=1}^{T} (P_x^{\exp}(t_k) - P_x(t_k, \theta))^2, \qquad (40) $$

where $\theta$ is the set of model parameters and $T$ the number of time points.

An optimization procedure was performed independently for each of the three experimental tests to identify the model parameters. The optimization was carried out using Powell's method, a derivative-free algorithm well-suited for non-linear problems with potential local minima. The identified values for the hyperelastic, viscoelastic, and Mullins softening components are reported in Table 1. Despite individual variability, the optimized parameters remained consistent across the three samples, confirming the robustness of the model. To assess the quality of the fit, the root mean square error (RMSE) was computed between model and experimental data, both over time and as a function of elongation. The RMSE values in time remain low for all three samples, indicating a good agreement between the model and the experimental data over the loading-relaxation-unloading sequence. However, the RMSE values computed against elongation are higher. This discrepancy is mainly due to the end of the unloading phase and the beginning of the subsequent reloading, where the experimental axial stress becomes slightly negative. This behavior reflects irreversible mechanisms (e.g., plasticity or microstructural damage) that are not accounted for in the present model, which assumes fully recoverable hyper-viscoelastic behavior. The discrepancy also partially originates from the loading phase, where the use of a more advanced hyperelastic model could improve the fit. Nevertheless, such a refinement would require a larger number of parameters to be identified, increasing the model's complexity. Interestingly, the optimization yields $c_1 = 0$ in all three cases. This result can be explained by the pronounced stress decay observed during the relaxation phase. As a consequence, the influence of the hyperelastic term—particularly $c_1$, which primarily affects the initial stress increase—is effectively masked by the dominant viscous contribution.

For example, a comparison between experimental curves and model predictions is presented in Fig. 5, showing good agreement between experimental and numerical data.

### 3.3. Sensitivity analysis

Following the definition of the mechanical model and the identification of optimal parameters, a global sensitivity analysis was carried out to evaluate the robustness and identifiability of the parameters over

<table><caption>Table 1
Optimized parameters for the three experimental tests.</caption><thead><tr><th>Parameter</th><th>Sample 1</th><th>Sample 2</th><th>Sample 3</th></tr></thead><tbody><tr><td>c<sub>1</sub> (kPa)</td><td>0.0</td><td>10.6</td><td>0.0</td></tr><tr><td>c<sub>2</sub> (kPa)</td><td>50.1</td><td>63.4</td><td>50.3</td></tr><tr><td>c<sub>3</sub> (kPa)</td><td>412.6</td><td>907.0</td><td>409.0</td></tr><tr><td>G<sup>(1)</sup> (kPa)</td><td>404.0</td><td>431.1</td><td>424.6</td></tr><tr><td>G<sup>(2)</sup> (kPa)</td><td>153.5</td><td>183.0</td><td>164.9</td></tr><tr><td>τ<sup>(1)</sup><sub>0,1</sub> (s)</td><td>1.1</td><td>10.2</td><td>1.7</td></tr><tr><td>τ<sup>(2)</sup><sub>0,2</sub> (s)</td><td>12.2</td><td>23.9</td><td>18.1</td></tr><tr><td>α</td><td>0.12</td><td>0.14</td><td>0.32</td></tr><tr><td>η</td><td>0.69</td><td>0.65</td><td>0.74</td></tr><tr><td>RMSE (in time)</td><td>0.040</td><td>0.028</td><td>0.025</td></tr><tr><td>RMSE (in elongation)</td><td>0.159</td><td>0.193</td><td>0.151</td></tr></tbody></table>

their full plausible range. This analysis, based on Sobol indices, quantifies how the variability in model outputs can be attributed to each parameter and to their interactions. Specifically, the first-order indices ($S_1$) measure the isolated effect of each parameter, while the total-order indices ($S_T$) account for both individual and interaction effects. This variance-based approach provides a comprehensive understanding of parameter influence, independently of the optimized values, and helps inform future experimental design and model improvement. For further methodological details, see Sobol (2001).

The full parameter space included nine model parameters: three hyperelastic constants ($c_1, c_2, c_3$), two time constants ($\tau^{(1)}, \tau^{(2)}$), two viscous moduli ($G^{(1)}, G^{(2)}$), and the two Mullins parameters ($\alpha, \eta$). The bounds were defined as follows:

$$ \cdot c_i, \tau^{(i)}, G^{(i)} \in [0, 100] \text{ for } i = 1, 2, $$

$$ \cdot \alpha, \eta \in [0, 1[. $$

The sampling was performed using a Sobol sequence with 1024 points to ensure appropriate coverage of the parameter space. This size was chosen based on a convergence analysis, which confirmed that increasing the sample size beyond 1024 did not significantly change the computed Sobol indices. In addition, the bounds for each parameter were also defined based on a preliminary convergence analysis to ensure that the selected ranges adequately captured the sensitivity of the model. Each sampled parameter set was used to simulate the mechanical response under a single loading-relaxation-unloading-reloading cycle.

Five scalar indicators were extracted from these simulations:

1. Initial stress slope (start of loading) $S_0$ (kPa),

2. Initial relaxation slope (start of relaxation) $S_r$ (kPa/s),

3. Unloading slope (start of unloading) $S_u$ (kPa),

4. Relative relaxation amplitude $RR$ (kPa),

5. Hysteresis area $HA$ (kPa) for all positive axial stresses.

Fig. 6 summarizes the indicator definitions.
<figure><img src="image_8.png" /><figcaption>Fig. 6. Definition of key indicators used for Sobol analysis.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. First-order and total Sobol indices $S_1$ and $S_T$ for each indicator.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Second-order Sobol indices (parameter interactions) for each indicator.</figcaption></figure>

Fig. 7 present the first-order and total-order Sobol indices for all parameters across the five retained indicators. Interaction terms are shown in the second-order coupling matrices (Fig. 8).

The results confirm that each indicator is particularly sensitive to a distinct subset of parameters:

• Initial slope (loading): dominated by $c_1$ ($S_1 = S_T = 0.66$), followed by $G_0^{(1)}$ and $G_0^{(2)}$ (~0.15) with negligible interaction.

• Hysteresis area: again $c_1$ dominates ($S_1 = 0.65, S_T = 0.65$), but with minor contributions from $\tau_0^{(i)}, G_0^{(i)}, i = 1,2$. Second-order
interactions are limited, with a mild effect from $G_0^{(i)}/\tau_0^{(i)}$, $i = 1, 2$ coupling.

* Initial relaxation slope: $\tau_0^{(1)}$ and $\tau_0^{(2)}$ are most influential ($S_T \approx 0.5$, $S_1 \approx 0.35$), with $G_0^{(1)}$ and $G_0^{(2)}$ playing smaller roles. The $G_0^{(i)}/\tau_0^{(i)}$ $i = 1, 2$ interactions are the most relevant couplings.

* Relaxation amplitude: $c_1$ is the main contributor ($S_1 = 0.3$, $S_T = 0.48$), followed by $c_2$ and $G_0^{(1)}$, $G_0^{(2)}$. Couplings are moderate and mostly involve $c_1$ with $\tau_0^{(i)}$, $G_0^{(i)}$ and $c_2$ with $\tau_0^{(i)}$, $i = 1, 2$

* Unloading slope: mainly influenced by Mullins parameters—$\alpha$ ($S_1 = 0.55$, $S_T = 0.8$) and $\eta$ ($S_1 = 0.1$, $S_T = 0.35$)—alongside $c_1$ and $c_2$ at lower levels. Interaction is strongest between $\alpha$ and $\eta$.

As a result, each group of model parameters (elastic, viscous, and Mullins-related) is mainly identifiable from a specific portion of the test. The loading phase allows the identification of elastic and viscous moduli with limited coupling; the relaxation phase is essential to determine the time constants, though with moderate $G_0^{(i)}/\tau_0^{(i)}$, $i = 1, 2$ interactions; and the unloading-reloading sequence provides strong sensitivity to the Mullins parameters. The generally good agreement between first-order and total-order Sobol indices confirms that parameter influence remains largely independent. Altogether, the results support the relevance of the proposed protocol to robustly identify all parameters of the visco-hyperelastic-Mullins model from a single experimental test.

## 4. Conclusion

This study introduces a unified thermodynamical framework that combines Mullins effect and visco-hyperelasticity to model the behavior of filled elastomers and soft tissues under large deformations. By integrating both hyperelastic and viscous components within a finite strain thermodynamic approach, the proposed model provides an accurate and comprehensive description of the material's behavior during repeated loading cycles.

A key feature of the model is its ability to capture the coupling between strain rate and stress, which is essential for reproducing the observed mechanical response in soft biological materials. Furthermore, the model accounts for history-dependent stress softening (Mullins effect), which progressively appears during unloading and reloading, as confirmed by the phase-wise sensitivity analysis.

The optimization procedure resulted in a good agreement between experimental and simulated curves. To further evaluate the robustness of the identified parameter set, a global sensitivity analysis was conducted. It revealed a clear separation in the influence of parameters across different phases of the loading protocol. This indicates that each parameter family (elastic, viscous, Mullins) can be effectively constrained by distinct segments of the experiment, thereby supporting the identifiability and relevance of the proposed protocol.

As we continue to refine and validate this model, its potential applications extend across various domains, including biomechanics, prosthetics, and biomaterials. The ability to predict the mechanical behavior of soft tissues and filled elastomers under complex loading conditions can aid in the design of more effective implantable devices, prosthetic materials, and biomechanical simulations. Future work will focus on extending the model to multiaxial protocols (potentially requiring the introduction of additional invariants (Anssari-Benam et al., 2018)) and in vivo tissue conditions, where the interplay between microstructure, damage, and rate effects becomes critical. In this context, efforts will also focus on going beyond the current assumption of linear viscous behavior by introducing non-linear rate-dependent effects.

In addition, we aim to incorporate anisotropic behavior, commonly observed in biological tissues, not only in the hyperelastic response but also in the Mullins effect, by introducing direction-dependent constitutive features. The thermodynamic framework also opens the way for the integration of entropic effects, which are essential to accurately capture temperature-sensitive responses in soft materials.

## Nomenclature

<table><thead><tr><th>Symbol</th><th>Description</th></tr></thead><tbody><tr><td>F</td><td>Transformation gradient tensor</td></tr><tr><td>C</td><td>Cauchy-Green strain tensor</td></tr><tr><td>E</td><td>Green-Lagrange strain tensor</td></tr><tr><td>S</td><td>Second Piola-Kirchhoff tensor</td></tr><tr><td>S<sub>iso</sub></td><td>Second isovolumic Piola-Kirchhoff tensor</td></tr><tr><td>S<sub>v</sub></td><td>Second Piola-Kirchhoff viscous tensor</td></tr><tr><td>P</td><td>First Piola-Kirchhoff tensor</td></tr><tr><td>1</td><td>Second order identity tensor</td></tr><tr><td>0</td><td>Second order zero tensor</td></tr><tr><td>p</td><td>Hydrostatic pressure</td></tr><tr><td>I<sub>1</sub></td><td>First invariant of C</td></tr><tr><td>λ<sub>i</sub>, i = 1..3</td><td>Stretch ratios</td></tr><tr><td>I<sub>1</sub><sup>*</sup></td><td>Maximum historical value of I<sub>1</sub></td></tr><tr><td>W<sub>e</sub></td><td>Hyperelastic strain energy density</td></tr><tr><td>F</td><td>Softening function (Mullins effect)</td></tr><tr><td>ψ</td><td>Helmholtz free energy</td></tr><tr><td>G</td><td>Thermodynamic potential</td></tr><tr><td>J</td><td>Jacobian determinant of F</td></tr><tr><td>μ<sup>(i)</sup></td><td>Viscosity coefficient of branch i</td></tr><tr><td>G<sup>(i)</sup></td><td>Shear modulus of branch i</td></tr><tr><td>τ<sup>(i)</sup></td><td>Relaxation time of branch i</td></tr><tr><td>R<sup>(i)</sup></td><td>Fourth-order compliance-like tensor of branch i</td></tr><tr><td>∂<sub>var</sub></td><td>Partial derivative with respect to the variable var</td></tr></tbody></table>

## CRediT authorship contribution statement

**Anne-Sophie Caro:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **André Chrysochoos:** Writing – review & editing, Writing – original draft, Methodology, Formal analysis, Conceptualization. **Sarah Iaquinta:** Writing – review & editing, Writing – original draft, Visualization, Investigation, Data curation. **Grégory Chagnon:** Writing – review & editing, Formal analysis, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Supplementary material

The Python codes containing the routines to reproduce the results presented in this article, as well as the experimental data, are provided in the Zenodo repository https://doi.org/10.5281/zenodo.15740736.

## Data availability

We add in supplementary material the link to a Zenodo repository for data and code.

## References

Alliliche, W., Renaud, C., Cros, J.-M., Feng, Z.-Q., 2023. Numerical simulation of mechanical tests on a living skin using anisotropic hyperelastic law. J. Mech. Behav. Biomed. Mater. 141, 105755.

Anssari-Benam, A., 2023. Continuous softening up to the onset of failure: A hyperelastic modelling approach with intrinsic softening for isotropic incompressible soft solids. Mech. Res. Commun. 132, 104183. http://dx.doi.org/10.1016/j.mechrescom.2023.104183.
Anssari-Benam, A., Hossain, M., 2023. A pseudo-hyperelastic model incorporating the rate effects fordiani isotropic rubber-like materials. J. Mech. Phys. Solids 179, 105347. http://dx.doi.org/10.1016/j.jmps.2023.105347.

Anssari-Benam, A., Hossain, M., 2024. A unified pseudo-elastic model of continuous and discontinuous softening in the finite deformation of isotropic soft solids. Int. J. Solids Struct. 290, 112670. http://dx.doi.org/10.1016/j.jmbbm.2024.106502.

Anssari-Benam, A., Saccomandi, G., 2024. Modelling the rate-dependent mechanical behaviour of the brain tissue. J. Mech. Behav. Biomed. Mater. 153, 106502. http://dx.doi.org/10.1016/j.jmbbm.2024.106502.

Anssari-Benam, A., Tseng, Y.-T., Bucchi, A., 2018. A transverse isotropic constitutive model for the aortic valve tissue incorporating rate-dependency and fibre dispersion: Application to biaxial deformation. J. Mech. Behav. Biomed. Mater. 85, 80-93.

Cantournet, S., Desmorat, R., J., B., 2009. Mullins effect and cyclic stress softening of filled elastomers by internal sliding and friction thermodynamics model. Int. J. Solids Struct. (46), 2255-2264.

Caro-Bretelle, A., Ienny, P., Leger, R., Corn, S., Bazin, I., Bretelle, F., 2016. Constitutive modeling of stress softening and permanent set in a porcine skin tissue: Impact of the storage preservation. J. Biomech. 49 (13), 2863-2869.

Chagnon, E., Gornet, L., Marckmann, G., Charrier, P., 2004. On the relevance of continuum damage mechanics as applied to the mullins effect: theory, experiments and numerical implementation. J. Mech. Phys. Solids 52, 1627-1650.

Diani, J., Fayolle, B., Gilormini, P., 2009. A review on the mullins effect. Eur. Polym. J. 45 (3), 601-612.

Gurtin, M.E., 1982. An Introduction to Continuum Mechanics. In: Mathematics in Science and Engineering, Academic Press, New York.

Halphen, B., Nguyen, Q.S., 1975. Sur les matériaux standard généralisés. Comput. Methods Appl. Mech. Engrg. 14 (1), 39-63.

Holzapfel, G.A., Gasser, T.C., 2001. A viscoelastic model for fiber-reinforced composites at finite strains: Continuum basis, computational aspects and applications. Comput. Methods Appl. Mech. Engrg. 190, 4379-4403.

Holzapfel, G.A., Ogden, R.W., 2020. A damage model for collagen fibres with an application to collagenous soft tissues. Proc. R. Soc. A: Math. Phys. Eng. Sci. 476 (2236), 20190821.

Joodaki, H., Panzer, M.B., 2018. Skin mechanical properties and modeling: A review. Proc. Inst. Mech. Eng. H 232 (4), 323-343.

Marckmann, G., Verron, E., 2006. Comparison of hyperelastic models for rubber-like materials. Rubber Chem. Technol. 79 (5), 835-858.

Mullins, L., 1969. Softening of rubber by deformation. Rubber Chem. Technol. 1 (42), 339-362.

Mullins, L., Tobin, N.R., 1965. Stress softening in rubber vulcanizates. J. Appl. Polym. Sci. 9, 2993-3009.

Ortún-Terrazas, J., Cegoñino, J., Santana-Penín, U., Santana-Mora, U., Pérez del Palomar, A., 2019. A porous fibrous hyperelastic damage model for human periodontal ligament: Application of a microcomputerized tomography finite element model. Int. J. Numer. Meth. Biomed. Eng. 35 (4), e3176.

Peña, E., Martins, P., Macarenhas, T., Natal Jorge, R.M., Ferreira, A., Doblaré, M., Calvo, B., 2011. Mechanical characterization of the softening behavior of human vaginal tissue. J. Mech. Behav. Biomed. Mater. 4 (3), 275-283.

Rebouah, M., Chagnon, G., 2014. Extension of classical viscoelastic models in large deformation to anisotropy and stress softening. Int. J. Non-Linear Mech. 61, 54-64.

Rickaby, S.R., Scott, N.H., 2014. Multicyclic modelling of softening in biological tissue. IMA J. Appl. Math. 79 (6), 1107-1125, Cited by: 6; All Open Access, Green Open Access.

Simo, J., 1987. On a fully three-dimensional finite-strain viscoelastic damage model: formulation and computational aspects. Comput. Methods Appl. Mech. Engrg. 60, 153-173.

Simo, J.C., Hughes, T.J.R., 1998. Computational Inelasticity. In: Interdisciplinary Applied Mathematics, vol. 7, Springer, New York.

Sobol, I., 2001. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math. Comput. Simulation 55 (1), 271-280.

Wenguang, L., 2016. Damage models for soft tissues: A survey. J. Med. Biol. Engng. 36 (3), 2199-4757.

Xiang, Y., Zhong, D., Wang, P., Yin, T., Zhou, H., Yu, H., Baliga, C., Qu, S., Yang, W., 2019. A physically based visco-hyperelastic constitutive model for soft materials. J. Mech. Phys. Solids 128, 208-218.

Yeoh, O., 1993. Some forms of the strain energy function for rubber. Rubber Chem. Technol. 66 (5), 754-771.