1 Introduction

Soft tissue is a general term that refers to various groups of cells in the human body. All tissues in the body that are neither bones nor organs are considered soft tissues. Soft tissues can be divided into connective tissues, such as tendons, ligaments, fascia, skin, fibrous tissues, fat, and synovial membranes, and non-connective tissues, such as muscles, nerves, and blood vessels [1]. The major physiological functions of soft tissues are to connect, support, and surround organs and other structures of the body. Some soft tissues, namely arterial smooth muscle, skeletal muscle and cardiac muscle, are stretched to generate a passive tension, but can also contract to generate an active tension. However, other soft tissues such as ligaments, tendons, and skins can be stretched to have passive tension only.

Soft tissues, especially arterial smooth muscles, ligaments, and tendons in joints of the human body can be injured or damaged by disease or excessive force applied during exercise, accidents, or surgery. For example, the human arterial inner wall can be damaged by high blood pressure, with subsequent plaque development. A tendon can be damaged or ruptured, as shown in Fig. 1. Knowledge of soft tissue damage, injury, or failure behavior is very helpful for artificial soft tissue design and fabrication.

Fig. 1
figure 1

Artery inner wall is damaged by high blood pressure of heart and achilles tendon is ruptured by accident. a artery (http://www.webmd.com/hypertension-high-blood-pressure/how-high-blood-pressure-damages-arteries) and b tendon (http://www.methodistorthopedics.com/achilles-tendon-problems)

Full size image

In addition to experiments, finite element analysis (FEA) is extensively used for soft tissue damage characterization when the tissue is in vivo state for surgery [2]. However, the damage models are based on existing in vitro measurements. Unfortunately, soft tissues usually exhibit nonlinear, heterogeneous, anisotropic, and viscoelastic behavior, making their damage modeling difficult. Thus, the understanding, physical description, and modeling of damage and failure in soft tissues have presented a tough challenge. It is necessary to assess existing soft tissue damage models and to take a forward look for further study.

Damage models developed for soft tissues can be traced back to the 1970s. These models can be divided into three categories: (1) deterministic models, in which a pseudo-elastic strain energy function with a few parameters or a strain energy function with a few damage variables of continuum damage mechanics is used to account for the softening/damage effect, and soft tissue can be either isotropic or anisotropic; (2) probabilistic models, in which the fiber recruitment effect, probabilistic damage process, or both are involved; in some models a damage variable of continuum damage mechanics is used; (3) microstructure-based damage models of collagen fibers, in which the individual collagen fibril damage behavior is characterized and then integrated into the whole tissue level. The last two kinds of model are anisotropic only. These damage models are summarised in Table 1.

Table 1 Summary of damage models for soft tissues
Full size table

The rest of this paper is organized as follows. Section 2 reviews damage models for isotropic soft tissues as they provide a basis for damage modeling, especially the model for rubber-like materials proposed by Ogden and Roxburgh [3]. In Sect. 3, a survey of damage models for anisotropic soft tissues, including the arterial wall, ligaments, and tendons, is given. Discussions and a few challenges for damage models are highlighted in Sect. 4. Finally, concluding remarks are given in Sect. 5.

2 Damage Models for Isotropic Materials

2.1 Rubber-Like Materials

The concept of continuum damage mechanics has been increasingly applied to predict damage in soft tissues. Damage mechanics involves the engineering predictions of the initiation, propagation, and fracture in a material using state variables, which represent the effects of damage caused from thermal or mechanical loading or aging on the stiffness and remaining life of the material [4, 5]. The state variables may be measurable variables or other physical variables.

The damage of soft isotropic materials is modeled in continuum damage mechanics as follows. A specimen of rubber or soft tissue often exhibits the Mullins effect in a simple tensile test under a cyclic load, as sketched in Fig. 2. Initially, the sample is loaded to b′ from a; the loading path is a b b′. If the sample is unloaded from b′, then the unloading path is b′ B a. If the sample is reloaded to point c′, then the loading path will be a B b′ c c′ and the new unloading path will be c′ C a. Such a path pattern is repeated until a total failure occurs in the sample under cyclic loading. This phenomenon is named the Mullins effect or stress softening.

Fig. 2
figure 2

Sketch of loading–unloading paths exhibiting the Mullins effect in simple tensile test under cyclic loading, adapted from [3]. Note that because unloading curve comes back to the origin, there is no plastic deformation in material

Full size image

The Mullins effect is interpreted as damage occurring in the sample at the microscopic level. This may be due to the bonds between the filler particles and the molecular chains being broken for rubber or collagen fiber being broken for soft tissues.

Since the lengths of chain links in rubber and collagen fibers vary, and they can break at different stretches as the damage process proceeds. Additionally, after damage, the reloading path is identical to the previous unloading path. This suggests that the energy consumed by damage is irrecoverable.

For a general biaxial simple tensile test, the following specific pseudo-elastic strain energy function was proposed by Ogden and Roxburgh [3] for an incompressible rubber-like material:

$$ \psi \left( {\lambda_{1},\lambda_{2},\eta } \right) = \eta \tilde{\psi }\left( {\lambda_{1},\lambda_{2} } \right) + \phi \left( \eta \right) $$
(1)

where λ 1 and λ 2 are two principal stretches, \( \tilde{\psi }\left( {\lambda_{1},\lambda_{2} } \right) \) is an Ogden-type strain energy without damage, expressed as:

$$ \tilde{\psi }\left( {\lambda_{1},\lambda_{2} } \right) = \mu \sum\limits_{i = 1}^{3} {\frac{{\mu_{i} }}{{\alpha_{i} }}\left( {\lambda_{1}^{{\alpha_{i} }} + \lambda_{2}^{{\alpha_{i} }} + \lambda_{1}^{{ - \alpha_{i} }} \lambda_{2}^{{ - \alpha_{i} }} - 3} \right)} $$
(2)

where α 1, α 2, α 3, μ 1, μ 2, and μ 3 are the material constants without damage and μ is the shearing modulus. \( \phi \left( \eta \right) \) is the damage function that satisfies the following equation:

$$ - \frac{d\phi }{d\eta } = m{\mathbf{erf}}^{ - 1} \left[ {r\left( {\eta - 1} \right)} \right] + \tilde{\psi }\left( {\lambda_{1m},\lambda_{2m} } \right) $$
(3)

where λ 1m and λ 2m are the principal stretches for the point where unloading has most recently been initiated from the primary loading path; erf −1 () is the inverse of the error function; m and r are positive material constants, where m is a parameter used to control the dependence of the damage on the extent of deformation, and r is a variable used to indicate the extent of damage relative to the virgin state. The variable η is expressed as:

$$ \eta = 1 - \frac{1}{r}{\mathbf{erf}}\left[ {\frac{1}{m}\left( {\tilde{\psi }\left( {\lambda_{1m},\lambda_{2m} } \right) - \tilde{\psi }\left( {\lambda_{1},\lambda_{2} } \right)} \right)} \right] $$
(4)

where η ∈ [0,1] and the error function is \( {\mathbf{erf}}\left( x \right) = \tfrac{2}{\sqrt \pi }\int\limits_{0}^{x} {e^{{ - t^{2} }} dt} . \) The boundary condition \( \phi \left( 1 \right) = 0 \) holds for Eq. (3). After damage, the principal stresses are estimated as:

$$ \sigma_{i} - \sigma_{3} = \eta \lambda_{i} \frac{{\partial \tilde{\psi }}}{{\partial \lambda_{i} }}, \quad i = 1,2 $$
(5)

For the simple tensile case, λ = λ 1, λ 2 = λ 3 = λ −1/2, and σ 2 = σ 3 = 0. The parameters μ, r, and m can be determined from a series of experimental stress–stretch curves. This model has been extended to the case with residual strain [6]. For industrial rubber, the Ogden model seems better than another damage model [7].

In [8], a general continuum damage mechanics model, the Ogden pseudo-elastic model just mentioned above, and Guo’s elastic model were compared against a series of rubber-like material experimental stress–strain curves under cyclic loads to evaluate model performance. It was shown that Ogden’s pseudo-elastic and Guo’s elastic models are better than the general continuum damage mechanics model.

2.2 Continuum Damage Mechanics Method for Isotropic Materials

In this section, the basic idea for modeling the damage effect in continuum damage mechanics is introduced. Considering an isotropic, homogenous, incompressible rubber-like material, the strain energy function in the damage state, \( \psi \left( {{\mathbf{F}},d} \right), \) can be expressed in terms of a strain energy function without damage, \( \tilde{\psi }\left( {\mathbf{F}} \right), \) and a scalar damage variable, d, as follows [8]:

$$ \psi \left( {{\mathbf{F}},d} \right) = \left( {1 - d} \right)\tilde{\psi }\left( {\mathbf{F}} \right) $$
(6)

where F is the deformation gradient tensor associated with the stress-free state, d is a continuous variable used to characterize the damage effect in the material, where d ∈ [0,1], d = 0, no damage; d = 1, complete damage; 0 < d < 1, a damage state in between. The second Piola–Kirchhoff stress with damage is calculated as:

$$ {\text{S}} = \left( {1 - d} \right)\frac{{\partial \tilde{\psi }\left( {\mathbf{F}} \right)}}{{\partial {\text{E}}}} = \left( {1 - d} \right){\tilde{\text{S}}} $$
(7)

where E is the Green–Lagrange strain tensor and \( {\tilde{\text{S}}} \) is the second Piola–Kirchhoff stress without damage. Usually, the scalar damage variable is an exponential function of the effective strain energy function α(t) at which damage occurs, i.e.:

$$ d\left( {\alpha \left( t \right)} \right) = d_{\infty } \left[ {1 - \exp \left( { - \frac{\alpha \left( t \right)}{\beta }} \right)} \right] $$
(8)

where d ∞ and β are the model parameters determined in experiments. α(t) can be determined as follows. Because damage is an irreversible process, the second law of thermodynamics should be applicable. In continuum mechanics, the second law of thermodynamics was expressed as the Clausius–Duhem inequality by Guo and Sluys [8]:

$$ - \frac{{d\psi \left( {{\mathbf{F}},d} \right)}}{dt} + {\text{S:}}\frac{{d{\text{E}}}}{dt} \ge 0 $$
(9)

Note that:

$$ \frac{{d\psi \left( {{\mathbf{F}},d} \right)}}{dt} = \frac{{\partial \psi \left( {{\mathbf{F}},d} \right)}}{{\partial {\text{E}}}}\frac{{d{\text{E}}}}{dt} + \frac{{\partial \psi \left( {{\mathbf{F}},d} \right)}}{\partial d}\dot{d} $$
(10)

According to Eq. (6), \( {{\partial \psi \left( {{\mathbf{F}},d} \right)} \mathord{\left/ {\vphantom {{\partial \psi \left( {{\mathbf{F}},d} \right)} {\partial d}}} \right. \kern-0pt} {\partial d}} \) can be written as:

$$ - \frac{{\partial \left( {{\mathbf{F}},d} \right)}}{\partial d} = \tilde{\psi }\left( {\mathbf{F}} \right) $$
(11)

Eventually, Eq. (9) is reduced to the very simple form:

$$ \tilde{\psi }\left( {\mathbf{F}} \right)\dot{d} \ge 0 $$
(12)

This suggests that the damage process is driven by the strain energy function (thermodynamic force) without damage, \( \tilde{\psi }({\mathbf{F}}). \) Accordingly, we can establish a damage criterion based on \( \tilde{\psi }\left( {\mathbf{F}} \right) \) from experiments under a series of loads versus time, i.e. \( \alpha \left( t \right) = {\text{Max }}\tilde{\psi }\left( {{\mathbf{F}}\left( t \right)} \right), \) which is the maximum strain energy function without damage. When \( \tilde{\psi }\left( {\mathbf{F}} \right) \) caused by a certain load profile equals α(t), damage will occur. This criterion is expressed mathematically as:

$$ \phi = \tilde{\psi }\left( {{\mathbf{F}}\left( t \right)} \right) - \alpha \left( t \right) \le 0 $$
(13)

If \( \phi < 0, \) there is no damage at all; otherwise, if \( \phi = 0, \) damage occurs. The variable α(t) can be used for damage evaluation with time. In this model, the parameters \( \tilde{\psi }\left( {\mathbf{F}} \right), \) α(t), d ∞, and β need to be determined from a series of experimental stress–stretch curves in various time-dependent loading courses.

Discontinuous and continuous damage evolution models were proposed by Miehe [7] for rubber-like materials based on the Ogden-type strain energy function. The models can deal with two damage problems: damage characterized by a function of maximum strain attained in a loading path and damage that is strain-rate-dependent (i.e., the viscoelastic effect).

In the former, there is no damage accumulation, and the maximum strain energy function without damage,\( {\text{Max }}\tilde{\psi }\left( {{\mathbf{F}}\left( t \right)} \right), \) during a loading path can serve as the failure criterion. In the latter, however, the damage accumulation exists during a cyclic loading path, and the maximum effective strain energy function, \( \int_{0}^{t} {\left| {{{d\tilde{\psi }\left( {{\mathbf{F}}\left( s \right)} \right)} \mathord{\left/ {\vphantom {{d\tilde{\psi }\left( {{\mathbf{F}}\left( s \right)} \right)} {ds}}} \right. \kern-0pt} {ds}}} \right|} ds, \) should be used as the damage criterion.

These two damage mechanisms [9, 10] were combined by Miehe [7] with two sets of damage variables in the continuum damage mechanics method shown above. Since the model includes the viscoelastic effect, which is beyond the scope of this review, it is not further discussed here.

2.3 Damage Model for Abdominal Aortic Aneurysms

A damage model for abdominal aortic aneurysms (AAAs) was proposed [11] based on the above continuum damage mechanics method. AAAs are considered as a compressible, homogenous, and isotropic material without any fibers. The strain energy function with damage is written as:

$$ \psi \left( {\mathbf{C}} \right) = \tilde{\psi }_{vol} \left( J \right) + \left( {1 - d} \right)\tilde{\psi }\left( {{\bar{\mathbf{C}}}} \right) $$
(14)

where \( {\mathbf{C}} = {\mathbf{F}}^{\mathbf{T}} {\mathbf{F}},\; J = \det \left( {\mathbf{F}} \right), \; \overline{\mathbf{C}} = J^{{ - {2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-0pt} 3}}} {\mathbf{C}}, \; {\mathbf{F}} = \partial {\mathbf{x}} /\partial {\mathbf{X}},\;{\mathbf{x}} \) is the current configuration, X represents the reference configuration, and d is the damage variable, which can be calculated as:

$$ d\left( {\alpha \left( t \right)} \right) = a\left[ {1 - \exp \left( { - b\alpha \left( t \right)} \right)} \right] $$
(15)

where \( \alpha \left( t \right) = \sqrt {2\tilde{\psi }\left( {{\mathbf{C}}\left( t \right)} \right)}, \) and a and b are model constants. This equation is the same as Eq. (8). Note that AAA tissues demonstrate anisotropic biomechanical properties [12–14], obviously, this isotropic damage may not be justified.

The energy limiter method is another alternative for dealing with the damage effect in isotropic materials such as rubber or rubber-like materials [15]. For solids, the energy limiter for damage/failure/rupture is equivalent to the bond energy, which can be measured using the strain energy function. The following constitutive law for AAA was proposed by Volokh [16] based on an isotropic strain energy plus an energy limiter:

$$ \psi \left( {\mathbf{C}} \right) = \phi - \phi \exp \left\{ { - {{\left[ {c_{1} \left( {tr{\mathbf{C}} - 3} \right) + c_{2} \left( {tr{\mathbf{C}} - 3} \right)^{2} } \right]} \mathord{\left/ {\vphantom {{\left[ {c_{1} \left( {tr{\mathbf{C}} - 3} \right) + c_{2} \left( {tr{\mathbf{C}} - 3} \right)^{2} } \right]} \phi }} \right. \kern-0pt} \phi }} \right\} $$
(16)

where \( \phi \) is the energy limiter, and c 1 and c 2 are material constants. These three model parameters need to be determined using uniaxial tension test data of AAAs. Note that there is no need to involve a material damage variable in this damage model, as described in Sect. 3.1.

3 Damage Models for Anisotropic Soft Tissues

3.1 Deterministic Damage Models

The artery wall is incompressible, nonlinear, and inhomogeneous, and exhibits hysteresis under a cyclic load. Its fibrous structure (i.e., collagen and elastin fibers) can be torn under a pressure higher than the physiological pressure. This micro-tearing is strain-related and contributes to the amount of damage. Like a rubber material, damage to the arterial wall is closely related to the maximum strain.

Under a steady load, the artery cannot be damaged until the maximum strain is achieved. Under a cyclic load, the loading and unloading stress–strain paths will remain unchanged until a previous maximum strain is exceeded. Such behavior is referred to as the Mullins or softening effect, which was originally used to describe rubber behavior.

In traditional methods, once the maximum von Mises stress or strain at a point in a material is beyond a criterion, the material is said to experience failure. However, such local failure does not lead to total failure in an artery. This means that traditional methods for predicting the failure of rubber materials may be unsuitable for arteries, and thus new methods are needed for predicting total failure in the arterial wall.

Arterial tissues are subject to the softening effect and have an S-shaped stress–stretch curve [17–20]. In [21], an anisotropic damage model was proposed to account for tensile and compressive damage based on continuum damage mechanics. This model has been applied to artery damage prediction [22]. However, fibers and the softening effect were excluded. It seems to be difficult to extend this model to tissues with fibers. Hence, this model is not further discussed here.

In the following paragraphs, we summarize a few kinds of fiber-based damage model. The important one is the artery biomechanical model proposed by Holzapfel et al. [23] updated with the damage effect.

The constitutive law is expressed by the strain energy function put forward by Holzapfel et al. [23] for arterial walls with an incompressible, homogenous matrix and two families of collagen fibers:

$$ \psi \left( {I_{1},I_{4},I_{6} } \right) = \frac{\mu }{2}\left( {I_{1} - 3} \right) + \frac{{k_{1} }}{{2k_{2} }}\sum\limits_{i = 4,6} {\left\{ {\exp \left[ {k_{2} \left( {I_{i} - 1} \right)^{2} } \right] - 1} \right\}} $$
(17)

where I 1, I 4, and I 6 are the stretch invariants, and I i  = f i •(Cf i ), i = 4, 6, where f i is the orientation vector of each family of fibers.

To accommodate the softening effect in the matrix material, the neo-Hookean model of the first term in Eq. (17) is updated as [24]:

$$ \psi \left( {I_{1},\phi } \right) = \phi - \phi \exp \left\{ { - \frac{\mu }{2\phi }\left( {I_{1} - 3} \right)} \right\} $$
(18)

The last two terms are modified as:

$$ \psi \left( {I_{4},I_{6},\xi_{4},\xi_{6},n_{4},n_{6} } \right) = \frac{{k_{1} }}{{2k_{2} }}\sum\limits_{i = 4,6} {\left\{ {\exp \left[ {k_{2} \left( {I_{i} - 1} \right)^{2} } \right] - 1 - \frac{{k_{2} }}{{2n_{i} + 1}}\left( {\frac{{I_{i} - 1}}{{\xi_{i}^{2} - 1}}} \right)^{{2n_{i} + 1}} } \right\}} $$
(19)

where \( \phi,\) \( \xi_{i}, \) and \( n_{i} \) are the damage parameters determined from experiments; the parameters \( \mu, \) \( k_{1}, \) and \( k_{2} \) are model constants without damage. This model is applicable for each layer of arterial walls.

In [25], the damage model described by Eqs. (18) and (19) was modified using a strain energy limiter, sharpness factor, and upper incomplete gamma function as follows:

$$ \psi \left( {I_{1},\phi } \right) = \frac{\phi }{m}\left\{ {\Upgamma \left( {\frac{1}{m},0} \right) - \Upgamma \left( {\frac{1}{m},\left[ {\frac{{\frac{\mu }{2}\left( {I_{1} - 3} \right)}}{\phi }} \right]^{m} } \right)} \right\} $$
(20)

and

$$ \psi \left( {I_{4},I_{6},\phi_{4},\phi_{6},m_{4},m_{6} } \right) = \sum\limits_{i = 4,6} {\frac{{\phi_{i} }}{{m_{i} }}\left\{ {\Upgamma \left( {\frac{1}{{m_{i} }},0} \right) - \Upgamma \left( {\frac{1}{{m_{i} }},\left[ {\frac{{\frac{{k_{1} }}{{2k_{2} }}\left( {e^{{k_{2} \left( {I_{i} - 1} \right)^{2} }} - 1} \right)}}{{\phi_{i} }}} \right]^{{m_{i} }} } \right)} \right\}} $$
(21)

where the parameters \( \phi,\) \( \phi_{4}, \) and \( \phi_{6} \) are the strain energy limiters for the matrix and two families of fibers, respectively, m, m 4, and m 6 are the sharpness factors for the matrix and two families of fibers, respectively, which can be determined from uniaxial or biaxial tensile test data, and Γ represents the upper incomplete gamma function, defined as \( \Upgamma \left( {s,x} \right) = \int\limits_{x}^{\infty } {t^{s - 1} e^{ - t} dt} . \)

Similarly, a damage model for soft tissues with fibers was proposed by Ehret and Itskov [26] to account for the softening effect. The generalized poly-convex strain energy function for an isotropic matrix and anisotropic fibers to meet the material stability criteria was adopted. The initial stiffness of fibers was reduced gradually to fit the experimental stress–stretch curves under cyclic loadings, suggesting that damage occurs only in the fibers.

Based on the continuum damage mechanics model described in Sect. 2.2, a damage model for artery walls was proposed by Peña et al. [27] and Balzani et al. [28]. It is considered that damage occurs in two families of fibers only. As a result, the strain energy function of the matrix material does not need to be modified. Only the last two energy functions for the fibers are updated to the following form:

$$ \psi \left( {I_{4} ,I_{6} ,d_{4} ,d_{6} } \right) = \frac{{k_{1} }}{{2k_{2} }}\sum\limits_{i = 4,6} {\left\{ {\exp \left[ {k_{2} \left\langle {\left( {1 - d_{i} } \right)\left( {\kappa I_{1} + \left( {1 - 3\kappa } \right)I_{i} } \right) - 1} \right\rangle^{2} } \right] - 1} \right\}} $$
(22)

where d i is the damage variable of fibers and κ is a constant associated with fiber orientation dispersion. Details can be found elsewhere [29]. This treatment of the damage effect is slightly different from that in Eq. (19).

Further, the following type of damage variable was applied to represent the damage process for the case where the maximum loading is fixed in a cyclic tension test:

$$ d_{i} \left( \beta \right) = d_{s} \left[ {1 - \exp \left( {\frac{{\ln \left( {1 - r_{s} } \right)}}{{\beta_{s} }}\beta } \right)} \right] $$
(23)

where \( \beta = \int_{0}^{t} {d\tilde{\psi }}, \) and r s and β s are model parameters, where β s is the variable at r s  = 0.99, \( r_{s} \in \left[ {0,1} \right] . \) The maximum damage variable d s is determined using:

$$ d_{s} = d_{\infty } \left[ {1 - \exp \left( {\frac{{\ln \left( {1 - r_{\infty } } \right)}}{{\alpha_{\infty } }}\alpha } \right)} \right] $$
(24)

where d ∞ denotes a predefined convergence limit for the overall damage value, \( d_{\infty } \in \left[ {0,1} \right] , \) and α ∞ is the variable at r ∞ = 0.99, \( r_{\infty } \in \left[ {0,1} \right] . \) The maximum strain energy function is \( \alpha \left( t \right) = {\text{Max }}\tilde{\psi }\left( {{\mathbf{F}}\left( t \right)} \right) . \) These two damage variables can be also found in [7].

In [27], based on discontinuous and continuous damage models [7], a damage model for the pig aorta under cyclic loads was proposed. The damage occurs in both the matrix and fibers. The damage variables are \( d_{i} = d_{i} \left( {\alpha_{i} \left( t \right)} \right) + d_{i} \left( {\beta_{i} \left( t \right)} \right) , \) i = m, 4, 6 for the matrix and two families of fibers, respectively, where \( \alpha_{i} \left( t \right) = {\text{Max}}\sqrt {2\tilde{\psi }_{i} \left( {{\mathbf{F}}\left( t \right)} \right)} \) and \( \beta_{i} \left( t \right) = \int_{0}^{t} {\left| {{{d\tilde{\psi }_{i} \left( {{\mathbf{F}}\left( s \right)} \right)} \mathord{\left/ {\vphantom {{d\tilde{\psi }_{i} \left( {{\mathbf{F}}\left( s \right)} \right)} {ds}}} \right. \kern-0pt} {ds}}} \right|} ds . \) The strain energy function with damage for the matrix and two families of fibers is written as:

$$ \psi = \psi_{vol} \left( J \right) + \left( {1 - d_{m} } \right)\tilde{\psi }_{m} + \left( {1 - d_{4} } \right)\tilde{\psi }_{4} \left( {I_{4} } \right) + \left( {1 - d_{6} } \right)\tilde{\psi }_{6} \left( {I_{6} } \right) $$
(25)

where \( \tilde{\psi }_{m} , \) \( \tilde{\psi }_{4} , \) and \( \tilde{\psi }_{6} \) are the strain energy functions of the matrix and two families of fibers without damage, respectively. The damage criterion for the discontinuous damage process is:

$$ \phi_{i} \left( {{\mathbf{F}}\left( i \right),\zeta_{i} } \right) = \sqrt {2\tilde{\psi }_{i} \left( {{\mathbf{F}}\left( t \right)} \right)} - \alpha_{i} \left( t \right) = \zeta_{i} \left( t \right) - \alpha_{i} \left( t \right) \le 0 $$
(26)

The damage variable evolution equation updated by Peña et al. [27] is expressed as:

$$d_{i} \left( {\zeta_{i} } \right) = \left\{ {\begin{array}{lll} 0 & \quad \zeta_{\text{i}} < \zeta_{i}^{\hbox{min}} \\ \frac{1}{2}\left[ 1+ \frac{2\xi_{i} \Uplambda_{i} \exp \left( 2\xi_{i} \left[2\Uplambda_{i} - 1 \right] \right) - 1}{2\xi_{i} \Uplambda_{i} \exp\left( 2\xi_{i} \left[ 2\Uplambda_{i} \right] \right) + 1} \right] & \quad \zeta_{\text{i}} \in \left[\zeta_{i}^{\hbox{min}},\zeta_{i}^{\hbox{max}} \right]\\ 1 & \quad \zeta_{\text{i}} >\zeta_{i}^{\hbox{max} } \end{array} } \right. $$
(27)

where the variable \( \Uplambda_{i} = {{\left( {\zeta_{i} - \zeta_{i}^{\hbox{min} } } \right)} \mathord{\left/ {\vphantom {{\left( {\zeta_{i} - \zeta_{i}^{\hbox{min} } } \right)} {\left( {\zeta_{i}^{\hbox{max} } - \zeta_{i}^{\hbox{min} } } \right)}}} \right. \kern-0pt} {\left( {\zeta_{i}^{\hbox{max} } - \zeta_{i}^{\hbox{min} } } \right)}} , \) and \( \zeta_{i}^{\hbox{min} } \) and \( \zeta_{i}^{\hbox{max} } \) are the strain energy function values for damage occurrence and total failure, respectively. The model parameter \( \xi_{i} \)  > 0. The continuous damage variable \( d_{i} \left( \beta \right) \) has the following form:

$$ d_{i} \left( \beta \right) = d_{i\infty } \left[ {1 - \exp \left( { - \frac{\beta }{{\gamma_{i} }}} \right)} \right] $$
(28)

where \( d_{i\infty } \) is the maximum possible continuous damage in the matrix and fibers, \( d_{i\infty } \in \left[ {0,1} \right] , \) and \( \gamma_{i} \) is the damage saturation parameter. This model requires ten experimental damage parameters and five constitutive law constants. Although the discontinuous and continuous damage models have good agreement with experiments, the plastic effect is not presented in the model in [27].

In [30], a damage model was proposed for arterial tissue that includes softening and plastic phenomena. The damage in both the matrix and collagen fibers was taken into account by introducing two damage variables. Note that the strain energy function is an exponential function rather than neo-Hookean type. The strain energy function is Eq. (17). The two damage variables yield Eq. (28). This model includes the softening and plastic effects, and thus has excellent agreement with observations.

Damage models for cerebral arterial tissue have been developed [31, 32]. It was considered that the damage occurs in elastin fibers only. The isotropic model for the elastin fibers is in terms of the following strain energy function:

$$ \psi_{elastin} \left( {I_{1} } \right) = \frac{{k_{1} }}{{2k_{2} }}\left\{ {\exp \left[ {k_{2} \left( {I_{1} - 3} \right)} \right] - 1} \right\} $$
(29)

The strain energy function for fibers can be that proposed by Holzapfel et al. [23] or Gasser et al. [29]. The strain energy function with elastin damage is written as:

$$ \psi \left( {I_{1} ,I_{4} ,I_{6} ,d} \right) = \left( {1 - d} \right)\tilde{\psi }_{elsastin} \left( {I_{1} } \right) + \tilde{\psi }_{fibre} \left( {I_{4} ,I_{6} } \right) $$
(30)

where the damage variable d can be determined using two approaches, one of which is:

$$ d\left( {I_{1} } \right) = \left\{ {\begin{array}{*{20}l} {0 \quad I_{1} - 3 < \alpha_{f} } \hfill \\ {1 \quad I_{1} - 3 \ge \alpha_{f}} \hfill \\ \end{array} } \right. $$
(31)

where α f is the experimental damage threshold for elastin fibers in cerebral artery; the other damage function for cyclic loadings is expressed as:

$$ \left\{ \begin{aligned} d &= 1 - \left( {1 - d_{1} }\right)\left( {1 - d_{2} } \right)\left( {1 - d_{3} } \right) \hfill\\ d_{i} &= \left\{ \begin{array}{ll} 0 &\alpha_{\text{i}} < \alpha_{si} \hfill \\ \frac{1 - \exp \left[{c_{i} \left( 1 - {{\alpha_{i}} \mathord{\left/ {\vphantom{{\alpha_{i} } {\alpha_{fi} }}} \right. \kern-0pt} {\alpha_{fi} }}\right)} \right]}{1 - \exp \left[ {c_{i} \left( {1 - {{\alpha_{i} }\mathord{\left/ {\vphantom {{\alpha_{i} } {\alpha_{si} }}} \right.\kern-0pt} {\alpha_{si} }}} \right)}\right]} & \alpha_{i} \in [\alpha_{si} ,\alpha_{fi} ) \hfill \\ 1 & \alpha_{i} \ge \alpha_{fi} \hfill \\ \end{array} \right. \end{aligned} \right. $$
(32)

where i = 1, 2, 3, which indicate three damage mechanisms of elastin fibers in the cerebral artery, namely the damage due to maximum strain \( \alpha_{1} , \) \( \alpha_{1} \left( t \right) = {\text{Max}}\sqrt {2\tilde{\psi }\left( {I_{1} } \right)} , \) the damage due to accumulated equivalent strain \( \alpha_{2} , \) \( \alpha_{2} = \int_{0}^{t} {\left| {{{d\sqrt {2\psi \left( {I_{1} } \right)} } \mathord{\left/ {\vphantom {{d\sqrt {2\psi \left( {I_{1} } \right)} } {dt}}} \right. \kern-0pt} {dt}}} \right|} dt , \) and the damage due to the haemodynamic shearing effect on the arterial wall \( \alpha_{3} , \) \( \alpha_{3} = f\left( \textit{shears\,tress} \right) , \) where \( \alpha_{si} \) and \( \alpha_{fi} \) are the damage start and complete failure thresholds for these damage mechanisms, which need to be determined from experiments. \( c_{i} \) is another model parameter.

The work for isotropic rubber-like materials [3] described in Sect. 2.1 was extended to materials with organized fibers, for example, for the artery wall by Peña and Doblaré [33], Weisbecker et al. [19], and Pierce et al. [34]. Like Eq. (1), the strain energy function of a material with damage is written as:

$$ \psi = \psi_{vol} \left( J \right) + \sum\limits_{{i = m,_{{f_{1} ,f_{2} }} }} {[\eta_{i} \tilde{\psi }_{i} } + \phi_{i} (\eta_{i} )] $$
(33)

where \( \psi \left( J \right) \) is the strain energy function for material volume change, which has nothing to do with damage. The damage function \( \phi_{i} \left( {\eta_{i} } \right) \) is given by:

$$ - \frac{{d\phi_{i} }}{{d\eta_{i} }} = \alpha_{i} {\mathbf{erf}}^{ - 1} \left( {\beta_{i} \left( {\eta_{i} - 1} \right)} \right) + \tilde{\psi }_{i}^{0} $$
(34)

with the boundary condition \( \phi_{i} \left( 1 \right) = 0 . \) \( \tilde{\psi }_{i}^{0} \) is the strain energy function at the primary loading path. The damage variable is like that in Eq. (4):

$$ \eta_{i} = 1 - \frac{1}{{\beta_{i} }}{\mathbf{erf}}\left( {\frac{{\tilde{\psi }_{i}^{0} - \tilde{\psi }_{i} }}{{\alpha_{i} + \gamma_{i} \tilde{\psi }_{i}^{0} }}} \right) $$
(35)

where α i , β i and γ i are positive material damage constants, which need to be determined from experiments. The minimum value of η i is determined as:

$$ \eta_{i}^{0} = 1 - \frac{1}{{\beta_{i} }}{\mathbf{erf}}\left( {\frac{{\tilde{\psi }_{i}^{0} }}{{\alpha_{i} + \gamma_{i} \tilde{\psi }_{i}^{0} }}} \right) $$
(36)

where \( \eta_{i} \in \left[ {\eta_{i}^{0} ,1} \right] . \) This damage model was applied to identify the material parameters of the porcine carotid artery subjected to a uniaxial cyclic test in the longitudinal and circumferential directions [35]. The strain energy function for fibers proposed by Holzapfel et al. [36] was used. The 13 model parameters were determined using an optimization method against a series of stress–strain experimental data under various cyclic loadings.

In [37], a damage model was developed for vaginal tissue that is composed of a homogenous matrix and one family of fibers. The model is based on the assumption that damage exists in both the matrix and fibers. Introducing two damage variables, d m and d f , the damage strain energy function is written as:

$$ \psi = \psi_{vol} \left( J \right) + \left( {1 - d_{m} } \right)\tilde{\psi }_{m} + \left( {1 - d_{f} } \right)\tilde{\psi }_{f} $$
(37)

where \( \tilde{\psi }_{m} \) and \( \tilde{\psi }_{f} \) are the strain energy functions of the matrix and fibers without damage, respectively. The damage criterion resembles Eq. (26):

$$ \left\{ \begin{array}{l} \sqrt {2\tilde{\psi }_{i} \left( t\right)} - \alpha_{i} \left( t \right) \le 0 \hfill \\ \alpha_{i}\left( t \right) = {\text{Max}}\sqrt {2\tilde{\psi }_{i} \left( t\right)} \hfill \\ \end{array} \right.\quad i = m,f $$
(38)

The damage variables d m and d f can be determined using the following empirical correlation:

$$ d_{i} = \left\{ \begin{array}{ll} 0 & \alpha_{\text{i}} \left( t\right) < \alpha_{i}^{\hbox{min} } \\ \xi^{2} \left[ {1 - \beta_{i}\left( {\xi^{2} - 1} \right)} \right] & \alpha_{i} \left( t \right)\in \left[ {\alpha_{i}^{\hbox{min} },\alpha_{i}^{\hbox{max} } }\right] \hfill \\ 1& \alpha_{\text{i}} \left( t \right) >\alpha_{i}^{\hbox{max}} \hfill \\ \end{array} \right.\quad i = m,f$$
(39)

where \( \xi = {{\left[ {\alpha_{i} \left( i \right) - \alpha_{i}^{\hbox{min} } } \right]} \mathord{\left/ {\vphantom {{\left[ {\alpha_{i} \left( i \right) - \alpha_{i}^{\hbox{min} } } \right]} {\left[ {\alpha_{i}^{\hbox{max} } - \alpha_{i}^{\hbox{min} } } \right]}}} \right. \kern-0pt} {\left[ {\alpha_{i}^{\hbox{max} } - \alpha_{i}^{\hbox{min} } } \right]}} , \) \( \alpha_{i}^{\hbox{min} } \) and \( \alpha_{i}^{\hbox{max} } \) are the damage variables indicating the damage start and complete failure, respectively, for the matrix and fibers, and β i is a material parameter, \( \beta_{i} \in \left[ { - 1,1} \right] . \)

The strain energy function for the matrix material is the well-known neo-Hookean type \( \tilde{\psi }_{m} = c\left( {I_{1} - 3} \right) . \) However, the strain energy function for the one family of fibers is slightly complicated:

$$ \tilde{\psi }_{f} = \left\{\begin{array}{ll} 0 & I_{4}< I_{40} \hfill \\ \frac{{k_{1} }}{{k_{2} }}\left\{ {\exp \left[{k_{2} \left( {I_{4} - I_{40} } \right)} \right] - k_{2} \left({I_{4} - I_{40} } \right) - 1} \right\}& I_{4} > I_{40} ,I_{4}< I_{4ref} \hfill \\ 2k_{3} \sqrt {I_{4} } + k_{4} \ln \left({I_{4} } \right) + k_{5} & I_{4} > I_{4ref} \end{array} \right.$$
(40)

where \( I_{4ref} \) is the stretch squared beyond which collagen fibers start to become straightened, \( I_{40} \) is the stretch squared at which the collagen fibers begin to engage a loading. The model parameters, c, k 1 through k 5, \( \alpha_{i}^{\hbox{min} } ,\; \alpha_{i}^{\hbox{max} } , \; \beta_{i} ,\;I_{40},\) and I 4ref need to be determined from a series of experimental stress–strain curves. This model has been used to identify the parameters of vaginal tissue [37] and the rectus and sheath [38].

The continuum damage mechanics models for matrix material and fibers in [39] are very similar to those in Eqs. (37)–(39). The only difference is in the damage variable formula, and thus these modes are not further discussed below.

In [40], another type of strain energy function for fibers and the continuum mechanics damage variable d were used in the damage model. Since this model just combines previous work, it is not discussed here.

The damage model proposed in [28] was updated in [41] by introducing the proteoglycan bridge damage of collagen fibrils. The bridge damage is modeled with statistical processes and related to the damage variable d. The statistical proteoglycan orientation (beta distribution) and bridge internal length (Gaussian distribution) contribute to the bridge damage process.

In the damage models mentioned above, a standard continuum damage formulation and damage variable are used, because the damage effect in a soft tissue is a natural result of cross-section reduction of a specimen of the tissue. This means that the material strength degradation is in a local sense. Such a local effect can result in an ill-posed problem and increased mesh refinement, especially for soft tissues usually subjected to a significant deformation [42]. To overcome this geometrically non-linear effect, a non-local gradient damage formulation has been presented [42]. In this formulation, two extra energy functions were added into the strain energy function, proposed in [29], to include the damage effect. The first extra energy function is the scalar product of the gradient of a non-local damage variable and a scalar field variable with respect to three coordinates in the current configuration; the second extra energy function is the penalty function of the squared diffidence between the scalar field variable and the usual damage variable. Based on the principle of minimum potential energy, a second-order partial differential equation of the non-local damage variable was established. The usual damage variable was the source term of that partial differential equation. Like in other damage models, the stress in the fibers is degraded by making use of the ordinary exponential damage function of the usual damage variable. Since an additional non-local damage variable has to be solved during damage simulation, this approach may be time-consuming.

3.2 Probabilistic Damage Models of Fibers

A damage model was proposed by Chu and Blatz [43] for cat mesentery. The stress–stretch curves of biaxial test specimens of cat mesentery exhibit hysteresis. This is mainly due the cumulative microdamage mechanism that occurs in the collagen fibers. To model this effect, it is assumed that the stretch in a region is statistically distributed among fibers and that the stresses are statistically distributed among the remaining unbroken fibers. Based on the Ogden-type strain energy function in Eq. (2), the curves were fitted by adjusting the property constants in the strain energy function and the parameters in the probability density function of stress. This may be the very first damage model for soft tissue.

In [39, 44, 45], a stochastic damage model for fibers was proposed and a comparison of the predictive capability between the model and the continuum damage mechanics model was made, with similar results obtained. In the stochastic damage model, the damage model for the matrix material is an ordinary one, like those mentioned above; but the fiber damage model is slightly different. However, the strain energy function of individual fibers in [44] is quite different from that in [39, 45]. The proper choice of function is not clear. In the following descriptions, the specific form of the strain energy function for a fiber is thus omitted.

In the loading-free state, a collagen fiber is wavy. With increasing loading, it starts to become straight until it fails. The total strain energy function of all fibers expressed by Rodriguez et al. [45] is:

$$ \tilde{\psi }_{f} \left( \lambda \right) = \left\{\begin{array}{ll} 0& \lambda < 1 \hfill \\ \int_{0}^{\lambda}\int_{{l_{\hbox{max} } }}^{l} {\sigma_{f} \left( {\xi ,x}\right)}p\left( x \right)dxd\xi & \lambda \ge 1 \hfill \\ \end{array} \right. $$
(41)

where σ f is the stress in a fiber, \( \sigma_{f} = {{\partial \tilde{\psi }_{1} } \mathord{/ {\vphantom {{\partial \tilde{\psi }_{1} } {\partial \lambda }}} \kern-0pt} {\partial \lambda }} , \; \tilde{\psi }_{1} \) is the strain energy function of a fiber that takes a form based on eight-chain model proposed by Arruda and Boyce [46] for rubber elastic materials, and p(x) is a beta probability density function with parameters m and n to account for stress variation among fibers:

$$ p\left( x \right) = \frac{1}{{l_{\lim } - l_{0} }}\frac{{\Upgamma \left( {m + n} \right)}}{\Upgamma \left( m \right)\Upgamma \left( n \right)}\left( {\frac{{x - l_{0} }}{{l_{\lim } - l_{0} }}} \right)^{n - 1} \left( {1 - \frac{{x - l_{0} }}{{l_{\lim } - l_{0} }}} \right)^{m - 1} \, x \in \left[ {l_{0} ,l_{\lim } } \right] $$
(42)

where l 0 is the reference length of fibers, \( l_{\hbox{max} } = \exp \left[ {\left( {{x \mathord{\left/ {\vphantom {x \delta }} \right. \kern-0pt} \delta }} \right)^{\theta } } \right]l_{0} \lambda_{\hbox{max} } ,\; \lambda_{\hbox{max} }\) is the maximum stretch of fibers without failure over the time history in a test (if \( \lambda \ge \lambda_{\hbox{max} } \) then a failure starts to occur), δ and θ are the model parameters, and l lim is the failure stretch limit (when a stretch reaches this value, the fiber breaks completely).

A ligament is a soft tissue that connects two bones in a joint. Ligaments stabilize joints and guide their motion when a tensile load is applied [47]. Ligaments are composed of collagens (approximately 75 % of the dry weight), proteoglycans (<1 %), elastin, other proteins (glycoproteins, such as actin, laminin, and integrins), and some water, which may be responsible for cellular function and viscoelastic behavior [47, 48].

The ligament microstructure comprises collagen bundles aligned along the long axis of the ligament and exhibits a wavy or crimp pattern along the length, which allows the ligament to elongate without sustaining damage after a load is applied.

Ligaments demonstrate passive nonlinear anisotropic biomechanical behavior only. At a low loading, they are relatively compliant because of crimped collagen fibers and the viscoelastic effect; at a high loading, however, they are much stiffer because fibers are recruited and straightened [47]. Ligaments are quite often damaged in traumatic joint injuries. The damage includes partial ligament failure or complete ligament tear [47]. The ligament stress–stretch curve has an S-shape.

Based on the recruitment models of collagen fibers in [49–51], see Appendix 1 for details, a probabilistic damage model for ligaments was developed by Liao and Belkoff [52]. To derive the damage model, several assumptions were made: (1) only the fibers in the ligament respond to a loading; (2) the interaction between the matrix and fibers and that among fibers are ignored; (3) the viscoelastic effect of a ligament is not taken into account; (4) the initial shape of fibers is wavy in the stress-free state and has a Gaussian distribution; (5) all the fibers are linearly elastic and have the same elastic modulus and limit strain; (6) fibers experience brittle failure and they fail in the same sequence in which they are recruited. The number of fibers recruited is given by the following expression:

$$ dn = \frac{n}{{s\sqrt {2\pi } }}\exp \left[ { - \frac{1}{2}\left( {\frac{{x - \overline{x} }}{s}} \right)^{2} } \right] $$
(43)

where n is the total number of fibers. The total force generated by all the fibers when they are subjected to a stretch is given by the equation:

$$ F\left( \lambda \right) = \frac{{nA_{i} E_{i} }}{{s\sqrt {2\pi } }}\int_{1}^{\lambda } {\left( {\frac{\lambda - x}{x}} \right)} \exp \left[ { - \frac{1}{2}\left( {\frac{{x - \overline{x} }}{s}} \right)^{2} } \right]dx $$
(44)

where A i and E i are the cross-section and Young’s modulus of a fiber, respectively. Then, the mean stress in the fibers is:

$$ \sigma \left( \lambda \right) = \frac{F\left( \lambda \right)}{{nA_{i} }} = \frac{{E_{i} }}{{s\sqrt {2\pi } }}\int_{1}^{\lambda } {\left( {\frac{\lambda - x}{x}} \right)} \exp \left[ { - \frac{1}{2}\left( {\frac{{x - \overline{x} }}{s}} \right)^{2} } \right]dx $$
(45)

At a breaking strain ε lim, or a breaking stretch \( \lambda_{\lim } = 1 + \varepsilon_{\lim}, \) some straightened fibers fail. The stress due to the contributions of fibers stretched beyond \( \lambda /\lambda_{\lim } \) is:

$$ \sigma \left( {{\lambda \mathord{\left/ {\vphantom {\lambda {\lambda_{\lim } }}} \right. \kern-0pt} {\lambda_{\lim } }} \ge 1} \right) = \frac{{E_{i} }}{{s\sqrt {2\pi } }}\int_{1}^{{{\lambda \mathord{\left/ {\vphantom {\lambda {\lambda_{\lim } }}} \right. \kern-0pt} {\lambda_{\lim } }}}} {\left( {\frac{\lambda - x}{x}} \right)} \exp \left[ { - \frac{1}{2}\left( {\frac{{x - \overline{x} }}{s}} \right)^{2} } \right]dx $$
(46)

The resultant stress after failure should be equal to \( \sigma \left( \lambda \right) - \sigma \left( {{\lambda \mathord{\left/ {\vphantom {\lambda {\lambda_{\lim } }}} \right. \kern-0pt} {\lambda_{\lim } }} \ge 1} \right) , \) and thus we have the stress after failure as:

$$ \sigma \left( \lambda \right) = \frac{{E_{i} }}{{s\sqrt {2\pi } }}\int_{{{\lambda \mathord{\left/ {\vphantom {\lambda {\lambda_{\lim } }}} \right. \kern-0pt} {\lambda_{\lim } }}}}^{\lambda } {\left( {\frac{\lambda - x}{x}} \right)} \exp \left[ { - \frac{1}{2}\left( {\frac{{x - \overline{x} }}{s}} \right)^{2} } \right]dx \, \quad \lambda \ge \lambda_{\lim }. $$
(47)

Before failure, i.e., \( \lambda \le \lambda_{\lim } , \) the stress can be estimated using Eq. (47). The parameters E i , \( \overline{x} , \) s, and λ lim need to be determined from experimental stress–stretch curves.

Another idea for modeling ligament rupture failure is that the damage of a ligament is a gradual reduction of its fiber stiffness at a randomly distributed stretch threshold rather than a constant limit strain or stretch. This type of damage model was proposed in [53, 54] based on a stretch threshold described by the Weibull distribution, which is frequently used to describe the random yield strength or fatigue life of a material [55].

In the model, the collagen fibers carry the load applied on a ligament, and the elastin contribution is neglected. The interaction between fibers and the matrix material is not taken into account. The collagen fibers are linearly elastic and their orientation is parallel to the loading direction. The fibers are wavy in the stress-free sate, but they become straightened with increasing loading until damage occurs. The fiber recruitment effect is ignored, however, and the collagen straightening stretch λ s and stretch threshold for damage, λ f , are considered as statistical variables specified by a Weibull distribution.

The Weibull probability distribution function for the collagen straightening stretch λ s at which fibers are straight and ready to engage a load is:

$$ P_{s} \left( \lambda_{s},\alpha_{s} ,\beta_{s} ,\gamma_{s}\right) = \left\{ \begin{array}{ll} 0 & for \; \lambda_{\text{s}}< \gamma_{s} \hfill \\ 1 - \exp \left[ - \left(\frac{\lambda_{s} - \gamma_{s}}{\beta_{s}} \right)^{{\alpha_{s}}} \right] & for \; \lambda_{\text{s}} < \gamma_{s} \hfill \\ \end{array} \right. $$
(48)

where the parameter \( \gamma_{s} = 1, \) and the other two positive parameters, α s and β s , are determined based on the microstructure of a specimen in the stress-free state or from the stress–stretch curve. The inverse function of Eq. (48) is:

$$ \lambda_{s} (P_{s} ;\alpha_{s} ,\beta_{s} ,\gamma_{s} ) = \gamma_{s} + \beta_{s} \left[ { - \ln \left( {1 - P_{s} } \right)} \right]^{{{1 \mathord{\left/ {\vphantom {1 {\alpha_{s} }}} \right. \kern-0pt} {\alpha_{s} }}}} $$
(49)

Let us assume that there are n fibers in a specimen. The initial length of the i-th fiber in the specimen, \( \lambda_{s}^{(i)} , \) can be calculated as:

$$ \lambda_{s}^{(i)} (P_{s}^{(i)} ;\alpha_{s} ,\beta_{s} ,\gamma_{s} ) = \gamma_{s} + \beta_{s} \left[ { - \ln \left( {1 - P_{s}^{(i)} } \right)} \right]^{{{1 \mathord{\left/ {\vphantom {1 {\alpha_{s} }}} \right. \kern-0pt} {\alpha_{s} }}}} $$
(50)

where \( P_{s}^{(i)} \) is a random number between 0 and 1.

It is assumed that fiber damage occurs in m fibers once a limit stretch is exceeded, which is a random variable described by a Weibull function. This means that damage consists of a series of sub-failures. Thus, the j-th sub-failure of the i-th fiber is expressed by the following Weibull distribution:

$$ \lambda_{d}^{(j)} (P_{d}^{(j)} ;\alpha_{d} ,\beta_{d} ,\gamma_{d} ) = \gamma_{d} + \beta_{d} \left[ { - \ln \left( {1 - P_{d}^{(j)} } \right)} \right]^{{{1 \mathord{\left/ {\vphantom {1 {\alpha_{d} }}} \right. \kern-0pt} {\alpha_{d} }}}} $$
(51)

Similarly, \( P_{d}^{(j)} \) is a random number between 0 and 1; \( \gamma_{d} \) is the mean limit stretch and can be determined easily from experiments. For example, for rat medial collateral ligaments, \( \gamma_{d} = 1.0514 \) [56]. The remaining positive parameters α d and β d need to be optimized from experimental data.

All fibers have the same Young’s modulus, k. When damage occurs, the fiber modulus will degrade gradually until a complete tear. It is assumed that the damage process conforms to an exponential law [54]. The damage variable d is related to the stress by the follow expression:

$$ \sigma^{(i)} \left( {k,d,\lambda^{(i)} ,\lambda_{s}^{(i)},\lambda_{d}^{(j)} } \right) = \left\{ \begin{array}{ll} 0 &\lambda^{(i)} \le \lambda_{s}^{(i)} \hfill \\ k\left({\lambda^{(i)} - \lambda_{s}^{(i)} } \right) & \lambda^{\left( i\right)} \in \left( {\lambda_{s}^{(i)} ,\lambda_{d}^{\left( j\right)} } \right) \hfill \\ d^{j} k\left( \lambda^{(i)} -\lambda_{d}^{(j)} \right) & \lambda^{\left( i \right)} \ge \lambda_{d}^{(j)} \hfill \\ \end{array} \right. $$
(52)

where the damage variable d ∈ (0,1) and modulus k need to be determined from experimental stress–stretch curves.

Finally, the stress in the specimen is defined as the mean of the stresses of n fibers, expressed by:

$$ \sigma \left( \lambda \right) = \frac{1}{n}\sum\limits_{i = 1}^{n} {\sigma^{(i)} \left( {k,d,\lambda^{(i)} ,\lambda_{s}^{(i)} ,\lambda_{d}^{(j)} } \right)} $$
(53)

The model above is subject to the set of parameters \( \left\{ {k,d,\alpha_{s} ,\beta_{s} ,\gamma_{s} ,\alpha_{d} ,\beta_{d} ,\gamma_{d,} n,m} \right\} . \) Since \( \gamma_{s} = 1, \) \( \gamma_{d} = 1.0514, \) n = 105, and m = 102 [54], the model requires six parameters: \( \left\{ {k,d,\alpha_{s} ,\beta_{s} ,\alpha_{d} ,\beta_{d} } \right\} . \)

The early version of the model above is also interesting [53]. In that model, there is no damage variable and it is supposed that once a limit stretch is exceeded, a fiber fails rather than breaking step by step. Moreover, it is assumed that the limit stretch itself is a random variable described by the Weibull function:

$$ \lambda_{f}^{(i)} = \gamma_{f} + \beta_{f} \left[ { - \ln \left( {1 - P_{f}^{\left( i \right)} } \right)} \right]^{{{1 \mathord{\left/ {\vphantom {1 {\alpha_{f} }}} \right. \kern-0pt} {\alpha_{f} }}}} $$
(54)

where \( \lambda_{f}^{\left( i \right)} \) is the failure stretch of the i-th fiber. Accordingly, the stress in that fiber is given by:

$$ \sigma^{(i)} \left( {k,\lambda^{(i)} ,\lambda_{s}^{(i)},\lambda_{f}^{(j)} } \right) = \left\{\begin{array}{ll} 0 & \lambda^{(i)} \le \lambda_{s}^{(i)} \hfill \\ k\left({\lambda^{(i)} - \lambda_{s}^{(i)} } \right) & \lambda^{\left( i\right)} \in \left(\lambda_{s}^{(i)} ,\lambda_{f}^{\left( j \right)}\right) \hfill \\ 0 & \lambda^{\left( i \right)} \ge \lambda_{f}^{(j)} \hfill \\ \end{array} \right. $$
(55)

This model is subject to five parameters \( \left\{ {k,\alpha_{s} ,\beta_{s} ,\alpha_{f} ,\beta_{f} } \right\} \) only and its performance is very satisfactory.

The tendon is a connecting tissue between muscles and bone and can be damaged after being excessively stretched. The tendon is mainly composed of collagen fibers, some proteoglycans, and fluid, and thus it can respond to a loading passively. The mechanical properties of tendons were investigated in vitro by Schechtman and Bader [57, 58]. A typical stress–strain relationship obtained from a simple tensile test [57] has an S-shape.

It is considered that a human tendon is a collagen-fiber-reinforced composite nonlinear material with a uniform matrix that may be compressible [59]. The fibers are wavy in the load-free state; if a load is applied in the physiological force direction, they are stretched but do not generate passive tension until straightened.

The biomechanical interaction between the matrix material and fibers as well as the viscous effect (strain-rate-dependent feature) are not taken into account. In [59], the damage model is based on the following strain energy function:

$$ \psi ({\mathbf{C}}) = U(J) + \psi_{m} (\bar{I}_{1} ,\bar{I}_{2} ) + \psi_{f} (I_{4} ) $$
(56)

where C is the right Cauchy-Green tensor, \( {\mathbf{C}} = {\mathbf{F}}^{\mathbf{T}} {\mathbf{F}} , \) \( \bar{I}_{1} \) and \( \bar{I}_{2} \) are the principal invariants associated with the iso-volumetric components of the right Cauchy-green tensor, \( {\bar{\mathbf{C}}} = J^{{ - {2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-0pt} 3}}} {\mathbf{C}} , \) and \( I_{4} \) is the squared stretch along the fiber orientation:

$$ I_{4} = a_{0} \cdot {\mathbf{C}}a_{0} = \lambda^{2} $$
(57)

where a 0 is the vector of fiber orientation and λ is the stretch.

The matrix material volume term U(J) is defined as:

$$ U(J) = K(J - 1)^{2} $$
(58)

where K is the modulus of the material (K = 1000 MPa) [59]. The iso-volumetric term of the matrix material is expressed as:

$$ \psi_{m} (\bar{I}_{1} ,\bar{I}_{2} ) = c_{1} (\bar{I}_{1} - 3) + c_{2} (\bar{I}_{2} - 3) $$
(59)

where the property constants \( c_{1} \) and \( c_{2} \) are the initial tangent shear modulus of the matrix material and need to be determined from experimental data. The strain energy function for fiber response is written as [59]:

$$ \tilde{\psi }_{f} (I_{4} ) = \frac{{k_{1} }}{{k_{2} }}\left[ {e^{{k_{2} (I_{4} - 1)}} - k_{2} (I_{4} - 1) - 1} \right] $$
(60)

where k 2 is the parameter associated with the initial crimp of the fiber and k 1 is the initial stiffness of the fiber. Equation (60) requires \( I_{4} \ge 1 . \)

For a tendon, damage occurs to fibers only (the matrix material is not damaged). The fiber damage function is associated with the last term \( \psi_{f} (I_{4} ) \) in Eq. (56). The stain energy function with the fiber damage effect takes the following form [59]:

$$ \psi ({\mathbf{C}},d) = U(J) + \psi_{m} (\bar{I}_{1} ,\bar{I}_{2} ) + g(d)\tilde{\psi }_{f} (I_{4} ) $$
(61)

The fiber damage function is related to the fiber stretch and expressed as [59]:

$$ g(d) = \frac{{1 - e^{{\beta (\lambda^{4} - \lambda_{\lim }^{4} )}} }}{{1 - e^{{\beta (\lambda_{{^{0} }}^{4} - \lambda_{\lim }^{4} )}} }} $$
(62)

where β is a parameter related to the initial wavy state of the fiber, β < 0, λ 0 is the maximum stretch without fiber damage, and λ lim is the fiber stretch at which all fibers are broken. If λ is between λ 0 and λ lim, then the fiber damage function, g(d), will be engaged in Eq. (62); otherwise, g(d) = 1.

The damage variable d is defined as the ratio of the number of damaged fibers, n b , to the total number of fibers, n, in a test sample of tissue:

$$ d = \frac{{n_{b} (\lambda )}}{n} $$
(63)

Since the crimped state of fibers varies with the sample, each fiber should start to be damaged at its own critical stretch. It is assumed that such an effect is in accordance with the Gaussian distribution [49]. The number of damaged fibers, \( n_{b} , \) is expressed as:

$$ n_{b} (\lambda ) = \frac{n}{{s\sqrt {2\pi } }}\int_{{\lambda_{0} }}^{\lambda } {\text{e}^{{ - \frac{{(x - \bar{\lambda })^{2} }}{{2s^{2} }}}} dx} $$
(64)

where \( \bar{\lambda } \) and s are the mean value and standard deviation of the critical stretch for damage of an individual fiber, respectively. They are determined from experimental data. Eventually, the damage variable d takes the following form:

$$ d(\lambda ) = \frac{1}{{s\sqrt {2\pi } }}\int_{{\lambda_{0} }}^{\lambda } {\text{e}^{{ - \frac{{(x - \bar{\lambda })^{2} }}{{2s^{2} }}}} dx} $$
(65)

It is related to the stretch during damage by a linear equation, i.e.:

$$ \lambda = \lambda_{0} + (\lambda_{1} - \lambda_{0} )d $$
(66)

If a series of experimental stretch data and the values of \( \lambda_{0} \) and \( \lambda_{1} \) are available, Eq. (66) can be used to calculate a series of \( d \); then, \( \bar{\lambda } \) and \( s \) can be determined by fitting the scattered \( d - \lambda \) plot with Eq. (65). Eventually, the relationship among \( g(d) , \) \( d , \) and \( \lambda \) can be obtained. The parameters \( c_{1} , \) \( c_{2} , \) and \( k_{1} \) are obtained from experimental stress–strain data. The parameters \( k_{2} \) and \( \beta \) are associated with the wavy state of fibers under the load-free condition. For a cyclic loading, k 2 = 10 and β = −2.8 (fibers are very wavy); for a steadily increasing loading, k 2 = 25 and β = −1.0 (fibers are less wavy) [59].

Theoretically, I 4 in Eq. (56) should be \( \bar{I}_{4} , \) which is associated with the iso-volumetric components of the right Cauchy-Green tensor, \( {\bar{\mathbf{C}}} , \) i.e., \( \overline{I}_{4} = a_{0} \cdot \overline{\mathbf{C}} a_{0} = \lambda^{2} \) [23]. Therefore, the justification of the damage models expressed with Eq. (61) needs to be clarified in the future.

3.3 Microstructure-Based Damage Model of Fibers

Studies have investigated constitutive models of vascular tissue and its damage modeling [60, 61] based on the microstructure of collagen fibrils in vascular wall tissue. The major contents of vascular tissue are elastin, collagen, and proteoglycans; in particular, collagen fibers play a very important role in determining the biomechanical properties of the tissue. It was shown that a series of proteoglycan (PG) bridges can be formed to generate force as soon as collagen fibrils become straightened (λ st  = 1) by stretch. The first Piola–Kirchhoff stress generated in a fibril is given as the following equation by employing a triangle probability density function [61]:

$$ T\left(\lambda \right) = \left\{ \begin{array}{ll} 0 & \lambda\in \left(0,\lambda_{\hbox{min}} \right] \hfill \\ \frac{2k}{{3\Delta \lambda^{2} }}\left(\lambda -\lambda_{\hbox{min}}\right)^{3} & \lambda \in \left(\lambda_{\hbox{min}} ,\overline{\lambda} \right] \hfill \\ k\left[\lambda - \tfrac{{2\left( {\lambda - \lambda_{\hbox{max} } }\right)^{3} }}{{3\Delta \lambda^{2} }} - \overline{\lambda}\right]&\lambda \in \left(\overline{\lambda } ,\lambda_{\hbox{max}} \right]\hfill \\ k\left( \lambda - \overline{\lambda} \right) & \lambda \in \left(\lambda_{\hbox{max} } ,\infty \right) \hfill \\ \end{array} \right. $$
(67)

where \( \Delta \lambda = \lambda_{\hbox{max} } - \lambda_{\hbox{min} } , \) \( \overline{\lambda } = 0.5\left( {\lambda_{\hbox{min} } + \lambda_{\hbox{max} } } \right) , \) and \( k \) is the stiffness of a collagen fibril. The Cauchy stress of a collagen fibril is \( \sigma \left( \lambda \right) = \lambda T\left( \lambda \right) \); usually, \( \lambda_{st} = \lambda_{\hbox{min} } = 1, \) and λ max = 2. The Cauchy stress in the fibers of tissue is expressed as:

$$ \sigma_{fibre} = \int\limits_{V} {\rho \left( {\text{N}} \right)} \sigma \left( {\lambda \left( {\text{N}} \right)} \right)dev\left( {{\text{n}} \otimes {\text{n}}} \right)dV $$
(68)

where \( \rho \left( {\text{N}} \right) \) is the orientation density function of fiber bundles, \( {\text{n}} \) is the direction vector of a fibril, and \( V \) is the total volume of a vascular tissue. The total Cauchy stress in the tissue is calculated as:

$$ \sigma = \sigma_{vol} + \sigma_{nH} + \sigma_{fibre} $$
(69)

where \( \sigma_{vol} \) and \( \sigma_{nH} \) are attributed to the volumetric energy function \( \psi_{vol} = K\left( {J - 1} \right)^{2} , \) and the neo-Hookean strain energy function. Such a treatment for the stress in fibers involves the collagen recruitment concept.

The damage is assumed to occur in fibers only due to a load [60]. A damage variable is involved in the second Piola–Kirchhoff stress of a collagen fibril:

$$ S = \left( {1 - d} \right)\tilde{S} = \left( {1 - d} \right)k\left( {{\lambda \mathord{\left/ {\vphantom {\lambda {\lambda_{st} - 1}}} \right. \kern-0pt} {\lambda_{st} - 1}}} \right) $$
(70)

where the damage variable \( d = 1 - \exp \left[ { - a\left( {{{\lambda_{st} } \mathord{\left/ {\vphantom {{\lambda_{st} } {\lambda_{st0} }}} \right. \kern-0pt} {\lambda_{st0} }}} \right)^{2} } \right] , \) \( a \) is a model parameter, \( \lambda_{st0} \) is the stretch of a fibril becoming straightened initially, and \( \lambda_{st} \) is the stretch of a fibril becoming straightened later. Because of plastic deformation, the relation \( \lambda_{st} \le \lambda_{st0} \) is kept. The damage criterion for fibers is written as [60]:

$$ \left\{ \begin{array}{ll} \tilde{S} < Y_{0} & {\text{ elastic deformation, no damage}} \hfill \\ \tilde{S} = Y_{0} + H & {\text{ plastic deformation, damage}} \hfill \\ \tilde{S} > Y_{0} + H & {\text{ complete damage}} \hfill \\ \end{array} \right. $$
(71)

where \( \tilde{S} \) is the second Piola–Kirchhoff stress without damage, \( Y_{0} \) is the elastic limit, \( H \) represents the hardening effect due to the slowly (viscous) sliding mechanism in proteoglycan bridges, \( H = \eta {{d\lambda_{st} } \mathord{\left/ {\vphantom {{d\lambda_{st} } {dt}}} \right. \kern-0pt} {dt}} , \) and \( \eta \) is an experimental coefficient. \( \lambda_{st} \) increases with time \( t . \) The Cauchy stress in a fibril is obtained as \( \sigma = J^{ - 1} S\lambda^{2} . \) Equations (67) to (71) represent the damage model for collagen fibers in vascular tissue based on their microstructure proposed by Gasser [61]. Strictly speaking, this damage model is probabilistic.

This damage model relies on the irreversible sliding damage of PG bridges across collagen fibrils. It was shown that PG bridges exist in cartilage and tendons [62–64]. Whether they exist in vascular tissue needs to be confirmed by microscopic observation. In addition, the sliding damage effect or plastic deformation needs to be clarified at microscopic experimental level.

4 Discussion and Challenges

In the reviewed damage models, there are a few parameters that are determined using a series of experimental data for a steady or cyclic load. Usually, they are determined mathematically by means of optimization or the least squares method under the condition that the error in the stress–stretch curve between measured and predicted values is the minimum.

The experimental data can be from uniaxial or biaxial tensile or shear tests or inflation measurements of a segment of organ or soft tissue. For organs, FEA is used to get the Cauchy stresses for a certain load, which is a time-consuming optimization procedure.

4.1 Experimental Data Curve-Fitting

In this section, a few important cases are given to illustrate the feasibility of damage models. In order to show the discontinuous damage effect, the results predicted by the damage model proposed in [27] are demonstrated. The arterial tissue of animals exhibits softening behavior during a uniaxial tensile test under cyclic loads [27]. This behavior has been investigated theoretically based on the damage model in Eqs. (22)–(28). The model parameters are μ = 0.072294 MPa, d m  = d 4 = d 6 = 0.00006, \( \zeta_{i}^{\hbox{min} } \) = 0.08, \( \zeta_{i}^{\hbox{max} } = 0.41, \) \( \xi_{i} = 0.61, \) \( d_{i\infty } = 0.507, \) and \( \gamma_{i} = 15.0 \) for the matrix, and k 1 = 0.000582 MPa, k 2 = 3.675678, \( \zeta_{i}^{\hbox{min} } = 0.01, \) \( \zeta_{i}^{\hbox{max} } = 0.465, \) \( \xi_{i} = 1.06, \) d 4 = d 6 = 0.00006, \( d_{i\infty } = 0.507, \) and \( \gamma_{i} = 2.55 \) for the fibers. A combination of discontinuous and continuous damage models can predict the softening/damage effect better than can either model alone.

The isotropic damage model in [3] has been extended to anisotropic cases. This extended damage model has potential applications in biomedical engineering [19, 33, 34]. Figure 3 shows the experimental and predicted Cauchy stress–stretch curves and damage variable variation under cyclic loads for arteries given by Weisbecker et al. [19]. The experimental data are fitted very well.

Fig. 3
figure 3

Experimental and predicted Cauchy stress-stretch curves and damage variable variation under cyclic loads. a, b Damage of collagen fibers, c, d damage of media, from [19]

Full size image

Figure 4 shows the predicted stress–stretch curve of the ligaments harvested from two groups of rabbits based on the probabilistic damage model proposed by Liao and Belkoff [52]. It can be seen that the abrupt failure behavior of ligaments is captured very well. Figure 4 also shows the performance of models proposed by De Vita and Slaughter [53] and Guo and De Vita [54] for freshly harvested rat medial collateral ligaments. Although both models produce an S-shape curve, the stress–stretch curve obtained by Guo and De Vita [54] is not smooth enough.

Fig. 4
figure 4

Comparison of experimental stress-stretch data of ligaments and their prediction made by damage models proposed respectively by a, b Liao and Belkoff [52], c Guo and De Vita [54], and d De Vita and Slaughter [53]. Symbols are experimental data and lines are model prediction results

Full size image

The damage model proposed by Natali et al. [59] was used to represent experimental data of human tendons before and after cyclic loadings. Two curves are shown in Fig. 5, where the model parameters in Eqs. (56)–(66) are K = 1000 MPa, c 1 = 1.0 MPa, c 2 = 2.0 MPa, k 1 = 4.0 MPa, k 2 = 10, λ 0 = 1.09, and λ 1 = 1.195 before cyclic loading, and k 2 = 25, λ 0 = 1.02, and λ lim = 1.08 after cyclic loading. The S-shape curve is well predicted by the model; however, the curve is not as sharp as the experimental data.

Fig. 5
figure 5

Stress–stretch curves for human tendon before and after cyclic loading, from [59]

Full size image

Based on the comparisons above, it can be concluded that existing damage models can represent experimental stress–stretch data precisely.

4.2 Computational Effort

The damage model proposed by Volokh [24, 25] is based on the strain energy function for the arterial wall proposed by Holzapfel et al. [23]. The softening effect is handled using the softening parameter \( \phi , \) ξ 4(=ξ 6), and n 4(=n 6) or strain energy limiters and sharpness factors. Therefore, a cyclic loading is unnecessary. Hence, this damage model can be readily applied to the analysis of simple tensile test results of soft tissue. Moreover, this model can be easily handled with MATLAB code for simple structures, such as a tube.

For damage models with damage variables, experimental stress–stretch curves of soft tissue are required for determining the damage parameters. Usually, this kind of damage model is so complex that a FORTRAN UMAT subroutine is needed for ABAQUS, ANSYS, or other FEA packages to perform a FEA in them. If so, the first-order derivatives of the Cauchy stresses with respect to stretch components and damage variables are desirable. Since these derivatives cannot be expressed analytically, a numerical increment formulation is required.

In the damage models for ligaments and tendons, the experimental stress–stretch curves can be well predicted by employing the fiber recruitment effect and probabilistic strain failure limit. The limitations in these models are that damage exists in fibers only and that the fibers are linearly elastic. To remove these drawbacks, proper strain energy functions for the matrix and fibers should be developed and the damage mechanisms should be introduced into the matrix and fibers as well. Moreover, probabilistic damage models require extremely long computational time and their application to complex structures may not be easy.

4.3 Application of Damage Models

The ultimate objective of generating soft tissue damage models is to apply them to disease diagnosis, surgery, surgeon training procedures, and the design and fabrication of artificial soft tissue. Even though many damage models have been available for soft tissues, their applications in biomedical engineering seem limited.

Figure 6 shows the damage variable distribution on a human arterial media inner surface presented by Schmidt et al. [41] based on their own damage model with parameters extracted from circumferential and longitudinal uniaxial tests for human carotid artery media. It can be observed that the media experiences a serious damage effect around the fibrous cap.

Fig. 6
figure 6

Damage variable distribution in human arterial wall under 80 kPa internal blood pressure. a 3D diseased arterial model and b damage variable distribution on wall, from [41]

Full size image

The strength, von Mises stress, rupture potential index, and damage variable distributions on 10 patient-specific AAA walls under various blood pressures based on a damage model were reported by Marini et al. [11]. The damage developed in the areas with a high von Mises stress and a large rupture potential index. Even though the damage variable is correlated to AAA rupture slightly more poorly than to the von Mises stress, the damage variable still provides a useful link between a mechanical stimulus and the response of an AAA wall.

The applications of damage models in biomechanical engineering are exciting and convincing. Hence, more trials should be conducted in the future.

4.4 Connection Between Damage and Fracture

Damage is closely related to fracture or crack extension/propagation in a soft tissue [4]. However, existing damage models are based on continuum damage mechanics and are macrostructure-based. The model parameters are usually obtained by fitting stress–stretch curves without any information about cracks. For soft tissues, a link between damage variables and crack propagation has not been established.

For fiber-reinforced soft tissues, a sub-failure or complete failure is driven by crack propagation in brittle failure or toughening in ductile failure inside a material. A visualization study of crack development or toughening during a simple tensile test needs to be conducted for soft tissues [65]. For a long-term objective, the crack characteristics or toughening behavior should be linked to the damage time-history profile and the damage behavior of soft tissues. For the visualization of cracking in soft tissue, damage patterns such as matrix cracking, fiber bridging, fiber rupture, fiber pull-out, and matrix/fiber de-bounding should be considered [66].

Numerical simulation of fracture propagation in soft tissue is equally important. Cohesive-zone law models are considered effective tools for tracking macro-crack propagation in a solid material under time-dependent loadings [67, 68]. As shown in Fig. 7, from Ferrara, Pandofi [69], with the cohesive-zone law model, crack generation and development in the arterial wall can be identified very clearly. With increasing inner blood pressure, more cracks are generated, and the existing cracks open widely and propagate deeply in the tissue.

Fig. 7
figure 7

Crack generation and development in arterial wall under increasing inner blood pressure. a 100, b 120, c 180, and d 260 mmHg, from [69]

Full size image

Another interesting study is the propagation of arterial dissection, which is frequently performed in clinical practice and can be caused by traffic accidents. A three-dimensional (3D) isotropic cohesive model was proposed by Gasser and Holzapfel [70] based on cohesive potential for human aortic media. The model involves cohesive tensile strength, two non-negative model parameters, and a damage variable that is the magnitude of the opening gap displacement of a crack. The model parameters are determined from the load-gap displacement curve obtained in a media dissection experiment. The radial Cauchy stress distributions predicted by the model during a dissection process of a two-dimensional (2D) human aortic media strip are shown in Fig. 8. This model can potentially be applied to the dissection of a 3D human aortic artery.

Fig. 8
figure 8

Predicted dissection process of 2D human aortic media obtained using isotropic cohesive model in [70], images from [70]

Full size image

The material point method (MPM), a numerical method, can potentially be applied to track the cracks in soft tissues and characterize the fracture failure behavior of the tissue. MPM has advantages over traditional FEA methods since it can adapt to complex geometry, large deformation, and fragmentations that may occur in the fracture failure of soft tissue [71, 72].

4.5 Other Issues

The viscoelastic effect may not be significant in damage to soft tissue [73]. However, a study showed that viscoelasticity is important for the damage modeling of elastic and viscoelastic materials [74]. Nevertheless, this problem needs to be clarified in detail. A method for estimating the dissipated energy via viscoelasticity in soft tissue has been proposed and compared with other methods [75].

It has been shown that discontinuous and continuous damage mechanisms are equally important for softening effect prediction in soft tissues [27, 28]. Nevertheless, more experimental evidence is required.

Soft tissue damage is relevant to tissue histological change and inflammation of cells. A well-defined relationship between the histological change and the inflammation of cells is unavailable. For smooth muscle, however, active stress is dominant but is not considered in any damage model.

It has been indicated that the soft tissue of left and right ventricles of animals also exhibits the softening effect under a cyclic loading [76, 77]. The softening effect with quite substantially plastic deformation was found for mouse skin [78] and ovine infrarenal vena cava tissue [79]. However, no damage model for myocardium or skin has been proposed.

Very soft tissues, such as those of the brain, liver, kidneys, and even skin exhibit a strong viscoelastic property and transversely isotropic behavior [73, 80–82]. Damage modeling for this kind of soft tissue is very important for automatic surgical tools and robots as well as surgeon training systems. When a surgical tool and robot is gasping a soft tissue with its gasper, the edges of the gasper can result in tissue injury because of stress concentration. It is assumed that once the peak stress is beyond a stress threshold, injury or damage can occur [2]. However, a proper in vivo stress threshold has not been reported. Damage models based on a strain energy function are necessary for very soft tissues at the moment.

5 Conclusion

A series of state-of-the-art damage models for soft tissues in animals and humans, especially those that do not consider the viscoelastic effect, was reviewed. The damage of fiber-reinforced soft tissues can be treated by using updated strain energy functions with the softening effect or by including the fiber recruitment effect with damage variables. Fiber damage can be handled by employing a strain or stretch failure criterion with a probability distribution function. Existing damage models can produce stress–stretch curves that are in very good agreement with observations, but they are less applied in healthcare, surgery, and biomedical engineering. Further, the interaction between the matrix and fibers is ignored. To develop micro-level structure damage models, microstructure visualization is necessary during the damage process in soft tissue. Crack generation and propagation in soft tissues need to be measured and simulated with suitable numerical methods. The application of damage models in biomedical engineering and clinical practice should be extended. Damage models for very soft tissues (e.g., liver, brain, kidneys, and skin) are unavailable. The viscoelasticity of soft tissue needs to be rechecked and should be considered in damage models more properly. Active stress should be taken into account in damage models for smooth muscle.