A Review on Damage and Rupture Modelling for Soft Tissues

Abstract

Computational modelling of damage and rupture of non-connective and connective soft tissues due to pathological and supra-physiological mechanisms is vital in the fundamental understanding of failures. Recent advancements in soft tissue damage models play an essential role in developing artificial tissues, medical devices/implants, and surgical intervention practices. The current article reviews the recently developed damage models and rupture models that considered the microstructure of the tissues. Earlier review works presented damage and rupture separately, wherein this work reviews both damage and rupture in soft tissues. Wherein the present article provides a detailed review of various models on the damage evolution and tear in soft tissues focusing on key conceptual ideas, advantages, limitations, and challenges. Some key challenges of damage and rupture models are outlined in the article, which helps extend the present damage and rupture models to various soft tissues.

1. Introduction

Soft tissue refers to non-mineralised tissue or organ in living creatures that are fibrous, which connects, support, or surrounds other structures and organs of the body [1,2,3]. Soft tissues are classified into two major groups: connective tissues and non-connective tissues. Tendons, skin, fat, ligaments fall under connective tissues, while muscles, nerves, and blood vessels fall under non-connective tissues [2,3,4]. Each tissue undergoes a specific mechanical loading according to its functionality. For example, blood vessels and arterial valves experience cyclic circumferential loading due to blood circulation; the eye tissue experiences constant tension due to intraocular pressure; the articular cartilages covering the ends of bones experience constant compression and friction.

In medical sciences, the study of damage in soft tissues is essential in both physiological and pathological conditions. Damage in soft tissues is intrinsic either due to excessive external mechanical loading [5,6] or pathological condition [7,8]. A pathological condition such as glaucoma in the human eye leads to excessive accumulation of aqueous humour that causes high intraocular pressure and damages the optic nerve head, leading to vision loss (Figure 1b) [9,10]. While damage due to external mechanical loading occurs when the applied load exceeds the physiological limit of the tissue (supra-physiological load). For example, ligament tear, which athletes encounter, occurs due to excessive or repetitive tensile loading during sport (Figure 1a) [5]. During extensive physical activities, at times, ligaments in the knee undergo supra-physiological loading causing flexed knee and lateral rotation of the tibia, resulting in anterior cruciate ligament tear [5,11]. Knowledge of the damage in soft tissues would help to come up with improved clinical practices.

Understanding the biomechanics of the soft tissues is essential for modelling damage. Experimental tests provide fundamental insights into underlying deformation and damage mechanisms. Since the inception of biomechanics, tissues from animals, cadavers, and volunteer living subjects have been used extensively for experimental tests [12,13,14,15,16,17,18]. Despite providing valuable information, the experimental approach suffers several limitations. These include stringent ethical and animal certifications, limited availability of subjects and non-standard experimental procedures. In addition, the tissue responses captured in the experiments are passive, i.e., the effect of muscle cells is neglected. Moreover, the risk of tissue degradation and tissue to tissue variations are always present in the experiments [19]. Therefore, the modelling and simulation of tissue damage is a cost-effective alternative in conceptualizing the improved clinical practices, tissue engineering, and medical device development with limited experimental validation. For instance, modelling the inelastic behaviour of arteries during angioplasty would help in optimising the surgical procedure. It is possible to understand the rupture phenomena of aortic aneurysms under in-vivo conditions using rupture simulations, which would help in its diagnosis. Furthermore, rupture simulations that simulate the in-vivo loading conditions would help in designing the artificial implants.

The inelastic deformation behaviour of soft tissue depends on the stress state and extent of the damage. The stress state is determined using the material constitutive relation: it governs the displacement response of the material to the mechanical loading. While the damage is characterised using initiation and evolution laws. And the failure is captured through rupture laws. Hence, the present study discusses various soft tissue constitutive, damage, and rupture models.

Soft tissue consists of cells, elastin, collagen, and a non-mineralized ground matrix. Collagen fibres are of high stiffness compared to the rest of the constituents of the tissue, and they majorly contribute towards the overall stiffness of the tissue [20,21]. Based on their constituents and structure, soft tissues exhibit non-linear, anisotropic, hyperelastic, and viscoelastic behaviour [21]. Constitutive behaviour varies from tissue to tissue, depending on the collagen fibre distribution and orientation inside the tissue. For example, in connective tissues such as ligaments and tendons, the collagen fibre orientation is regular and unidirectional [22,23]. While in non-connective tissues, such as arterial walls, collagen fibres orientation is irregular and multi-directional [6,20,21]. In the arterial wall, collagen fibre distribution is arranged in a double-helical pattern [24]. Hence due to the tissue-specific arrangement of the collagen and other constituents, ligaments possess high tensile strength, i.e., 50–100 MPa, while the arteries possess low tensile strength, i.e., 0.3–0.8 MPa [2,23]. However, both the tissues exhibit anisotropic hyperelastic behaviour.

Figure 2 represents a typical stress-strain curve of skin subjected to uniaxial tension. The skin constitutes three layers: epidermis, dermis, and hypodermis [25]. Wherein the epidermis is the thin outermost cellular layer. The dermis is the middle layer with elastin and collagen fibres embedded in the extrafibrillar matrix, providing strength. The hypodermis is the innermost supportive layer that primarily consists of adipose tissue [25,26]. The skin majorly consists of collagenous fibres that account for 60 to 80% of dry weight, and the fibres are woven into a rhombic shaped pattern [2,25]. In general, most collagenous soft tissues manifest a typical J-shaped stress-strain behaviour, similar to that of the skin, as illustrated in Figure 2 [2,25].

In the current review, the mechanical behaviour of the skin is explained with a macroscopic perspective. However, the nanomechanics of collagen fibres where hydrogen and covalent bonds in the protein backbone plays a key role also contributes to the macroscopic mechanical behaviour [26,27]. The deformation behaviour of the skin can be divided into three phases:

• The initial region of the stress-strain response, i.e., phase-I (toe region). The soft tissue’s mechanical behaviour in this region is similar to a soft isotropic rubber sheet. The collagen fibres are in a relaxed state, and they appear wavy and crimped. Therefore, very low stress is required for attaining large deformation without stretching the collagen fibres. As a result, the mechanical behaviour in phase-I is approximately linear, and the elastic modulus is low (0.1–2 MPa) [25,28].
• In phase-II (heel region), the tissue exhibits a highly non-linear mechanical behaviour [25]. The collagen fibres get uncrimped as they elongate with the increase in the load. The elongated fibres slide into the matrix and align themselves to the direction of load, thereby increasing the load-carrying capacity.
• In Phase-III (linear region), the tissue exhibits stiffer and linear behaviour. Most of the fibres get aligned to the loading direction; hence, no crimp pattern is observed. The aligned and straightened fibres resist the load, making the tissue stiffer and linear in mechanical behaviour [23,28]. Beyond phase III, ultimate tensile strength is reached, resulting in tissue rupture.

Damage in the soft tissues is defined as “injury or harm that reduces value or usefulness” [29], and rupture results from accumulated damage, which is catastrophic. Damage models phenomenologically capture the damage in the soft tissue due to pathological or supra-physiological conditions. On the other hand, fracture mechanics concepts alone cannot quantify the soft tissue rupture that exhibits a toughening mechanism leading to high defect tolerance [30,31]. Hence rupture in soft tissues needs a damage model along with fracture mechanics. Earlier review works [6,7,32,33] presented damage and rupture separately, wherein this paper reviews both damage and rupture in soft tissues. The earlier reported reviews are either tissue-specific or confined to phenomenological models. The current article reviews the recently developed damage models and rupture models that considered the microstructure of the tissues. In Section 2, continuum kinematics and soft tissue material constitutive model is presented that helps in reviewing the damage and rupture models in a unified manner. The damage models for soft tissues are reviewed in Section 3, which are classified into three groups, namely (1) continuum damage mechanics (CDM), (2) pseudo-elasticity, and (3) softening hyperelasticity. In Section 4, soft tissue rupture models are reviewed and classified into three groups: (1) extended finite element method (XFEM), (2) cohesive zone modelling (CZM), and (3) crack phase-field method (CPFM). A summary highlighting all the damage and rupture models and their respective challenges is discussed in Section 5, and the final concluding remarks are given in Section 6.

2. Kinematics and Constitutive Model

This section is the preliminary to provide a unified representation for all the damage and rupture models. The essential kinematics of the continuum mechanics needed for defining the damage, rupture and material models are presented.

Figure 3 shows a body that occupies the domain ${\mathrm{B}}_0$ in reference configuration and occupies domain $\mathrm{B}$ in the current configuration. Let $\mathit{X}$ denote the material point in the reference configuration, and $\mathit{x}$ denote the same material point in the current configuration. Deformation map $\mathit{F}=\frac{\partial x}{\partial X}$ maps reference configuration to final configuration and Jacobian $J=\det \left(\mathit{F}\right)>0$. Right Cauchy–Green tensor is $\mathit{C}={\mathit{F}}^T\mathit{F}$. Modified right Cauchy–Green tensor $\overline{\mathit{C}}=J^{-\frac{2}{3}}\mathit{C}$ is used for pure distortion. ${\mathit{a}}_1$ and ${\mathit{a}}_2$ be the orientation of two fibre families in the reference configuration. ${\mathit{a}}_1^{\prime }$ and ${\mathit{a}}_2^{\prime }$ be their corresponding orientations in the current configuration. They are related as ${\mathit{a}}_1^{\prime }=\mathit{F}{\mathit{a}}_1$ and ${\mathit{a}}_2^{\prime }=\mathit{F}{\mathit{a}}_2$.

Invariants associated with Cauchy–Green deformation tensor are

$$I_1=tr\left(\mathit{C}\right), I_2=tr\left(cof\left(\mathit{C}\right)\right)\mathrm{and} I_3=det\left(\mathit{C}\right)$$

(1)

Invariants related to fibres are given as

$$I_4=\mathit{C}:\left({\mathit{a}}_1\otimes {\mathit{a}}_1\right)={\mathit{a}}_1^{\prime }.{\mathit{a}}_1^{\prime }, I_6=\mathit{C}:\left({\mathit{a}}_2\otimes {\mathit{a}}_2\right)={\mathit{a}}_2^{\prime }.{\mathit{a}}_2^{\prime }$$

(2)

The majority of damage and rupture models reviewed have used the HGO (Holzapfel Gasser and Ogden) model for defining the material constitutive response [34]. Hence, the HGO model is briefly discussed here. In the HGO model, the isotropic and anisotropic parts of the strain energy function is decomposed as:

$$\mathsf{\Psi }={\mathsf{\Psi }}_{vol}+{\mathsf{\Psi }}_{dev}+{\mathsf{\Psi }}_{ani}$$

(3)

where

$${\mathsf{\Psi }}_{vol}=\frac{1}{2}K{\left(J-1\right)}^2$$

(4)

$${\mathsf{\Psi }}_{dev}=\frac{\mu }{2}\left({\overline{I}}_1-3\right)$$

(5)

$${\mathsf{\Psi }}_{ani}={\displaystyle \sum }_{i=4,6}\frac{k_1}{2k_2}\left[\exp \left(k_2\langle {\overline{E}}_i^{ }\rangle {}^2\right)-1\right]$$

(6)

where the first two components are dilatational $\left({\mathsf{\Psi }}_{vol}\right)$ and deviatoric $\left({\mathsf{\Psi }}_{dev}\right)$ parts of the isotropic response and last component correspond contribution of the two families of fibres $\left({\mathsf{\Psi }}_{ani}\right)$. The invariants are defined for the modified right Cauchy–Green tensor. $K$ is the bulk modulus, $\mu$ is the shear modulus, $k_1$ and $k_2$ are the parameters that define the contribution of fibres. The Macaulay bracket is defined as ($\langle ■\rangle =\left(\left|■\right|+■\right)\frac{1}{2}$) that assumes the fibres can only support tension. ${\overline{E}}_i^{ }$ is a strain-like quantity that represents the deformation of the family of fibres ${\mathit{a}}_l \left(l=1, 2\right)$ and it is defined as

$${\overline{E}}_i^{ }\stackrel{\mathrm{def}}{=}\kappa \left({\overline{I}}_1-3\right)+\left(1-3\kappa \right)\left({\overline{I}}_i-1\right), i=4, 6$$

(7)

This model assumes that the collagen fibres orientation and distribution for the respective family of fibres, which are described along a preferred mean direction. Wherein, $\kappa$ is the radial dispersion parameter that defines the radial symmetry dispersion of the fibre orientation. $\kappa$ is defined based on the fibre orientation density function $\rho$ and the fibre orientation, $\mathsf{\Theta }$, as

$$\kappa =\frac{1}{4}{{\displaystyle \int }}_0^{\pi }\rho \left(\mathsf{\Theta }\right){\sin }^3\mathsf{\Theta }d\mathsf{\Theta }$$

(8)

The radial dispersion parameter $\kappa \in \left[0,\frac{1}{3}\right]$, where $\kappa =0$ describes perfectly aligned fibres without dispersion and $\kappa =\frac{1}{3}$ describes random distribution where the material behaves as isotropic. The strain, as with quantity ${\overline{E}}_i^{ }$, becomes

$${\overline{E}}_i^{ }=\left({\overline{I}}_i-1\right)\mathrm{for}\kappa =0$$

(9)

$${\overline{E}}_i^{ }=\frac{\left({\overline{I}}_1-1\right)}{3}\mathrm{for}\kappa =\frac{1}{3}$$

(10)

Further, the anisotropic strain energy function (6) for each family of fibres can be additively decomposed as ${\mathsf{\Psi }}_{ani}={\mathsf{\Psi }}_{f1}+{\mathsf{\Psi }}_{f2}$ and for aligned fibres the strain energy is given as

$${\mathsf{\Psi }}_{af1}=\frac{k_1}{2k_2}\left[\exp k_2{\left({\overline{I}}_4-1\right)}^2-1\right]$$

(11)

and

$${\mathsf{\Psi }}_{af2}=\frac{k_1}{2k_2}\left[\exp k_2{\left({\overline{I}}_6-1\right)}^2-1\right]$$

(12)

Similarly, for distributed fibres, the strain energy density (${\mathsf{\Psi }}_{ani}={\mathsf{\Psi }}_{df1}+{\mathsf{\Psi }}_{df2}$) is given as

$${\mathsf{\Psi }}_{df1}=\frac{k_1}{2k_2}\left[\exp k_2\langle \kappa \left({\overline{I}}_1-3\right)+\left(1-3\kappa \right){\left({\overline{I}}_4-1\right)\rangle }^2-1\right]$$

(13)

and

$${\mathsf{\Psi }}_{df2}=\frac{k_1}{2k_2}\left[\exp k_2\langle \kappa \left({\overline{I}}_1-3\right)+\left(1-3\kappa \right){\left({\overline{I}}_6-1\right)\rangle }^2-1\right]$$

(14)

3. Damage Models

Damage in the soft tissues occurs when the applied load goes beyond the physiological range, as shown in Figure 2. The physiological limit of the soft tissue lies in phase-II [28]; when the applied load exceeds this limit, the damage at the microscopic level initiates, resulting in microscopic failure. Damage in a material is a progressive physical phenomenon that leads to failure, as shown in Figure 2. Experimental studies have reported four phenomena associated with the damage in soft tissues: (a) Mullins effect, (b) hysteresis, (c) permanent set, and (d) rupture [6,35]. Under the cyclic loading of soft tissue, the stress required after the first cycle for the same stretch reduces for the consecutive cycles (Figure 4a). This stress softening phenomenon in materials under cyclic loading is known as the Mullins effect [36]. The maximum load governs the Mullins effect during the loading cycles. In the literature, the Mullins effect is extensively studied by several authors on arteries [37,38], vaginal tissue [39], etc. The continuous softening of the material subjected to cyclic loading under constant load is known as hysteresis (Figure 4b). The softening phenomenon continues until a saturation point is reached [6,35,40,41,42]. Hysteresis is used to precondition the soft tissues at lower loads to overcome the shape effects in experimental tests [43,44]. Hysteresis studies are often reported for arteries by Pena et al. [39] and Balzani et al. [37,42]. The inelastic behaviour that occurs due to accumulated strain in the soft tissue due to load is known as the permanent set [35]. This phenomenon is also studied extensively for soft tissues as with arteries [45,46] and bioprosthetic heart valves [47]. Lastly, the accumulated microscopic damage in the tissue leads to a macroscopic failure known as rupture. The tissue rupture may arise due to the failure of the matrix or rupture in fibres. In the literature, rupture studies are reported in arteries [48,49], corneas [50,51], skin [52], ligaments [53], etc.

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, and the experimental conclusions that can be drawn.

The aforementioned phenomena require a damage model along with a constitutive response to capture the deformation behaviour of the tissue. In general, the damage models for soft tissues available in the literature can be broadly classified into three categories: (a) models based on continuum damage mechanics (CDM), (b) theory of pseudo-elasticity, and (c) the softening hyperelasticity approach [6]. The CDM gives a continuum level description for damage phenomena, i.e., it provides damage initiation, propagation, and microscopic failure [54,55]. It is based on an irreversible thermodynamics process where the Clausius–Duhem inequality is used for defining the internal state variables and internal dissipation. In CDM, the damage is modelled using the state variables that defines the onset of the damage and govern the degraded material response. For soft tissues, CDM has been used to analyse the Mullins effect, permanent set, and tissue rupture [37,56,57]. In CDM, the damage initiation is characterised based on the concept of equivalent strain for undamaged material introduced by Simo and Ju [54]. Later, Miehe used the same concept to define discontinuous damage characterised by the maximum strain in the loading path and continuous damage characterised by strain-rate-dependent [7]. Comparatively, CDM requires a higher number of material parameters to define the damage.

Phenomenologically, pseudo-elasticity models the damage in the soft tissues using two elastic material models, one during loading and another while unloading. In this approach, the stress-strain behaviour in cyclic loading is defined using fewer material parameters than CDM. The pseudo-elasticity approach has been used to model arterial walls [58]; Mullins effect in rubber [59,60]; Mullins effect and permanent set in brain tissue and arteries [17,38,61]. The softening hyperelasticity technique was introduced as a substitute to CDM and pseudo-elasticity by Volokh [62], wherein the constitutive response is incorporated with strain softening using the material constants named as energy limiters. In this approach, the internal variables to quantify damage, initiation condition, and the evolution equation for damage are not required. Therefore, it is more straightforward than CDM and pseudo-elasticity. This approach has been adopted by Li and Lou [63] in modelling human and animal skin. The present section systematically reviews the modelling techniques such as CDM, pseudo-elasticity, and softening hyperelasticity. A comparison of the three damage modelling techniques is summarised in

Table 1. Summary of damage modelling techniques.

Modelling Considerations	Continuum Damage Mechanics	Pseudo-Elasticity	Softening Hyperelasticity
Strain energy density	$\mathsf{\Psi }={\mathsf{\Psi }}_{vol}+\left(1-d\right)\left({\mathsf{\Psi }}_{dev}+{\mathsf{\Psi }}_{ani}\right)$	$\mathsf{\Psi }={\mathsf{\Psi }}_{vol}+{\mathsf{\Psi }}_{dev}+\eta {\mathsf{\Psi }}_{ani}+\varphi \left(\eta \right)$	$\mathsf{\Psi }=\varphi \left[1-\exp \left(-\frac{W}{\varphi }\right)\right]$
Damage parameter	$\mathrm{Kachanov}\mathrm{reduction}\mathrm{factor}\left(1-d\right)$	$\eta$–damage variable $\varphi$–damage function	$\varphi$–energy limiter
Damage initiation	$$\begin{gathered} \mathrm{Equivalent}\mathrm{strain}\mathrm{based}\mathrm{on}\mathrm{Simo} \\ \sqrt{2{\mathsf{\Psi }}_{\alpha }}-r_{\alpha }\le 0 \end{gathered}$$ $r_{\alpha }$ –damage threshold at the current time.	Based on the critical stretch in fibres, i.e., $I_i\ge I_{ic}, i=4.6$	Strain softening incorporated using energy limiters.
Damage evolution	Based on the model. Discontinuous damage modelled with the maximum strain in the loading path. Continuous damage is strain-rate-dependent	$\eta =1-\frac{1}{r}\mathrm{erf}\left[\frac{{\mathsf{\Psi }}_{ani}^{max}-{\mathsf{\Psi }}_{ani}}{m}\right]$ $\varphi \left(\eta \right)=m\eta \log \eta +\left(1-\eta \right)\left[m+{\mathsf{\Psi }}_{ani}\left(I_{ic}\right)\right]$
Thermodynamic consideration	$\mathrm{Clausius}\text{–}\mathrm{Duham}\mathrm{inequality}\mathrm{used}\mathrm{to}\mathrm{define}\mathrm{damage}\mathrm{threshold}(r_{\alpha }$), which is maximum strain energy without damage.	$\mathrm{In}\mathrm{the}\mathrm{primary}\mathrm{loading}\mathrm{curve}\eta =1$$\mathrm{and}\varphi \left(1\right)=0$$\mathrm{and}\mathrm{for}\mathrm{the}\mathrm{subsequent}\mathrm{unloading}\mathrm{and}\mathrm{reloading},\eta <1$$\mathrm{and}\varphi \left(\eta \right)$ evolves, which is consistent with Clausius–Duham inequality.	Energy limiters activate the irreversible damage and dissipation that ensure the thermodynamic stability of the model.

.

3.1. Continuum Damage Mechanics (CDM)

In this section, various damage models based on continuum damage mechanics are reviewed chronologically, and all the reviewed CDM based damage models are summarized in

Table 2. Summary of various damage models based on continuum damage mechanics (CDM).

References	Tissue	Tissue Structure	Damage	Modelling Features
Balanco et al. [64]	Soft tissue with fibres	Isotropic matrix and collagen fibres	Matrix and fibres	(1) Anisotropic, incompressible, (2) HGO strain energy function, (3) three parameters to define the damage, (4) continuum damage based on Simo and Ju [55].
Comellas et al. [65]	Rectus sheath	Isotropic matrix	Matrix	(1) Isotropic, incompressible, (2) strain energy function: neo-Hookean and Ogden, (3) one parameter related to the softening effect.
Polindara et al. [66]	Blood vessel	Isotropic matrix and collagen fibres	Fibres	(1) Anisotropic, incompressible, (2) HGO strain energy function, (3) two parameters to define the damage, (4) continuum damage based on Simo and Ju [55].
Ferreira et al. [67]	Arteries	Isotropic matrix and collagen fibres	Matrix and fibres	(1) Anisotropic, incompressible, (2) HGO strain energy function, (3) seven parameters to define the damage, (4) continuum damage based on Simo [68].
Rausch et al. [68]	Soft tissue with fibres	Isotropic matrix and collagen fibres	Matrix and fibres	(1) Anisotropic, incompressible, (2) HGO strain energy function, (3) two parameters to define the damage, (4) continuum damage based on Simo [68].
Fathi et al. [69]	Soft tissue with fibres	Isotropic matrix and collagen fibres	Matrix and fibres	(1) Anisotropic, incompressible, (2) HGO strain energy function, (3) six parameters to define the damage, (4) continuum damage based on Simo and Ju [55].
Gao et al. [70]	Annulus fibrosus	Isotropic matrix and collagen fibres	Matrix and fibres	(1) Anisotropic, incompressible, (2) HGO strain energy function, (3) four parameters to define the damage, (4) damage model is based on [38].
Mousavi et al. [71]	Ascending thoracic aortic aneurysms	Smooth muscle cells and collagen fibres distributed in the elastin matrix	Elastin matrix and collagen fibres	(1) Anisotropic, incompressible, (2) HGO strain energy function for fibres, (3) three parameters to define the damage, (4) damage model is based on linear softening by [66]
Ghasemi et al. [35]	Arteries	Elastin and collagen fibres in an isotropic matrix	Elastin fibres and collagen fibres	(1) Anisotropic, incompressible, (2) HGO strain energy function, (3) four parameters to define the damage, (4) continuum damage based on [58].

.

Blanco et al. [64] have developed a continuum damage model for soft tissues by including damage in fibre and matrix. The model was aimed to define the relationship between the mesoscopic structural mechanisms and macroscopic material parameters that dominate the inelastic behaviour in the soft fibrous tissues. The strain energy function, including the damage, is defined with the help of Neo-Hookean for the ground matrix and the Holzapfel Gasser Ogden (HGO) model for collagen fibres [24]. The strain energy function is defined as:

$$\mathsf{\Psi }=\frac{1}{2}{\mathsf{\Psi }}_{vol}+\left(1-d_{dev}\right){\mathsf{\Psi }}_{dev}+\left(1-d_{f1}\right){\mathsf{\Psi }}_{af1}+\left(1-d_{f2}\right){\mathsf{\Psi }}_{af2}$$

(15)

And $\left\{\left(1-d_{dev}\right),\left(1-d_{af1}\right),\left(1-d_{af2}\right)\right\}$ are the reduction factors that introduce the inelastic material damage that occurs in the matrix and every set of fibres, complying with $0\le d_{\alpha }\le 1$ for $\alpha =\left\{dev,af1,af2\right\}$. The damage parameters are expressed as:

$$d_{\alpha }=1-\frac{q_{\alpha }\left(r_{\alpha }\right)}{r_{\alpha }},\alpha =\left\{dev,af1,af2\right\}$$

(16)

Here $q_{\alpha }$ is defined as internal stress such as a variable whose variation defines the softening effect and damage threshold in the current configuration is defined by stress such as a variable, $r_{\alpha }$ satisfying $r_{\alpha }\le r_{\alpha }^0$, with $r_{\alpha }^0$ the initial damage threshold.

The equivalent strain criterion of Simo and Ju [54] was used to define the energy norm of strain tensor as, $\tau :=\sqrt{2{\mathsf{\Psi }}_{\mathsf{\alpha }}}$ (i.e., the equivalent strain). The internal variable ${\mathrm{r}}_{\mathsf{\alpha }}$ can be integrated over time in the closed interval as:

$$r_{\alpha }\left(t\right)=\begin{array}{c}\max \\ \mathrm{T}\in \left[0,t\right]\end{array}\left[r_{\alpha }^0,\sqrt{2{\mathsf{\Psi }}_{\alpha }}\right]$$

(17)

In the equation above, $\mathrm{T}$ is time. The group of damage initiation criteria at any point of loading is given as:

$$\sqrt{2{\mathsf{\Psi }}_{\alpha }}-r_{\alpha }\le 0$$

(18)

$q_{\alpha }$ is defined using the following hardening rule:

$$\dot{q_{\alpha }}=H_{\alpha }\left(r_{\alpha }\right)\dot{r_{\alpha }},H_{\alpha }\left(r_{\alpha }\right)=\frac{\partial q_{\alpha }}{\partial r_{\alpha }},q_{\alpha }^0=r_{\alpha }^0$$

(19)

where $H_{\alpha }$ is the hardening/softening modulus given as:

$$H_{\alpha }\left(q_{\alpha }\left(t\right),h\right)=-\frac{{\left(q_{\alpha }^0\right)}^{2-\chi }}{\left(2-\chi \right)}\frac{1}{G_{\alpha }^f}{q}_{\alpha }^{\chi }\left(t\right)h$$

(20)

Here $\chi$ is a parameter that defines the rate of softening, material parameter $G_{\alpha }^f$ specifies the surface density of dissipated energy, and $h$ is the finite element characteristic size, which makes the model mesh dependent.

Blanco et al. [64] derived the macroscopic damage model for soft tissues using mesoscopic parameters. These macroscopic parameters include the constitutive model parameters of the collagen fibres, namely $k_1$ and $k_2$, and inelastic parameters such as elastic stress threshold ${\sigma }_f^u$ and surface density of dissipated energy ${𝒤}_f^u$, wherein the elastic stretch threshold were defined using wavy fibril nature into consideration. The geometrical properties of mesoscopic fibril are modelled as staggered arrays of tropocollagen molecules using a two-dimensional Hodge–Petruska model. The mechanical characteristics of each fibril constituent were established by identifying modes of failures and their associated weak planes. The macroscopic behaviour was obtained by homogenizing the obtained properties of fibrils and proteoglycan rich matrix. The hierarchical soft tissue model considered by Blanco et al. [64] is shown in Figure 5.

Comellas et al. [65] developed a generalized damage model for quasi-incompressible hyperelasticity in a total Lagrangian finite-strain framework. The deviatoric part of the hyperelastic constitutive model was incorporated with the Kachanov-like reduction factor (1–d), which defines the softening effects. The damage model was implemented on Neo-Hookean and Ogden hyperelastic models. It adopts additively decomposed volumetric and deviatoric parts of Helmholtz free energy.

$$\mathsf{\Psi }\left(\mathit{C},D\right)={\mathsf{\Psi }}_{vol}\left(J\right)+\left(1-d\right){\mathsf{\Psi }}_{dev}\left(\overline{\mathit{C}}\right)$$

(21)

The finite-strain based Kachanov effective stress using second Piola-Kirchhoff stress tensor ($\mathit{S}$) is given as:

$$\mathit{S}={\mathit{S}}_{vol}+\left(1-d\right){\overline{\mathit{S}}}_{dev}\mathrm{with}{\mathit{S}}_{vol}=-pJ{\mathit{C}}^{-1}\mathrm{and}{\overline{\mathit{S}}}_{dev}=2\frac{\partial {\mathsf{\Psi }}_{dev}}{\partial C}$$

(22)

Here $p$ is the hydrostatic pressure given by $p=\frac{\partial {\mathsf{\Psi }}_{\mathrm{vol}}}{\partial \mathrm{J}}$.

The damaged surface $\mathcal{F}=G\left(\tau \right)-G\left({\tau }^{max}\right)=0$ determine the limits of the initiating point of non-linear behaviours. This model allows the usage of different energy-based norms. The damage evolution law $G\left(\tau \right)$ is defined in terms of the norm, and $G\left({\tau }^{max}\right)$ is a scalar function where ${\tau }^{max}$ is the damage threshold. The criterion of Simo and Ju [54] was used to define the strain energy norm as $\tau =\sqrt{2{\mathsf{\Psi }}_{dev}}$. They considered two types of damage evolution laws: (1) linear and (2) exponential forms. The damage variable $d$ is expressed:

for linear softening,

$$d=G\left(\tau \right)=\frac{1-\frac{S_0^d}{\tau }}{1+H}$$

(23)

for exponential softening,

$$d=G\left(\tau \right)=1-\frac{S_0^d}{\tau }\exp \left[A\left(1-\frac{\tau }{S_0^d}\right)\right]$$

(24)

where

$$H=\frac{-{\left(S_{dev}^d\right)}^2}{2g_f^d}, A={\left[\frac{g_f^d}{{\left(S_0^d\right)}^2}-\frac{1}{2}\right]}^{-1}$$

(25)

Here $H$ and $A$ are the dissipation parameters for linear and exponential softening, respectively. $S_0^d$ is the basic damage threshold stress and $g_f^d$ represents rupture energy per unit volume.

The energy norm differentiation of the evolution law is essential for evaluating the constitutive tangent tensor; the differentials are given as:

for linear softening,

$$\frac{\partial G\left(\tau \right)}{\partial \tau }=\frac{-S_0^d}{{\tau }^2\left(1+H\right)}$$

(26)

for exponential softening,

$$\frac{\partial G\left(\tau \right)}{\partial \tau }=\frac{S_0^d+A\tau }{{\tau }^2}\exp \left[A\left(1-\frac{\tau }{S_0^d}\right)\right]$$

(27)

Using the decoupled Helmholtz free energy Equation (21) along with Equation (27) with $d=G\left(\tau \right)$, the damage incorporated and additively decomposed material elastic tangent constitutive are given as:

$$\begin{gathered} {ℂ}^{tan}={ℂ}_{vol}^{tan}+{ℂ}_{dev}^{tan} \\ with \left\{\begin{array}{c}{ℂ}_{vol}^{tan}=2\frac{\partial {\mathit{S}}_{vol}}{\partial \mathit{C}}=2p\frac{\partial \left(J{\mathit{C}}^{-1}\right)}{\partial \mathit{C}}+2J{\mathit{C}}^{-1}\otimes p\frac{\partial p}{\partial \mathit{C}}\\ {ℂ}_{dev}^{tan}=2\frac{\partial }{\partial \mathit{C}}\left[\left(1-d\right){\overline{\mathit{S}}}_{dev}\right]=\left(1-d\right){ℂ}_{0dev}^{tan}-\frac{\partial G\left(\tau \right)}{\partial \tau }{\overline{\mathit{S}}}_{dev}\otimes {\overline{\mathit{S}}}_{dev}\end{array}\right. \end{gathered}$$

(28)

where ${\widetilde{ℂ}}_{0dev}^{tan}=\left(2,\partial ,{\mathit{S}}_{\_dev}^{\circ }\right)/\partial \mathit{C}$ is deviatoric part of the damage free hyperelastic model, and $\partial G\left(\tau \right)/\partial \tau$ is damage dissipation defined in Equations (26) and (27) for linear and exponential softening, respectively. The damage parameters are defined using Matlab curve-fitting on experimental data of Martins et al. [72]. The quasi-incompressible volumetric part of the tangent modulus is unaffected by the damage, while the deviatoric part is influenced by the induced damage [65]. This model was applied for the rectus sheath to reproduce the experimental results numerically using their in-house finite element code (PLCd).

Polindara et al. [66] developed a damage model in line with Waffenschmidt et al. [73] to simulate balloon angioplasty in residually stressed blood vessels [66,74]. The model incorporates an anisotropic hyperelastic constitutive model defined by HGO [25] in the strain energy function for inclusion of the damage is given as:

$$\mathsf{\Psi }={\mathsf{\Psi }}_{vol}+{\mathsf{\Psi }}_{dev}+f_d\left(\chi \right)\left({\mathsf{\Psi }}_{df1}+{\mathsf{\Psi }}_{df2}\right)$$

(29)

To account for the fibre softening, a simple exponential damage function was introduced and is given by:

$$f_d=\exp \left({\eta }_d\left[{\chi }_d-\chi \right]\right)$$

(30)

where the evolution of local damage is controlled by the rate of change of damage variable $\chi$, ${\chi }_d$ is the damage threshold that initiates the evolution, and ${\eta }_d$ is an exponential saturation parameter that controls the rate of damage evolution. Polindara et al. [66] have assumed damage parameters based on the stiffness degressive behaviour found in the inflation test experiments by Holzapfel [75]. Using an 8-noded hexahedral Q1Q1P0 element in Abaqus, they implemented the above damage model through UEL for non-local gradient-enhancement. Further, Polindara et al. [66] have used the same methodology and extended it to incompressible material with three damage variables [74]. These damage variables evolve independently from each other, accounting for damage in the matrix and the two fibre families.

Ferreira et al. [67] provided a general framework for inducing damage in hyperelastic materials. The computational framework to locally model the anisotropic damage is considered in the non-linear geometry. This model assumes that the stretch patterns in the soft tissues result in pathological conditions. Further, they cause the stable degradation of the collagen fibres and the ground matrix of the soft tissue. The fully anisotropic hyperelastic material in the form of the strain energy density was defined as

$$\mathsf{\Psi }={\mathsf{\Psi }}_{vol}\left(J\right)+\left(1-d_{dev}\right){\mathsf{\Psi }}_{dev}\left(\overline{\mathit{F}}\right)+\left(1-d_{ani}\right){\mathsf{\Psi }}_{ani}+{\mathsf{\Phi }}_{dev}\left(d_{dev}\right)+{\mathsf{\Phi }}_{ani}\left(d_{ani}\right)$$

(31)

where $d_{dev},d_{ani}\in \left[0,1\right]$, are the internal damage variables for matrix and fibres, and functions ${\mathsf{\Phi }}_{dev}$ and ${\mathsf{\Phi }}_{ani}$ represent the damage propagation in the material before the tear propagation in matrix and fibres, respectively.

The proposed model adopts the Cauchy stress and effective tangent moduli tensor by an additive composition of each contribution.

$$\mathit{\sigma }={\mathit{\sigma }}_{vol}+\left(1-d_{dev}\right){\mathit{\sigma }}_{dev}+\left(1-d_{ani}\right){\mathit{\sigma }}_{ani}$$

(32)

$$\mathit{c}={\mathit{c}}_{jr}+{\mathit{c}}_{vol}+{\mathit{c}}_{dev}+{\mathit{c}}_{ani}$$

(33)

The parameters $d$ for the damage evolution for matrix and fibres are represented by an irreversible equation. The damage parameters defined in terms of the reduction factor is given as:

$$d=1-\overline{\mathrm{g}}$$

(34)

The reduction factor is obtained in terms of the equivalent strain ${\Xi }_s$ at time $s$ as defined by Simo [68]. The maximum value evolved during the deformation history till the current time is used for the evaluation of the reduction factor $\overline{\mathrm{g}}$.

$${\Xi }_{\mathrm{t}}=\begin{array}{c}max\\ s\in \left(0,t\right)\end{array}\sqrt{2\overline{U}}$$

(35)

In material degradation, the law is given by:

$$\overline{\mathrm{g}}=\left\{\begin{array}{c}1, {\Xi }_t<{\Xi }_{min}\\ \frac{1-\exp \left(\beta \left({\Xi }_{\mathrm{t}}-{\Xi }_{\max }\right)\right)}{1-\exp \left(\beta \left({\Xi }_{\min }-{\Xi }_{\mathrm{t}}\right)\right)}, {\Xi }_{\min }<{\Xi }_t<{\Xi }_{max}\\ 0, {\Xi }_t>{\Xi }_{max}\end{array}\right.$$

(36)

The thresholds ${\Xi }_{\min }$ and ${\Xi }_{max}$ are defined based on the tensile experiments. In the above equation, $\beta$ is the averaging operator that takes a bell-shaped form and is expressed as:

$$\beta \left(x,\xi \right)=\frac{1}{\mathsf{\Omega }\left(x\right)}{\left(1-\frac{{(x-\xi )}^2}{l_r^2}\right)}^2$$

(37)

Here $x$ is the point in the material, $\mathsf{\Omega }$ is the defined finite volume containing a set of points $\xi$ and $l_r$ is the characteristic length.

The developed damage formulation was tested for hyperelastic constitutive laws such as Neo-Hookean, Moony–Rivlin, Ogden, Humphrey–Martins, and HGO [67]. Wherein the damage parameters are set by admitting the degradation of matrix and fibres. The finite element simulations for internally pressurized healthy artery is implemented using various Abaqus user subroutines.

Rausch et al. [68] developed a soft tissue damage model by combining the continuum damage theory with smoothened particle hydrodynamics (SPH). The material response was assumed to be non-linear hyperelastic with a single family of fibres. The considered strain energy function for the fibre orientation was defined by HGO [25], as follows

$${\mathsf{\Psi }}_0=\frac{\mu }{2}\left(I_1-3\right)-K\left(J-1\right)+\frac{k_1}{4k_2}\left[\exp \left(k_2{\left(I_4-1\right)}^2\right)-1\right]$$

(38)

The damage initiation was based on the formulation of Simo [76]. ${\mathsf{\Psi }}_0$ is the strain energy free from damage $d\in \left[0,1\right]$. The concentration of $d$ increases with the increase in the maximum principal stretch in the material beyond the critical stretch ${\lambda }_{crit}$.

$$\mathsf{\Psi }=\left(1-d\right){\mathsf{\Psi }}_0$$

(39)

An exponential evolution equation is considered for the damage function,

$$d=\left\{\begin{array}{c}\exp \left({\left(\frac{\left({\lambda }_m-{\lambda }_{crit}\right)}{\tau }\right)}^2\right)-1,{\lambda }_m\ge {\lambda }_{crit}\\ 0, {\lambda }_m<{\lambda }_{crit}\end{array}\right.$$

(40)

The maximum principal stretch (${\lambda }_m$) from the previous increment is evaluated as ${\lambda }_m\left(t\right)=\underset{s\in \left[-\infty ,t\right]}{\max }\lambda \left(s\right)$ at a rate determined through $\tau$. Damage in the material initiates once ${\lambda }_m$ exceeds ${\lambda }_{crit}$ and it is irreversible. These damage parameters are defined from sensitivity analysis of damage parameters on material behaviour. The study of Rausch et al. [68] is aimed to explore the applicability of the Normalized Total Lagrangian-SPH method in soft tissue mechanics. The Rausch et al. [68] damage model was used to simulate material discontinuities in a single edge notch soft tissue specimen under uniaxial tension. The SPH simulations are carried using MATLAB, and for validating the SPH results, similar FEM simulations are carried using open-source FEbio.

Fathi et al. [69] had implemented a non-local integral-damage model to overcome the numerical artefacts such as mesh dependency and spurious localization in soft tissues. This damage model, which considers large deformation, was used to simulate aortic dissections. The strain energy equation, consisting of the damage variables, is given by:

$$\begin{gathered} \mathsf{\Psi }=K\ln J^2+\left(1-d_m\right)\left\{c_1\left(1-{\overline{I}}_1\right)+c_2\left(1-{\overline{I}}_2\right)\right\} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left(1-d_f\right)\frac{k_1}{2k_2}\left\{\exp \left(k_2\kappa \langle {\overline{I}}_1+\left(1-3\kappa \right){\overline{I}}_{fi}-I_{0i}{\rangle }^2\right)-1\right\}, i=1,2 \end{gathered}$$

(41)

where $c_1$, and $c_2$, are material parameters, and ${\overline{I}}_{f1}={\overline{I}}_4$ and ${\overline{I}}_{f2}={\overline{I}}_6$ are deviatoric invariants. The Macaulay bracket ($\langle ■\rangle =\left(\left|■\right|+■\right) /2$) assumes the fibres can support tension only. The state of folding fibres and their mobilization is defined using dimensionless parameter $I_{0i}$.

The damage function is defined as:

$$\tilde{G}\left({\Xi }_{\mathrm{t}}^{\mathrm{k}}\right)=1-d_k=\left\{\begin{array}{cc}1& if {\Xi }_{\mathrm{t}}^{\mathrm{k}}<{\Xi }_{min}^{\mathrm{k}} \\ \frac{{\Xi }_{\max }^{\mathrm{k}}-{\Xi }_{\mathrm{t}}^{\mathrm{k}}}{{\Xi }_{\max }^{\mathrm{k}}-{\Xi }_{\min }^{\mathrm{k}}}\exp \left[{\beta }_k\left(\frac{{\Xi }_{\max }^{\mathrm{k}}-{\Xi }_{\mathrm{t}}^{\mathrm{k}}}{{\Xi }_{\max }^{\mathrm{k}}-{\Xi }_{\min }^{\mathrm{k}}}\right)\right]& if {\Xi }_{min}^{\mathrm{k}}\le {\Xi }_{\mathrm{t}}^{\mathrm{k}}\le {\Xi }_{max}^{\mathrm{k}}\\ 0& if {\Xi }_{\mathrm{t}}^{\mathrm{k}}>{\Xi }_{max}^{\mathrm{k}}\end{array}\right.$$

(42)

where the total and initial equivalent strain are defined using ${\Xi }_{max}^{\mathrm{k}}$ and ${\Xi }_{min}^{\mathrm{k}}$, respectively and the exponential coefficient ${\beta }_k$ represents the damage saturation. These damage parameters are defined based on Martins et al. [72] for the abdomen rectus sheath and Weiss et al. [77] for the ligament.

The strain evolution parameter is defined in similar lines with Simo and Ju [55].

$${\Xi }_{\mathrm{t}}^{\mathrm{k}}=\begin{array}{c}\max \\ s\in \left(-\infty ,t\right)\end{array}\sqrt{2{\mathsf{\Psi }}_0^{\mathrm{k}}\left(s\right) }; k=m,f_1,f_2$$

(43)

where the subscript $0$ represents the intact material, and the superscript $k$ represents the particular component of the strain energy, i.e., $m$, $f_1$ and $f_2$ corresponds to the matrix, and the first and second family of fibres, respectively.

The damage model proposed by Fathi et al. [69] was simulated with an integral-type non-local scheme that can be implemented on geometrically nonlinear soft tissue by overcoming mesh dependency. The role of the collagen fibres is also studied with the damage model and was found to be predominant. The mechanical response of the soft tissue varies with the position and the direction of the fibre. The same study also proposed soft tissue tear simulations with the combination of damage model and XFEM or meshless methods.

An anisotropic multi-physics damage model to capture damage initiation and propagation in annulus fibrosus was proposed by Gao et al. [70]. This model was derived based on the damage model proposed by Balzini et al. [37] by integrating continuum mixture theory and CDM. The constitutive strain energy equation is given by:

$${\mathsf{\Psi }}_{\mathrm{AF}}={\mathsf{\Psi }}_m+{\displaystyle \sum }_{f=4,6}{\mathsf{\Psi }}_f$$

(44)

with

$${\mathsf{\Psi }}_m=\frac{\mu }{2}\left(I_1-3\right)-\mu \ln J+\lambda {\left({\varphi }_0^w\right)}^2\left[\frac{J-1}{{\varphi }_0^w}-\ln \left(\frac{J-1}{{\varphi }_0^w}+1\right)\right]$$

(45)

and

$${\mathsf{\Psi }}_f=\frac{k_1}{2k_2}\left\{H\left(I_f-1\right)\exp \left[I_f-1\right]-1\right\}, f=4,6$$

(46)

where the intact strain energy density function of the matrix is given by ${\mathsf{\Psi }}_{\mathrm{m}}$ and fibre is given by ${\mathsf{\Psi }}_f$. ${\varphi }_0^w$ represents the reference state water volume fraction, $\lambda$ is material constants, $I_1$ is the first principal invariant of $\mathit{C}$, and $H$ is the Heaviside step function.

After introducing the damage parameters, the strain energy function becomes:

$${\mathsf{\Psi }}_{\mathrm{AF}}=\left(1-d_m\right){\mathsf{\Psi }}_m+{\displaystyle \sum }_{f=4, 6}\left(1-d_f\right){\mathsf{\Psi }}_f$$

(47)

In the reported study [71], the damage in the matrix was neglected, and only the damage contribution of the fibre was solely considered. The damage evolution variable was defined as:

$$d_f=d_{\max }\left[1-\exp \left(-\frac{{\beta }_i}{\gamma }\right)\right]$$

(48)

$d_{\max }$ and $\gamma$ is the damage parameters and the internal damage variable ${\beta }_i$ defined in terms of function ${\mathsf{\Psi }}_f$ as:

$${\beta }_i=\begin{array}{c}\max \\ 0\le \tau \le t\end{array}{\mathsf{\Psi }}_f\left(\tau \right)-{\hat{\mathsf{\Psi }}}_f$$

(49)

where ${\hat{\mathsf{\Psi }}}_f$ is the limiting strain energy density of fibre beyond which damage occurs. Damage parameters are approximated from the experimental data of Pezowicz [78].

The study of Gao et al. [69] was aimed to develop a multi-physics damage model to study damage in annulus fibrosus. The continuum mixture theory was used to develop a coupled field problem that considers the water content along with the soft tissue material. While CDM is used to model the damage in the coupled field problem. Numerical studies of annulus fibrosis for damage under compression found the results are sensitive to ${\hat{\mathsf{\Psi }}}_f$ then the other damage parameters.

A damage model to predict damage growth before the rupture in ascending thoracic aneurysms was reported by Mousavi et al. [71]. The proposed model is layer-specific and uses a constrained mixture theory that inherently considers internal/residual stresses. As per constrained mixture theory, distinct strain energy density was proposed for every individual constituent concerned with the contribution of its mass function. The damage was assumed to occur in fibre constituents, i.e., elastin and collagen, and the strain energy function was given as:

$$\mathsf{\Psi }=2{\mathsf{\Psi }}_{vol}+{\rho }_m{\mathsf{\Psi }}_m+\left(1-d_e\right){\rho }_{e}2{\mathsf{\Psi }}_{dev}+{\displaystyle \sum }_{i=1}^n\left(1-d_{ci}\right){\rho }_{ci}{\mathsf{\Psi }}_{ci}$$

(50)

where subscripts $e$, $ci,$ and $m$ represent elastin fibre, collagen fibres, and smooth muscle cells, respectively; $d_j$ are respective damage parameters and ${\rho }_j$ are specific mass fractions. The constitutive strain energy density in tension and compression for collagen and smooth muscle cells, distinctly defined as:

$${\mathsf{\Psi }}_{ci}=\frac{k_1^{ci,t}}{2k_2^{ci,t}}\left\{\exp \left[k_2^{ci,t}{\left({\overline{I}}_4^{ci}-1\right)}^2\right]-1\right\},\mathrm{under}\mathrm{tension}$$

(51)

$${\mathsf{\Psi }}_{ci}=\frac{k_1^{ci,c}}{2k_2^{ci,c}}\left\{\exp \left[k_2^{ci,c}{\left({\overline{I}}_4^{ci}-1\right)}^2\right]-1\right\},\mathrm{under}\mathrm{compression}$$

(52)

and

$${\mathsf{\Psi }}_m=\frac{k_1^{m,t}}{2k_2^{m,t}}\left\{\exp \left[k_2^{m,t}{\left({\overline{I}}_4^m-1\right)}^2\right]-1\right\},\mathrm{under}\mathrm{tension}$$

(53)

$${\mathsf{\Psi }}_m=\frac{k_1^{m,c}}{2k_2^{m,c}}\left\{\exp \left[k_2^{m,c}{\left({\overline{I}}_4^m-1\right)}^2\right]-1\right\},\mathrm{under}\mathrm{compression}$$

(54)

A linear softening law is used to evaluate the damage variable, as described by Comellas et al. [65]. The apparent density of damage is assumed to evolve with mechanical loading.

$$d_j=G_j\left({\psi }_j\right)=\frac{1-\frac{{\psi }_j^0}{{\psi }_j}}{1+H_j}$$

(55)

here

$$H_j=-\frac{{\left({\psi }_j^0\right)}^2}{2{\omega }_j}, {\omega }_j=\frac{{\mathsf{\Omega }}_j}{L_o}\mathrm{and}{\psi }_j=\sqrt{2{\mathsf{\Psi }}_j}$$

(56)

where $j\in \left\{e,ci\right\}$, ${\psi }_j^0$ is the basic damage threshold and ${\omega }_j$ represents the rupture energy per unit volume that is defined based on the maximum fracture energy dissipated per unit area, ${\mathsf{\Omega }}_j$, and $L_o$ characteristic element length. Buldge inflation experiments were performed, and with the help of curve fitting to experimental data, damage parameters are defined. Numerical simulations of the developed damage model were carried in a commercial finite element software (Abaqus) using a user-defined subroutine (UMAT). Human aortic aneurysm specimens were also simulated under uniaxial tension, patient-specific over-pressurization, and bulge inflation tests.

A damage model to accurately capture the phenomena of passive damage in arteries was proposed by Ghasemi et al. [35]. These passive damage phenomena include Mullins effects, hysteresis, permanent set, and the effect is captured up to the level of rupture. In the model of Ghasemi et al. [35], the damage is assumed to take place in elastic fibre and collagen fibre only. The constitutive equation with damage variables incorporated is defined as:

$$\begin{gathered} \mathsf{\Psi }={\mathsf{\Psi }}_{vol}+{\mathsf{\Psi }}_{dev}+{\displaystyle \sum }_{M_{ef}=M_{4ef},M_{6ef}}\frac{k_1^{ef}}{2k_2^{ef}}\left\{\exp \left[k_2^{ef}\left(1-d^{ef}\right){\left({\overline{I}}_{M_{ef}}-1\right)}^2\right]-1\right\} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum }_{M_{cf}=M_{4cf},M_{6cf}}\frac{k_1^{cf}}{2k_2^{cf}}\left\{\exp \left[k_2^{cf}\left(1-d^{cf}\right){\left(\kappa {\overline{I}}_1+\left(1-3\kappa \right){\overline{I}}_{M_{cf}}-1\right)}^2\right]-1\right\} \end{gathered}$$

(57)

where material parameters $k_1^{ef}$& $k_2^{ef}$ and $k_1^{cf}\& k_2^{cf}$ corresponds to elastin fibre and collagen fibre, respectively. The deviatoric invariants ${\overline{I}}_{M_{ef}}$ and ${\overline{I}}_{M_{cf}}$ represent the square of stretch in elastin fibre and collagen fibre, respectively. The damage of elastin and collagen fibres are represented by $d^{ef}$ and $d^{cf}$, respectively.

The damage functions capture the continuous and discontinuous softening using two internal variables ${\beta }_i$ and ${\gamma }_i$ respectively, and defined as:

$${\beta }_i=\langle {\tilde{\beta }}_i-{\tilde{\beta }}_i^{ini}\rangle$$

(58)

Here $i=ef$ for elastin fibres and $i=cf$ for collagen fibres, ${\tilde{\beta }}_i^{ini}$ the initial parameter that ensures the damage evolution. Macaulay brackets, ($\langle ■\rangle =\left(\left|■\right|+■\right) /2$) ensures only positive values. ${\tilde{\beta }}_i$ is constructed as per the changes in the pseudo-invariant over the complete loading history, as:

$${\tilde{\beta }}_i={{\displaystyle \int }}_0^t\langle {\overline{I}}_{M_i}^{*}\rangle ds$$

(59)

where ${\gamma }_i$ is defined based on the maximum value of ${\overline{I}}_{M_i}^{*}$ evolved during the loading history till the current loading state, and it is given as:

$${\gamma }_i={}_{S\in \left[0, S\right]}{}^{max}\overline{I}\langle {\overline{I}}_{M_i}^{*}-{\overline{I}}_{M_i}^{*ini}\rangle$$

(60)

The damage variable is constructed as:

$$d^i=d_{\infty }^i\left[1-\exp \left(-\frac{{\gamma }_i}{{\gamma }_{{\infty }_i}}\right)\right]\left[1-\exp \left(-\frac{{\beta }_i}{{\beta }_{S_i}}\right)\right], d^i\in \left[0,1\right]$$

(61)

where the predefined variable $d_{\infty }^i$ limits the overall damage in elastin fibres and collagen fibres, while the parameters ${\beta }_{S_i}$ and ${\gamma }_{{\infty }_i}$ represent the continuous and discontinuous softening of the soft tissues, respectively. Ghasemi et al. [35] have used an inverse FE algorithm to draw damage parameters using experimental results of uniaxial tension tests.

An anisotropic microsphere-based approach to model the damage in soft vascular tissue was developed by Saez et al. [79]. A microsphere-based damage approach for modelling damage by considering the microstructure is initially developed by Miehe et al. [80] and Dal and Kaliske [81] for rubbers. Saez et al. [79] have extended the model for anisotropic soft tissues by neglecting fibre crosslinks and sliding between fibres and the surrounding matrix. However, the fit between the experimental stress-stretch data and that predicted by the microstructural model of Saez et al. [79] was reported to be not satisfactory by Pena et al. [82]. Particularly, the correlation between the experimental data and the prediction by the microstructural damage model was found to be worse, as compared to the phenomenological model by Pena [83]. Hence detailed discussion on the microstructural damage model by Saez et al. [79] is not included in this review.

3.2. Pseudo-Elasticity

In 2015, Pierce et al. [84] had modelled the material response of human thoracic and abdominal aortic tissues using the damage model proposed by Weisbecker et al. [41]. Pierce et al. damage model differentiated between physiological and supraphysiological loading. They compared the damage response in healthy and diseased vascular tissues using it. In particular, for the diseased tissues, abnormal aortic aneurysms were considered (abnormal swelling or bulge in the wall of the artery is known as an aortic aneurysm). The elastic damage model postulated in terms of the strain energy density is given as:

$$\mathsf{\Psi }={\mathsf{\Psi }}_{vol}+{\mathsf{\Psi }}_{dev}+{\displaystyle \sum }_{i=4,6}\left[{\eta }_{fi}\left({\mathsf{\Psi }}_{dfi}\right)+{\mathsf{\Phi }}_{fi}\left({\eta }_{fi}\right)\right]$$

(62)

where ${\mathsf{\Phi }}_{fi}$ is the smooth damage function of the fibres, and ${\eta }_{fi}\in \left[0, 1\right]$ represents the damage variables of the two fibre families given as:

$${\eta }_{fi}=1-\frac{1}{r_f}\mathrm{erf}\left[\frac{1}{m_f}\left({\mathsf{\Psi }}_{dfi}^{max}-{\mathsf{\Psi }}_{dfi}\right)\right]$$

(63)

where ${\mathsf{\Psi }}_{dfi}^{max}$ is the maximum strain energy evolved during the deformation history, $r_f>1$ defines the maximum allowable damage in the fibres subjected to loading, and $m_f>0$ defines the accumulation of softening in the fibres. The minimum of the damage variable characterizes the damage induced in the fibres as:

$${\eta }_{fi}^{min}=1-\frac{1}{r_f}\mathrm{erf}\left(\frac{1}{m_f}{\mathsf{\Psi }}_{dfi}^{max}\right)$$

(64)

The damage model presented includes the effect of damage on fibres only, and the damage parameters are defined using curve-fitting on uniaxial tension experimental data. Pierce et al. [84] study involve constitutive modelling for damage, damage experiments, and statistical data analysis to identify the material and damage parameters. A similar pseudo-elasticity model is used by He et al. [85] to simulate damage of the artery during stent deployment.

Holzapfel and Ogden [86] have proposed a progressive damage model for the collagen fibres based on the pseudo-elastic model. Their model considers both the Mullins effect and cross-linking between the collagen fibres. The same model is validated for an experimental behaviour of rat tail tendon fibres. The damage induced strain energy density is defined as:

$$\mathsf{\Psi }={\mathsf{\Psi }}_{dev}+\eta {\mathsf{\Psi }}_{af1}+\varphi \left(\eta \right)$$

(65)

where $\eta$ is a dimensionless damage variable that introduces the damaging effect and $\varphi \left(\eta \right)$ is some damage measure. The damage variable and damage function are given by

$$\eta =\exp \left(-\frac{{\mathsf{\Psi }}_{af1}-{\mathsf{\Psi }}_{af1}^c}{m}\right)$$

(66)

$$\varphi \left(\eta \right)=m\eta \log \eta +\left(1-\eta \right)\left[m+{\mathsf{\Psi }}_{af1}^c\right]$$

(67)

Here parameter $m>0$ has the same dimension as $\mathsf{\Psi }$ and $I_{4c}={\lambda }_c^2$ is the critical stretch value of $I_4$, which is responsible for the initiation of the damage, i.e., $\eta$ decreases from 1 for the stretch $\lambda >{\lambda }_c$ down to ${\eta }_f$ for a failure occurs at $\lambda ={\lambda }_f$.

The collagen fibres cross-links are included with the unit vectors ${\mathit{L}}^{+}$ and ${\mathit{L}}^{-}$ around the collagen fibre direction ${\mathit{a}}_1$. The unit vectors ${\mathit{L}}^{+}$ and ${\mathit{L}}^{-}$ represents the fibre cross-links are logically symmetric about ${\mathit{a}}_1,$ and the operation of the deformation gradient $\mathit{F}$ on them given as:

$${\mathit{L}}^{\pm }=\pm c_0{\mathit{a}}_1+s_0{\mathit{a}}_2$$

(68)

$$\mathit{F}{\mathit{L}}^{\pm }=\pm c_0\mathit{F}{\mathit{a}}_1+s_0\mathit{F}{\mathit{a}}_2$$

(69)

For conciseness, the representation of $s_0=\sin {\alpha }_0$ and $c_0=\cos {\alpha }_0$ were used by Holzapfel and Ogden [86], where ${\alpha }_0$ defines the relative orientation of fibre cross-link vectors (${\mathit{L}}^{+}$ and ${\mathit{L}}^{-}$) with reference to the direction of ${\mathit{a}}_1$ (Figure 6). To model the effect of the collagen cross-links, Holzapfel and Ogden [86] introduced a couple of pseudo invariants $I^{\pm }$ and $I_8^{\pm }$. Wherein $I^{\pm }$ represents the squares of the stretches in the cross-link directions and $I_8^{\pm }$ describes the coupling between the collagen fibre and cross-links.

$$I^{\pm }=c_0^2{I}_4\pm 2s_0{c}_0\left(\mathit{C}{\mathit{a}}_1\right)·{\mathit{a}}_2+s_0\left(\mathit{C}{\mathit{a}}_2\right)·{\mathit{a}}_2$$

(70)

$$I_8^{\pm }=\pm c_0{I}_4+s_0\left(\mathit{C}{\mathit{a}}_1\right)·{\mathit{a}}_2$$

(71)

where $\mathit{C}$ is the right Cauchy-Green tensor.

For uniaxial tension where the stretch $\lambda$ in the fibre direction ${\mathit{a}}_1$ gives $\mathit{F}{\mathit{a}}_1=\lambda {\mathit{a}}_1$ and $\mathit{F}{\mathit{a}}_2={\lambda }^{-\frac{1}{2}}{\mathit{a}}_2$. The deformation gradient acting on the cross-link vectors gives,

$$\mathit{F}{\mathit{L}}^{\pm }=\pm c_0\lambda {\mathit{a}}_1+s_0{\lambda }^{-\frac{1}{2}}{\mathit{a}}_2$$

(72)

and additionally, gives the cross-link directions and quantities as

$$I\equiv I^{\pm }=c_0^2{\lambda }^2+s_0^2{\lambda }^{-1}, I_8=I_8^{+}=c_0{\lambda }^2\mathrm{and}{I}_8^{-}=-I_8^{+}$$

(73)

The specific strain energy function with isotropic strain energy, anisotropic strain energy along with the quadratic terms correlating the cross-links and fibre/cross-link density interactions given as:

$$\mathsf{\Psi }={\mathsf{\Psi }}_{dev}+\eta {\mathsf{\Psi }}_{af1}^c+\frac{1}{2}\nu {\left(I-1\right)}^2+\frac{1}{2}\kappa {\left(I_8-c_0\right)}^2$$

(74)

where the stress-like parameters $\nu$ and $\kappa$ correlate the cross-links and interactions, respectively. Wherein $\nu$ represents the density of cross-links, and $\kappa$ represents the measure of the interaction energy. For instance, the damage variable for uniaxial tension is given by:

$$\eta =\exp \left[-\frac{k_1}{2m{k}_2}\left\{\exp \left[k_2{\left(\lambda -1\right)}^2\right]-\exp \left[k_2{\left({\lambda }_c-1\right)}^2\right]\right\}\right]$$

(75)

The Cauchy stress $\sigma$ becomes:

$$\begin{gathered} \sigma =\mu \left({\lambda }^2-{\lambda }^{-1}\right)+2k_1\eta {\lambda }^2\left({\lambda }^2-1\right)\exp \left[k_2{\left({\lambda }^2-1\right)}^2\right]+4\nu \left(I-1\right)\left(c_0^2{\lambda }^2-s_0^2{\lambda }^{-1}\right) \\ +4\kappa \left(I_8-c_0\right)c_0{\lambda }^2 \end{gathered}$$

(76)

The same model was applied for planar deformations for the case of simple shear, wherein both collagen fibres and cross-links are assumed to be lying in a plane [86,87]. Wherein the critical stretch (damage parameter) is defined using least square curve-fitting on uniaxial tension experimental data. The proposed model focuses on damage at the collagen fibre level and does not consider the fibrils and proteoglycans structure.

3.3. Hyperelastic Softening

An invariant-based constitutive model that accounts for damage for skin was developed by Li and Luo [63]. The skin was assumed to have two symmetric families of fibre, and their structures were constant across its depth. Their damage model was developed based on the HGO type strain energy function and the Volokh damage model [89,90]. The strain energy function that incorporates the damage is given by:

$$\mathsf{\Psi }=\frac{\mu }{2}\left[\left(I_1-3\right)-\frac{{\left(I_1-3\right)}^{m+1}}{\left(m+1\right){\left(\zeta -3\right)}^m}\right]+\frac{k_1}{k_2}\left\{\exp \left(k_2{A}^2\right)-1-\frac{2k_2{A}^{n+2}}{\left(n+2\right){\left({\xi }^2-1\right)}^n}\right\}$$

(77)

where $A={\lambda }_f^2-1$, with ${\lambda }_f=\sqrt{\kappa I_1+\left(1-3\kappa \right)I_4}$ is the fibre stretch, and $m, n, \zeta ,$ and $\xi$ are phenomenological damage parameters to induce softening. In particular, when damage occurs, $m$ defines the stretch curve sharpness, and $\zeta$ represents the value of $I_1$ associated with the matrix. $n$ is the equivalent contrary to $m,$ and $\xi$ represents the onset of damage of the fibre in terms of ${\lambda }_f$.

The softening hyperelasticity model is extended to human artery adventitia, where failure is considered as part of the constitutive model [88]. This model is developed using the HGO material model [34] and the energy limiters. The limiting value for the strain energy and failure energy is enforced using the energy limiters that restrict the stresses in the constitutive equations [90]. Volokh [88] proposed the strain energy function for adventitia as

$$\psi =\frac{\mu }{2} \left(I_1-3\right)H\left({\zeta }_4\right)H\left({\zeta }_6\right)+H\left({\zeta }_4\right)\left\{{\psi }_4^f-H\left({\zeta }_4\right){\psi }_4^e\right\}+H\left({\zeta }_6\right)\left\{{\psi }_6^f-H\left({\zeta }_6\right){\psi }_6^e\right\}$$

(78)

where the step function $H\left({\zeta }_i\right)$ is defined as $H\left({\zeta }_i\right)=0$ for ${\zeta }_i<0$, else $H\left({\zeta }_i\right)=1$, and strain energies ${\psi }_i^f$ & ${\psi }_i^e$ of collagen fibres are defined as:

$${\psi }_i^f={\mathsf{\Phi }}_i{m}_i^{-1}\mathsf{\Gamma }\left(m_i^{-1}, 0\right), {\psi }_i^e={\mathsf{\Phi }}_i{m}_i^{-1}\mathsf{\Gamma }\left(m_i^{-1}, W_i^{m_i}{\mathsf{\Phi }}_i^{-m_i}\right), i=4,6$$

(79)

Here the gamma function is defined as $\mathsf{\Gamma }\left(s,x\right)={{\displaystyle \int }}_x^{\infty }{a}^{s-1}\exp \left(-a\right)da$, and ${\mathsf{\Phi }}_i$, ${\mathrm{m}}_i$ material failure parameters. In particular, ${\mathsf{\Phi }}_i$ is the energy limiter that describes the average bond energy. The strain energy function $W_i$ for the aligned intact fibres is given as,

$$W_i={\mathsf{\Psi }}_{afi}=\frac{k_1}{2k_2}\left[\exp k_2{\left(I_i-1\right)}^2-1\right] , i=4,6$$

(80)

In this model, fibres contribute to strain energy only in tension, i.e., $I_4>1$ and $I_6>1$.

In Equation (78), ${\zeta }_i\in \left(-\infty , 0\right]$ is a switch parameter, and its evolution is defined as:

$$\dot{{\zeta }_i}=-H\left({ϵ}_i-\frac{{\psi }_i^e}{{\psi }_i^f}\right), {\zeta }_i\left(t=0\right)=0, i=4, 6$$

(81)

where ${ϵ}_i$ is a dimensionless precision constant that is defined as $0<{ϵ}_i\ll 1$.

In the proposed model, the material response is hyperelastic when ${\psi }_i^e<{\psi }_i^f$. The strain energy remains constant ($W_f={\psi }_i^f)$ and prevents healing in the material & enables energy dissipation. The switch parameters differentiate the elastic and damage response, i.e., ${\zeta }_i=0$ for elastic response and ${\zeta }_i<0$ for irreversible damage, and strain is dissipated. The step multipliers assume that the damage in either family of fibres results in whole tissue failure. Volokh has defined the damage parameters using least square curve-fitting on uniaxial tension experimental results [88].

4. Rupture Modelling

Propagation of crack in soft tissues is considered as tear propagation/rupture. Soft biological tissue exhibits a complex rupture phenomenon due to the presence of collagen fibres. These collagen fibres in the soft tissue make them resistant to defects [30]. Crimped collagen fibres in the vicinity of the crack tip gradually straighten, providing resistance to crack propagation [16,31]. A tear propagation approach should be capable of capturing such complex phenomena. Various rupture modelling approaches used for soft tissues are reviewed in the present section.

The above-discussed damage approaches operate within the standard continuum mechanics approach wherein the displacement fields are continuous. In contrast, soft tissue rupture is macroscopic damage, which is considered as the point where the crack initiates, i.e., where discontinuity in the material occurs. Such a macroscopic damage phenomenon is dealt with by the fracture mechanics approach [55,91]. Fracture mechanics (FM) involves discontinuous displacement fields and uses techniques such as extended finite element method (XFEM), cohesive zone modelling (CZM), crack phase-field modelling (CPFM), etc. In rupture modelling, XFEM and CZM are the numerical techniques to simulate fracture, while CPFM is a mathematical model. The initiation and propagation criteria need to be defined for modelling the propagating crack. The modelling of the crack propagation in the conventional finite element method (FEM) requires re-meshing, which is a cumbersome task.

In contrast, XFEM can model discontinuities and their propagation by overcoming the problem of re-meshing [92]. XFEM is an extension of FEM, which is based on the concept of partition of unity that introduces the enrichment functions associated with additional degrees of freedom [93]. The Heaviside functions and asymptotic crack-tip fields represent the enriched discontinuous displacement fields and crack-tip singularity [94]. The XFEM has been extensively used for surgical cutting simulations [95,96,97], rupture simulations of soft-hydrated tissues [98,99], and arterial dissections [100,101].

In CZM, a cohesive surface is placed in the intact region of the material where the crack propagates, as shown in Figure 7. The cohesive surface is modelled with special elements called cohesive elements. The crack propagation is modelled with the help of traction separation law, i.e., when the opening displacement reaches a limiting value, the traction along the surface disappears [102,103]. Pandolfi and co-workers have extensively applied CZM on arteries, also proposed an anisotropic extension to irreversible cohesive law [104,105,106]. Gasser and Holzapfel [107] have simulated arterial strip dissection by adopting the modality of combined CZM, and XFEM approaches introduced by Moes and Belytchko [108].

In contrast to XFEM and CZM, CPFM models discontinuity as a separate continuous field, and it is coupled with the deformation field to model the crack propagation. The minimum energy variational principle is used for numerical analysis to solve the coupled field. In line with XFEM, a damage initiation criterion is needed for the crack propagation, i.e., the crack phase-field evolution initiates. Miehe and co-workers made a phenomenal development of CPFM with the thermodynamically consistent and algorithmically robust formulation [109,110]. Gultekin introduced the application of CPFM in biomechanics, which was later used for simulating aortic dissections [48,111,112,113]. Raina and Miehe [114] proposed an anisotropic failure criterion for computing driving forces during damage growth in soft biological tissues. Furthermore, numerical aspects associated with aortic dissections were investigated in the studies by Raina and Miehe [114] and Gultekin et al. [48,112].

4.1. Extended Finite Element Method (XFEM)

In conventional FEM, the displacement field is interpolated using shape functions ($N^I$) and nodal degrees of freedom (${\underset{\_}{\mathit{u}}}_c^I$). Additionally, in XFEM, to model the crack and its propagation, the displacement field is incorporated with Heaviside function ($\mathcal{H}$) and enriched degrees of freedom (${\underset{\_}{\mathit{u}}}_e^I$) as follows,

$$\underset{\_}{\mathit{u}}={\displaystyle \sum }_{I=1}^{n_{el}}{N}^I{\underset{\_}{\mathit{u}}}_c^I+\mathcal{H}{\displaystyle \sum }_{I=1}^{n_{el}}{N}^I{\underset{\_}{\mathit{u}}}_e^I+{\displaystyle \sum }_{I=1}^{n_{el}}{\displaystyle \sum }_{\alpha =1}^4{\mathit{F}}_{\alpha }{N}^I{\underset{\_}{\mathit{u}}}_{\alpha }^I$$

(82)

where the index ($I$) runs from $1$ to $n_{el}$ (number of nodes per element). ${\mathit{F}}_{\mathit{\alpha }}$ asymptotic crack tip functions and ${\underset{\_}{\mathit{u}}}_{\alpha }^I$ nodal enriched degrees of freedom. The third term in the displacement characterizes the stress singularity at the crack tip. The Heaviside function ($\mathcal{H}$) is used to model strong discontinuities such as cracks; in other words, it represents the displacement jump at the crack surface. The asymptotic crack tip functions (${\mathit{F}}_{\alpha }$) represent the singularity near the crack tip while ${\underset{\_}{\mathit{u}}}_{\alpha }^I$ provides an additional degree of freedom to model and compute the crack propagation.

The partition of unity is an important characteristic of XFEM. As mentioned earlier, a damage model is required to model crack propagation using XFEM. In XFEM, the damage is modelled using either a cohesive law or fracture mechanics approach; both approaches are discussed below.

4.1.1. XFEM Using Fracture Mechanics

Wang et al. [115] have developed a model for tear propagation in two-dimensional arteries using the energy-based approach, in which the linear elasticity based Griffith energy balance principle is extended to fibre-reinforced materials. The crack propagation criteria are based on the energy release rate $G$ (ERR). It is defined as the variation in the total potential energy per unit propagation of the tear. Wang et al. [115] have numerically calculated ERR based on the variation in the global energy, where the numerical approximation of $G$ is defined using the potential energy $\Pi$ and change in the crack length $\delta a$ as:

$$G=-\frac{\delta \Pi }{\delta a}=-\frac{{\Pi }_{a+\delta a}-{\Pi }_a}{\delta a}$$

(83)

Here, potential energy $\Pi \left(\mathrm{a}\right)={\mathrm{U}}_e-W$, $U_e$ is the equilibrium strain energy of the tissue, and $W$ represents the work done due to load. The anisotropic hyperelastic tissue response is incorporated using strain energy density defined by the HGO model [34]. To obtain $\Pi \left(a\right)$, a boundary value problem is solved numerically for increasing crack lengths at intermittent points and energy values so obtained are interpolated with a cubic spline polynomial. This interpolation smoothly approximates $\Pi \left(a\right)$, which can be used to estimate ERR. The obtained values of ERR are used in simulating the two-dimensional crack propagation in arteries. This model ignores the plastic effects of tear propagation.

Karimi et al. [116] have modelled the initiation and propagation of a crack in the coronary artery to study the relation between rupture in the coronary artery and atherosclerosis. To achieve it, crack initiation and propagation in the healthy and atherosclerotic human coronary arterial walls were simulated with cracks placed circumferentially along the luminal and in the radial direction. Their model makes use of the virtual crack extension method (VCE) of XFEM and elastoplastic fracture mechanics criteria for crack initiation and propagation. In particular, they used an energy release rate criteria (J-integral) using nonlinear fracture mechanics of LS-DYNA [117]. The J-integral using the VCE is defined using a continuous weighting function $q$. It is defined on the surface of the body as zero and unity on the crack front nodes and interior of the body surface. In VCE, J-integral is given as:

$$J=-\frac{1}{\Delta A_c}{{\displaystyle \int }}_A^{ }\left[\left(W{\delta }_{1i}-{\sigma }_{ij}{u}_{j,1}\right)q_{,i}+{\sigma }_{ij}{\epsilon }_{ij,1}^{0}q\right]dA$$

(84)

where $\Delta A_c$ is the virtual increase in crack area, $W$ is the strain energy density, ${\delta }_{1i}$ is Kronecker delta, ${\sigma }_{ij}$ is stress in the defined area, ${\epsilon }_{ij,1}^0$ is the initial strain, and $u_i$ is the displacements in the area of interest. The yield parameters for rupture simulation were drawn from the uniaxial tension experiments. Their study made use of the standard library of LS-DYNA, which was originally developed by Lindström et al. [117]. The material constitutive response for the healthy and atherosclerotic coronary arteries is assumed to be linear elastic.

4.1.2. XFEM Using Cohesive Law

In XFEM, a traction law defines the criterion for crack initiation and propagation using cohesive law. The studies of Jayendran and Ruimi [118] and Wang et al. [119] used the linear traction separation law developed by Ferrara and Pandolfi [105]. Both studies simulate crack propagation in arteries, whereas the Wang et al. [119] study dealt with crack propagation in the residually stressed artery. The crack propagation is governed using linear traction separation law given as:

$$G_c=\frac{1}{2}{T}_c\Delta u_c$$

(85)

where $G_c$ is the separation energy, $T_c$ is maximum traction before the damage and $\Delta u_c$ is the maximum displacement jump. $G_c$ and $\Delta u_c$ are the material parameters. When the maximum principal stresses (${\sigma }_p$) reaches $T_c$, the displacement jump ($\Delta u$) is evaluated, and the crack propagates when $\Delta u>\Delta u_c$.

The study of Jayendran and Ruimi [118] aims to investigate the state of stresses in the artery during crack propagation, which would help study the mechanics of aortic dissections. In their study, a three-layer artery model simulated a radial tear in the intima layer and a circumferential tear in the media layer. Wang et al. [119] have used a two-layer arterial model with residual stresses to study their effect on propagating arterial dissections. The initial tear was placed circumferentially in the middle of the media loaded with internal pressure. Both the studies have used the anisotropic hyperelastic model of HGO [34] for the arteries. Additionally, both the studies authors have defined crucial parameters for CZM from uniaxial tension results by Holzapfel [75].

4.2. Cohesive Zone Modelling

In XFEM with cohesive law, a traction separation law is defined along the surface where crack propagates, which is not known as apriori. While in CZM, the traction separation defines the onset of crack and its propagation along the predefined crack propagation path. In CZM, a layer of the surface is placed between two bulk materials where the crack propagates. This layer of the surface is modelled with special elements that vanish when the crack propagates based on the traction separation and evolution. In this section, CZMs of the soft tissue is discussed.

Badel et al. [120] applied CZM to an atherosclerotic coronary artery to study dissection mechanisms triggered due to angioplasty. The simulation uses a two-dimensional model of the artery, i.e., partially embedded in the myocardium and epicardium. And the artery is modelled with two layers of medial with a plaque incorporated. Medial layer material response is modelled with a Neo-Hookean model, and epicardium and myocardium are modelled as a linear elastic material. The damage initiation criteria ($DIC$) for the onset of the material degradation in the interface is defined as:

$$DIC=\max \left[\frac{{\delta }_n}{{\delta }_n^0},\frac{{\delta }_t}{{\delta }_t^0}\right]=1$$

(86)

where ${\delta }_n^0$ and ${\delta }_t^0$ are the maximum separation limits that define the damage initiation in normal and tangential directions, respectively.

The overall damage is characterized by $D,$ a scalar damage variable that is specified on the onset of damage in the interface. $D$ is a monotonically increasing variable from zero for without damage to one where the crack propagates. The effect of $D$ on the contact stress components is given as:

$${\sigma }_n=\left\{\begin{array}{c}\left(1-D\right)\overline{{\sigma }_n}, if \overline{{\sigma }_n}>0\\ \overline{{\sigma }_n}, if\overline{{\sigma }_n}<0\end{array}\right.$$

(87)

$${\sigma }_t=\left(1-D\right)\overline{{\sigma }_t}$$

(88)

where $\overline{{\sigma }_n}$ is normal and $\overline{{\sigma }_t}$ are tangential components of contact stress evaluated using elastic traction separation response without damage given as $\overline{{\sigma }_i}=Q_i{\delta }_i$. Here ${\delta }_i$ is the separation in the $i$ direction and $Q_i$ is stiffness parameter (units is MPa/mm) that specifies the separation and interfacial stress. The damage variable $D$ is defined as:

$$D={{\displaystyle \int }}_{{\delta }_m^0}^{{\delta }_m^f}\frac{{\sigma }_{n}d{\delta }_n+{\sigma }_{t}d{\delta }_t}{G_c-G_0}$$

(89)

where $G_c$ and $G_0$ are the critical fracture energy and elastic energy at damage initiation, respectively. And ${\delta }_m^0$ & ${\delta }_m^f$ are effective separation at the initiation of damage and propagation of the crack. Wherein the effective separation defined as ${\delta }_m=\sqrt{{\delta }_n^2+{\delta }_t^2}$. The critical parameters required for the traction separation are defined from the literature.

Leng et al. [18,121] have used CZM in two different studies applied to the artery. In one study, the fibrous cap delamination process is simulated to investigate its underlying process resulting in the delamination. Furthermore, the other study uses CZM to model the delamination in medial layers of the artery to study the failure mode. In both these studies, HGO [34] model is used to model the anisotropic hyperelastic response. Additionally, in Leng et al. [121] study, viscoelasticity is also considered. Effective displacement jump, $\delta ,$ and effective traction, $t$ are used to define the CZM, and they are given as:

$$\delta =\sqrt{{\lambda }^2\left({\delta }_{s1}^2+{\delta }_{s2}^2\right)+{\delta }_n^2}$$

(90)

$$t=\sqrt{{\lambda }^{-2}\left(t_{s1}^2+t_{s2}^2\right)+t_n^2}$$

(91)

while ${\delta }_{s1}$ & ${\delta }_{s2}$ represents the shearing and tearing displacements in tangential directions of the cohesive surface, ${\delta }_n$ is the opening displacement, and $\lambda$ is a scalar parameter that assigns weights to the displacements. Similarly, $t_n$, $t_{s1}$, and $t_{s2}$ are the tractions in normal and two shear directions across the cohesive surface.

In loading conditions, by using exponential CZM [106], the effective traction is given by

$$t=e{\sigma }_c\frac{\delta }{{\delta }_c}\exp \left(-\frac{\delta }{{\delta }_c}\right), if \delta \le {\delta }_{max} or \dot{\delta }\ge 0$$

(92)

$$t_n=K{\delta }_n, if {\delta }_n<0$$

(93)

Additionally, the effective traction during the unloading condition is given by:

$$t=\left(\frac{t_{max}}{{\delta }_{max}}\right)\delta , if \delta <{\delta }_{max} or \dot{\delta }<0$$

(94)

where $e=\exp \left(1\right)\approx 2.71828$, ${\sigma }_c$ is cohesive strength of the material, the maximum effective displacement is defined as ${\delta }_c=\frac{G_c}{e{\sigma }_c}$, $G_c$ is critical energy release rate, which is a material constant, $K$ is the penalty stiffness of the penetration resistance, ${\delta }_{max}$ is the maximum effective displacement during one delamination cycle, and $t_{max}$ is the corresponding effective traction.

A scalar damage parameter $d$ is defined on the onset of the damage, which monotonically increases from zero for without damage to one where the crack propagates. The damage parameter based on the displacement jump function is defined as:

$$d=1-\left(1+\frac{{\delta }_{max}}{{\delta }_c}\right)\exp \left(-\frac{{\delta }_{max}}{{\delta }_c}\right)$$

(95)

The crack initiates when $\delta ={\delta }_c$ resulting in the loss of the effective load-carrying capability of the cohesive elements. When $\delta ={\delta }_{sep}$, the material does not carry any load, and the crack propagates with the element deletion. Critical parameters required for their CZM are derived from the delamination experimental results.

Noble et al. [122] used CZM to simulate catheter-induced dissections (CID) in the porcine aorta. The tissue’s constitutive response was modelled with Ogden hyperelastic model given by:

$$\mathsf{\Psi }={\displaystyle \sum }_{p=1}^2\frac{{\mu }_p}{{\alpha }_p}\left({\overline{\lambda }}_1^{{\alpha }_p}+{\overline{\lambda }}_2^{{\alpha }_p}+{\overline{\lambda }}_3^{{\alpha }_p}-3\right)+\frac{9}{2}K{\left(J^{\frac{1}{3}}-1\right)}^2$$

(96)

where ${\mu }_p$ is shear modulus and ${\alpha }_p$ are the dimensionless constants such that $\mu =\frac{1}{2}{\displaystyle {\displaystyle \sum }_{p=1}^2}{\mu }_p{\alpha }_p$. ${\overline{\lambda }}_1, {\overline{\lambda }}_2, {\overline{\lambda }}_3$ represents the principal stretches.

The cohesive zone given by Bosch et al. [123] was employed to account for the large deformations. Since the opening and traction vectors are evaluated globally, no distinction was made between normal and tangential directions. The traction vector $t=te$ and separation vector $\delta =\delta e$ are related by:

$$t=\frac{G_c}{{\delta }_c} \left(\frac{\delta }{{\delta }_c}\right)\exp {\left[-\frac{\delta }{{\delta }_c}\right]}^{ }$$

(97)

where ${\delta }_c$ represents critical opening displacement, $G_c$ represents the critical energy release rate, and the unit vector $\mathit{e}$ is specified along the line bounded by all points opposing the interface.

Critical traction $t_{max}$ for the material at the critical opening point (i.e., $\delta ={\delta }_c$), $t_{max}$ is given by:

$$t_{max}=\frac{G_c}{\delta \exp \left(1\right)}$$

(98)

where $G_c$ is obtained from the experimental force-displacement graphs of the tissue dissection tests, whereas ${\delta }_c$ or ($t_{max}$) is found through the initiation zone of the same experimental force-displacement graphs [49]. Noble et al. [122] have conducted wedge dissection experiments for defining rupture parameters for the CZM. In their study, the crack was propagated by element deletion once the critical traction condition was achieved.

The CZM developed by Maiti and Geubelle [124] was applied to different soft tissue for tear studies. Fortunato et al. [125] applied this model to simulate tear propagation in arterial tissue during uniaxial testing. Ferrer et al. [126] used the model of Maiti and Geubelle to simulate tear propagation in tendons to study the effect of localized tendon remodelling. The traction separation law is normal ($t_n$) and tangential ($t_t$) direction across the cohesive surface is defined as:

$$t_n=\frac{d}{1-d}\frac{{\sigma }_{max}}{d_{ini}} \left(\frac{{\delta }_n}{{\delta }_{nc}}\right)$$

(99)

$$t_t=\frac{d}{1-d}\frac{{\tau }_{max}}{d_{ini}} \left(\frac{{\delta }_t}{{\delta }_{tc}}\right)$$

(100)

where ${\sigma }_{max}$ & ${\tau }_{max}$ are the maximum normal and shear stress of the cohesive zone, respectively, ${\delta }_{nc}$ & ${\delta }_{tc}$ are the critical displacement jumps corresponding to normal and tangential direction, ${\delta }_n$ & ${\delta }_t$ are the displacement jumps corresponding to the normal and tangential direction and $d$ is the monotonically decreasing damage variable which is scalar and is defined as:

$$d=\min \left[d_{ini}, 1-\delta \right]$$

(101)

here $d_{ini}$ is initially defined as damage parameter, Macauley bracket is defined as $1-\delta =\max \left(0, 1-\delta \right)$ and the magnitude of displacement jump is defined as $\delta :=\sqrt{{\delta }_n^2+{\delta }_t^2}$. Both these studies used the HGO model to define an anisotropic hyperelastic material response and CZM parameters from the literature. The study of Ferrer et al. [126] has used only traction in the normal direction since the tissue is under uniaxial tension.

4.3. Crack Phase-Field Modelling

The crack phase-field model (CPFM) defines the discontinuity with a special field equation along with the balance of linear momentum equation for the continuous field describing the elastic response of the material. A Ginzburg-Landau equation for defining the crack phase was introduced by Hakim and Karma [127]. Primarily two field variables are used in CPFM, i.e., the deformation map ($\mathit{\phi }$) and crack phase-field ($d$), as shown in Figure 8a,b, respectively. The crack phase-field ($d$) is solved with the crack evolution equation (102), while the deformation field ($\mathit{\phi }$) is solved with linear momentum balance (103).

The governing equations of CPFM problem to model fracture in anisotropic hyperelastic solid given as,

$$J div \left(J^{-1}\mathit{\tau }\right)+{\rho }_0\overline{\mathit{\gamma }}=0$$

(102)

$$\left(d-l^2\Delta d-\left(1-d\right)\overline{\mathcal{H}}\right)\dot{d}=0$$

(103)

where the Jacobian is defined with deformation gradient $\mathit{F}$ as $J:=\det \mathit{F}$, ${\rho }_0$ is the density, $\overline{\mathit{\gamma }}$ is the prescribed body forces, $l$ is the length scale parameter, and Kirchhoff stress tensor for the multi-field problem is defined as $\mathit{\tau }:=g\left(d\right){\mathit{\tau }}_0$. Where ${\mathit{\tau }}_0$ is the stress tensor for the rupture free material and $g\left(d\right)$ is the monotonically diminishing quadratic function given by:

$$g\left(d\right)={\left(1-d\right)}^2$$

(104)

with the boundary conditions that describe the evolution of the phase-field results in the degradation of the tissue as defined by:

$$g^{\prime }\left(d\right)\le 0$$

(105)

The degradation is ensured by the above condition, with $g\left(0\right)=1, g\left(1\right)=0$ acts as the limits for the flawless and torn state of the material, $g^{\prime }\left(1\right)=0$ illustrates a saturation as $d\to 1$.

The dimensionless crack driving function is given as:

$$\overline{\mathcal{H}}=\frac{{\mathsf{\Psi }}_0}{g_c/l}$$

(106)

Two distinct energy-based failure measures for ground matrix and fibres are considered. Accordingly, the isotropic and anisotropic strain energies of the dimensionless crack driving function is decomposed as:

$$\overline{\mathcal{H}}={\overline{\mathcal{H}}}^{iso}+{\overline{\mathcal{H}}}^{ani}$$

(107)

$${\overline{\mathcal{H}}}^{iso}=\frac{U_0+{\mathsf{\Psi }}_0^{iso}}{g_c^{iso}/l}, {\overline{\mathcal{H}}}^{ani}=\frac{{\mathsf{\Psi }}_0^{\mathrm{ani}}}{g_c^{ani}/l}$$

(108)

where $g_c^{iso}/l$ is critical fracture energy of the ground matrix per length scale, similarly $g_c^{ani}/l$ is for the fibres. They are defined from the experimental results of Sommers et al. [128]. For example, the anisotropic hyperelastic constitutive model for the intact artery is additively decomposed into the isotropic part with a neo-Hookean, and an exponential form considers the contribution of fibres for the anisotropic part [24].

$${\mathsf{\Psi }}_0=U_0+{\mathsf{\Psi }}_{dev}+{\mathsf{\Psi }}_{af1}+{\mathsf{\Psi }}_{af2}$$

(109)

where

$$U_0=\kappa \left(J-\ln J-1\right)$$

(110)

The deformation and crack phase-field are decoupled to subproblems with one-pass operator-splitting to solve the multi-field problem. The non-convex multi-field problem is divided into two convex sub-problems that are numerically simple to simulate when compared to the monolithic scheme. Gultekin et al. [34] have extended the approach by considering the fibre distribution to incorporate the anisotropy in the crack phase-field. Further, this model was applied to the artery peel test with different stress-based and energy-based criteria, and they showed that energy-based criteria are well suited for soft biological tissues [112].

5. Discussion

CDM was used to study all the damage phenomena, namely the Mullins effect, hysteresis, and permanent set. The CDM approach has been widely used to model damage phenomena across various tissues: artery, rectus sheath, ligament, annulus fibrosis, blood vessels, and ascending thoracic aortic aneurysm. Depending upon the tissue HGO model was used to define the collagen fibre structure, and appropriate strain energy was used to define the matrix material. Damage was introduced in the deviatoric part of the strain energy density in the form of Kachanov [129] based internal variable. Then, the equivalent strain concept of undamaged material by Simo and Ju [54] is used to define damage initiation. Damage propagation is phenomenologically modelled based on the behaviour of the tissue under supraphysiological loading. Based on the modelling of the targeted phenomenon, each reviewed model follows a specific approach.

To capture tear phenomena in soft tissues, a combination of CDM and SPH was used by Rausch et al. [67], and CDM combined with XFEM was proposed by Fathi et al. [69]. The model proposed by Ghasemi et al. [35] also captured soft tissue rupture. In addition, to circumvent the mesh dependency in CDM, a local gradient enhancement [66,74] and nonlocal integrals were applied [70,76]. While most of the reviewed models used a stretch-based damage evolution, Mousavi et al. [71] have used a fracture energy-based damage evolution. CDM models were developed in conjunction with various constitutive models, namely neo-Hookean, Ogden, and HGO. Additionally, Gao et al. [69] have developed a multi-physics model that considers water content in the constitutive model. While the damage model of Blanco et al. was developed from the mesoscopic level to model softening phenomena.

In pseudo-elasticity, the study of Pierce et al. [84] on the diseased and healthy aortic aneurysms using damage experiments have found a good agreement of their damage model. Wherein it demonstrates the capability of the pseudo-elasticity approach in modelling the softening and permanent. Further, the study of Holzapfel and Ogden [86] extends the pseudo-elasticity approach to consider the microstructural effect of collagen cross-links. Mainly, the effect of cross-links density and their interactions are considered, which can be physically interpreted as a soft tissue property. In addition to the softening and permanent set, Holzapfel and Ogden model [86] was able to capture the stiffening effect of the fibrous tissues with an increase in the density of the cross-links, while hyperelastic softening is a recently developed approach, and its application is limited to skin and arteries. The study of Li and Luo [63] has adopted the skin model and defined the parameters for swine, human, rabbit, and bovine skins. In a recent study, Volokh [88] has extended the initial model for two families of fibres, which can be applied to various soft tissue.

All damage modelling approaches are summarized in

Table 3. Summary of damage modelling applications.

Damage Approach	Capabilities	Tissues	Benefits
CDM	Mullins effect Hysteresis Permanent set	Artery, Rectus sheath, ligament, annulus fibrosus, ascending aortic aneurysm, thoracic aneurysm.	Rupture can be simulated. Applied to multiphysics models. Applied across different constitutive models. Mesh-size independent.
Pseudo-elasticity	Mullins effect Hysteresis Permanent set	Aortic aneurysms, brain tissue	Less number of material parameters. Parameters with a physical meaning.
Softening hyperelasticity	Permanent set	Skin, artery	Damage is incorporated in the constitutive model. A simple model does not involve any internal variables and damage evolution equations.

with their capabilities, application to tissues, and benefits. In general, a stretch based criteria defined by Simo [67] is used to define damage initiation and its evolution. Mesh dependency can be surpassed by using a local gradient enhancement or nonlocal integral in damage evolution. In addition, CDM can be used for multi-physics problems and microstructural based problems. However, the number of damage parameters and their evolution equations increases the computation. The hyperelastic constitutive response is inherent when the HGO model is used for soft tissues, and it further elevates by damage. The fundamental behaviour of the damage model in pseudo-elasticity makes it a good fit for continuous and discontinuous softening, and it is also extended for the permanent set. Even though the model is straightforward for implementing numerical simulations, its application is limited to arteries and brain tissues. Particularly, Holzapfel and Ogden [86] damage model phenomenologically captures the damage with an additional physical parameter, i.e., cross-links. Pseudo-elasticity damage in the material is controlled by the maximum strain attained, making it numerically simple. Lastly, in softening hyperelasticity, the damage is incorporated in the constitutive model and does not involve any damage variables and their evolution equations. This approach is applied to model the permanent set in the artery and skin. As the model is still evolving, its capabilities and limitations are not fully explored.

Table 4. Damage models applied to tissues and their validation method.

Reference	Tissue	Mechanism	Validation
Blanco et al. [64]	Soft tissue	Mullins effect	Numerical simulation of tropocollagen failure by Buehler et al. [132]
Comellas et al. [65]	Rectus sheath	Mullins effect	Based on the uniaxial tension experiments of Martins et al. [72]
Polindara et al. [66]	Blood Vessel	Permanent set	Wedge geometry simulation for balloon angioplasty was validated with analytical of neo-Hookean tube tests [75].
Ferreira et al. [67]	Arteries	Mullins effect	The damage model is not validated.
Rausch et al. [68]	Arteries	Permanent set Rupture	Damage model with results of Stepmer et al. [133] and tear simulations with Tong et al. [134] and Sommer et al. [128]
Fathi et al. [69]	Rectus sheath Ligament	Mullins effect Rupture	Uniaxial tension experiments of Martins et al., [72] and numerical results of Waffenschmidt et al. [73] for rectus sheath. For ligament, the model is validated with experimental results of Weiss [77] and numerical results of Calvo et al. [135]
Gao et al. [72]	Annulus fibrosus	Permanent set	Simulation results validated with the experimental results of Ebara et al. [136] and Skaggs et al. [137]
Mousavi et al. [71]	Ascending thoracic aortic aneurysms	Permanent set	Buldge inflation test with graft size of 45 × 45 mm2 and inflation of circular area of diameter 30 mm.
Ghasemi et al. [35]	Arteries	Mullins effect Hysteresis Permanent set	Experiments of uniaxial tension tests and cyclic loading in uniaxial tension.
Pierce et al. [84]	Thoracic aortic tissues Abdominal aortic tissues	Permanent set Mullins effect	Experimental results from tissues under uniaxial tension and cyclic loading.
Holzapfel and Ogden [86]	Soft tissue	Mullins effect	A reduced model with uniaxial fibres is validated with the rat tail tendon experiment results of Pins and Silver.
Li and Luo [63]	Skin	Permanent set	Experimental results of Annaidh et al. [138] for human skin under uniaxial tension. Additionally, validate with porcine skin.
Volokh [88]	Artery adventitia	Permanent set	Uniaxial tension tests of artery adventitia in longitudinal and circumferential directions.

summarises the reviewed damage models, mechanisms and their validation methods. The authors of the reviewed models have validated with basic mechanical tests, which may differ from the physiological condition. For instance, damage models for rectus sheath are validated with the uniaxial tests by Martins et al. [72]. However, the biaxial tension would better represent the physiological loading condition in the rectus sheath. The same was mentioned as a limitation of the study by Martins et al. [72].Conventional mechanical tests can give the tissue essential mechanical properties. However, they may not behave similarly for the intended application [130]. Damage models developed for arterial tissues can be validated with internal pressure tests by Perez et al. [131]. The results are applicable in the domains such as balloon angioplasty and aneurysms. A damage model validated with physiological loading conditions would enable clinically translatable simulations.

For rupture modelling, most of the studies focused on the arterial tissues, while few were focused on the tendon. In XFEM, except for the study of Karimi et al. [116], rupture is simulated using a two-dimensional model with plane strain approximation. This approximation would be computationally efficient but oversimplifies the problem. While in the study of Karimi et al. [116], the virtual crack method of XFEM is used to simulate crack propagation in arteries. The CZM has been extensively used to study the delamination in arterial tissues. A good experimental agreement was reported in the reviewed studies for arterial delamination. Particularly, the CZM by Maiti and Geubelle [124] was used to study arterial rupture under uniaxial tension by Fortunato et al. [125], and the same model was applied for studying rupture in tendons by Ferrer et al. [126]. While XFEM is computationally expensive, CZM requires the crack path apriori. CPFM overcomes both the limitations of XFEM and CZM by dealing with rupture as a multi-field problem. A comparison of the three rupture approaches was reported by Gultekin et al. [33]. In CPFM, various stress-based and energy-based criteria can be used for crack initiation and propagation. However, as CPFM was recently adopted for soft tissues, its application is limited to arteries. Since CPFM uses multiple families of fibres and is given the freedom of using various crack initiation and propagation criteria, its implementation can be extended to other soft fibrous tissues.

The modelling of the damage parameter considers representing the damage in an inactive tissue by neglecting all the biological aspects of tissue [6,32]. In the discussed damage models, damage initiation is considered based upon the loading condition, in which the damage is initiated after reaching a particular load or particular stretch in the tissue caused by the load. Such models work to simulate the soft tissues under supra-physiological loading, i.e., when the tissues undergo loading higher than the physiological limit. In general, supra-physiological loading occurs due to external loading, for instance, in anterior cruciate ligament tear, catheter induced dissections, balloon angioplasty, etc. The damage model for diseased tissue requires knowledge of supra-physiological loads as well as pathological effects. However, these effects are not considered while modelling damage in diseased tissue.

Therefore, numerical simulation of the damage due to pathological conditions demands a model where its parameters represent the mechanical and physiological changes due to disease [6,139]. Such a model would be able to describe both the anatomical and physiological changes. Developing the aforementioned damage model requires constitutive model parameters with physical interpretation [32] and robust experiments in conjunction with tissue engineering to study the effect of disease on the tissue [139]. However, the active response of the soft tissues in damage and rupture models can be developed by introducing the growth and remodelling [140,141,142]. The damage models by Ghasemi et al. [36] considers damage mechanisms at mesoscopic scales, and Holzapfel and Ogden [86] consider the cross-links between the collagen fibres, where both the models aim towards defining the constitutive parameters with physical interpretation. So far, the studies reported are based on the experiments of healthy tissue and diseased tissue, which gives the constitutive response of the tissue in different conditions. However, the effect of the disease progression on the constitutive response is still at large. The advancements in tissue engineering and evolving in-vitro disease modelling [143] can enable the experiments to study the effect of disease progression [144] in the tissues.

6. Conclusions

A series of state-of-the-art damage and rupture models for modelling soft tissue failure were reviewed. The damage models were classified based on the approach of employing the damage in the soft tissues. Similarly, rupture models are grouped based on the method to deal with discontinuity during the rupture process. However, the present study has two limitations; firstly, it does not cover the modelling aspects of plastic phenomena (residual strains) and the capabilities of the discussed models to capture these phenomena. Secondly, the nanomechanics of the soft tissues to model damage and rupture were not reviewed.

In damage modelling, CDM and pseudoelasticity are widely applied to various tissues to model the damage under mechanical loading. Even though few microstructure-based damage models have been developed, there is a definite need for models that can consider the physiological effects of disease progression simulations. CDM based modelling uses reduction factor based on Kachanov [129] and equivalent strain concept of Simo and Ju [54], where it requires an initiation condition, evolution function for damage initiation and progression, respectively. However, in the pseudoelasticity, the damage variables are defined based on the maximum strain that occurred in the loading along with some damage control parameters. Lastly, the softening hyperelasticity approach uses an energy limiter incorporated in the strain energy density to model the softening effect. Further, it neither needs internal variables nor threshold conditions.

In rupture modelling, along with classical XFEM and CZM, recently developed CPFM is reviewed. CPFM overcomes complications of computations associated with XFEM and path dependency associated with CZM. However, the application of CPFM is limited to the artery, and its application to other tissues needs to be explored. The damage and rupture models reviewed here do not consider viscoelasticity, fibre recruitment, etc. These are considered as material effects and physiological effects. However, the reviewed damage models demonstrated their capability to capture different damage phenomena, and they were applied to various tissues in both humans and animals under supra-physiological loading. These damage and rupture models need to be extended to other tissues that can find biomedical and clinical research applications. For instance, simulation of the sutured condition of the vascular or skin grafts to evaluate the near in-vivo fracture toughness would help design artificial grafts. Further, the damage model of Gao et al. [70] can be extended to simulate ocular infections as this model captures the softening and tissue hydration. Moreover, corneal ectasia, where localised progressive softening occurs, can be studied by extending the model of Volokh to the cornea [88].