An anisotropic visco-hyperelastic constitutive model was introduced and characterized for the solid matrix of ligament and articular cartilage. The relevant numerical procedure was developed for the matrix and implemented into the commercial finite element software ABAQUS (Simulia, RI, USA) so that the constitutive model can be used for the stress analysis of general boundary-value problems.

### General formulation of the constitutive model

The deformation gradient is denoted by **F**(**X**) and the right Cauchy-Green deformation tensor by **C** = **F**
^{T}
**F**. The Helmholtz free energy function can be written as a function of the current deformation, **C**, the history of deformation, **Ξ**, and the material structure tensor, **N**
_{0}

$$ \Psi =\Psi \left(\mathbf{C},\boldsymbol{\Xi}, {\mathbf{N}}_0\right) $$

(1)

The history **Ξ** is referred to as internal variable, because it is associated with the intrinsic material properties. The structure tensor is the tensor product of the unit vector **n**
_{0} in the primary material direction, e.g. the fiber direction, **N**
_{0} = **n**
_{0} ⊗ **n**
_{0}. This tensor characterizes the anisotropy of the material. We can further introduce short and long-term internal variables to describe the history of deformation. The short-term internal variable describes the rate and history of deformation during the loading phase and can be chosen to be the time derivative of the right Cauchy-Green deformation tensor, **Ċ**. The long-term internal variable, **Γ**, determines the viscoelastic behavior of the material at larger time scales, e.g. during stress relaxation. The Helmholtz free energy function can therefore be written in terms of the short and long-term internal variables:

$$ \Psi =\Psi \left(\mathbf{C},\dot{\mathbf{C}},\boldsymbol{\Gamma}, {\mathbf{N}}_0\right) $$

(2)

It is further assumed that the elastic and viscous terms can be decoupled following previous studies [10,18,19]:

$$ \Psi \left(\mathbf{C},\dot{\mathbf{C}},\boldsymbol{\Gamma}, {\mathbf{N}}_0\right)={\Psi}^e\left(\mathbf{C},{\mathbf{N}}_0\right)+{\Psi}_s^v\left(\mathbf{C},\dot{\mathbf{C}},{\mathbf{N}}_0\right)+{\Psi}_l^v\left(\mathbf{C},\boldsymbol{\Gamma}, {\mathbf{N}}_0\right) $$

(3)

The superscripts *e* and *v* stand for elastic and viscoelastic responses respectively; and the subscripts *s* and *l* represent the short and long-term viscous responses respectively.

The viscous response is normally determined by the combined effects of several internal variables. For simplicity, however, the short and long-term internal variables were assumed to be decoupled in the present study: the viscous response during loading phase is solely determined by the short-term internal variable (\( {\varPsi}_l^v=0 \) during loading); and the stress relaxation is only determined by a set of evolutionary mechanisms (\( {\varPsi}_s^v=0 \) during relaxation).

The evolution equation commonly used for viscoelastic materials is in the differential form of viscous stress, **S**
^{v}, and elastic stress, **S**
^{e} [21-23]:

$$ {\dot{\mathbf{S}}}^v+\frac{{\mathbf{S}}^v}{\tau_i}={g}_i{\dot{\mathbf{S}}}^e $$

(4)

where *τ*
_{
i
} are time constants and *g*
_{
i
} are ratios of short-term versus equilibrium stresses (*i* = 1, 2, …). Upon solving this equation, the viscous stress can be obtained as:

$$ {\mathbf{S}}^v(t)={g}_i{\displaystyle \underset{0}{\overset{t}{\int }} \exp \left[-\left(t-T\right)/{\tau}_i\right]\ }{\dot{\mathbf{S}}}^e\ dT,\kern1em \left(0\le t\le \delta \right) $$

(5)

The time derivative of the elastic stress, **Ṡ**
^{e}, is non-zero during loading phase (0 ≤ *t* ≤ *δ*). Therefore, the viscous response depends on the rate of loading reflected in the elastic stress rate as well as the exponentially reduced relaxation function. Eq (5) is a general form of the viscous stress, so we can use it to determine the short-term viscous stress, \( {\mathbf{S}}_s^v(t)\ \left(0\le t\le \delta \right) \). When the loading phase is short, only one time constant *τ*
_{
i
} is needed, i.e. one internal variable is used for the short-term response. The long-term viscous stress (*t* > *δ*) may require a few internal variables (time constants) to account for. Each variable contributes a particular weight, *w*
_{
i
}, that decays at a specific rate associated with a time constant, *τ*
_{
i
} . The long-term viscous stress is assumed to depend on the short-term viscous stress at the end of loading phase, \( {\mathbf{S}}_{\delta}^v={\mathbf{S}}_s^v\left(t=\delta \right) \) as determined by Eq. (5). Therefore, we introduce the following form of evolution equation:

$$ {\dot{\mathbf{S}}}_l^v+\frac{{\mathbf{S}}_l^v}{\tau_i}={w}_i{\mathbf{S}}_{\delta}^v,\kern1em \left(t\ge \delta \right) $$

(6)

Here, \( {\mathbf{S}}_l^v \) can be considered as the part of the long-term viscous stress contributed by the *i*
^{th} internal variable. The long-term viscous response can then be obtained as follows:

$$ {\mathbf{S}}_l^v(t)={\displaystyle \sum_i{w}_i}{\mathbf{S}}_{\delta}^v{\displaystyle \underset{\delta }{\overset{t}{\int }} \exp \left[-\left(t-T\right)/{\tau}_i\right]\ }dT,\kern1em \left(t\ge \delta \right), $$

(7)

This equation shows the decay of the long-term viscous stress (*t* > *δ*) from the peak stress at the end of the loading phase, \( {\mathbf{S}}_{\delta}^v \). It should be noted that *g*
_{
i
} in Eq. (4) represents the magnitude of stress relaxation. In the standard linear solid model of viscoelasticity, it can be interpreted as the ratio of the stiffness of the Maxwell body to the stiffness of the elastic body. On the other hand, a *w*
_{
i
} reflects the contribution associated with an internal variable to the total response (dimension is 1/time as shown in Eq. (6), *Σw*
_{
i
} = 1).

In consistence with Eq. (3), the total stress response of the material can be written as the summation of the elastic stress (**S**
^{e}), short \( \left({\mathbf{S}}_s^v\right) \) and long-term \( \left({\mathbf{S}}_l^v\right) \) viscous stresses using the same index convention:

$$ \mathbf{S}(t)={\mathbf{S}}^e(t)+{\mathbf{S}}_s^v(t)+{\mathbf{S}}_l^v(t) $$

(8)

The long-term response is naturally contributed from multiple mechanisms with different time constants *τ*
_{
i
} and weights *w*
_{
i
}. In the examples to follow, the time constants *τ*
_{
i
} were considered to be independent of strain. Therefore, a quasi-linear form of viscoelasticity was obtained. However, the time constant can also be defined as a function of strain leading to a fully nonlinear description of viscoelasticity [24].

### Particular formulation for ligaments

The extracellular matrix of ligament is mainly composed of a dense network of collagen fibers mostly aligned in the longitudinal direction of the tissue. This structure forms a transversely isotropic characteristic to the tissue. In the present study, ligament was modeled as a composite material consisting of an isotropic non-fibrillar matrix reinforced by the fiber network. As a ligament is physiologically subjected to tension when the fibers sustain most of the load, only the fibrillar matrix was considered viscoelastic and the non-fibrillar matrix was modeled as elastic. The energy of the non-fibrillar matrix, as indicated by the subscript *m*, is uniquely determined by the deformation, **C**. The energy of the fibrillar matrix, as indicated by the subscript *f*, is also associated with the structure tensor **N**
_{0} and the internal variables. Therefore, the Helmholtz free energy for ligaments can be written in the decoupled form as

$$ \Psi \left(\mathbf{C},\dot{\mathbf{C}},\boldsymbol{\Gamma}, {\mathbf{N}}_0\right)={\Psi}_m^e\left(\mathbf{C}\right)+{\Psi}_f^e\left(\mathbf{C},{\mathbf{N}}_0\right)+{\Psi}_{s,f}^v\left(\mathbf{C},\dot{\mathbf{C}},{\mathbf{N}}_0\right)+{\Psi}_{l,f}^v\left(\mathbf{C},\Gamma, {\mathbf{N}}_0\right) $$

(9)

The non-fibrillar matrix was considered to be Neo-Hookean (first term), while an exponential function was adopted for the hyperelastic part of the fibrillar matrix (second term) [19]. The exponential function describes the stiffening of the tissue at larger strains arisen by the gradual recruitment of collagen fibers during further deformation [25]. Although the collagen fibers show a great stiffness in tension, they cannot sustain compression due to their slenderness. Therefore, a piecewise energy function was employed for the second term in equation (9) that determines whether the collagen fibers contribute to the load support based on their deformation state:

$$ \begin{array}{l}{\Psi}^e={\Psi}_m^e+{\Psi}_f^e\\ {}\kern1em =\left\{\begin{array}{l}{a}_1\left({I}_1-3\right),\kern10em \mathrm{if}\kern0.75em {I}_4\le 1\kern0.5em \left(\mathrm{compression}\right)\hfill \\ {}\hfill {a}_1\left({I}_1-3\right)+\frac{a_2}{2{a}_3} \exp \left[{a}_3{\left({I}_4-1\right)}^2\right],\kern0.75em \mathrm{if}\kern0.75em {I}_4>1\kern0.5em \left(\mathrm{tension}\right)\hfill \end{array}\right.\end{array} $$

(10)

Here, *I*
_{1} is the first invariant of **C** defined as *I*
_{1} = **C** : **I** = *C*
_{
ij
}
*I*
_{
ij
} = *tr*{**C**}, and *I*
_{4} is an invariant defined as *I*
_{4} = **C** : **N**
_{0}. Therefore, *I*
_{4} is actually the square of the stretch ratio *λ*
_{
f
} in the fiber direction, i.e. \( {I}_4={\lambda}_f^2 \). As can be seen in this equation, the contribution of the collagen fibers is zero when the tissue is under compression, i.e., \( {\varPsi}_f^e=0 \) when *I*
_{4} ≤ 1. The material constants, *a*
_{1}, *a*
_{2} and *a*
_{3}, must be positive to ensure the convexity of the function.

The short term viscous function was proposed as follows [19]:

$$ {\Psi}_{s,f}^v=\left\{\begin{array}{l}0\\ {}0.5\kern0.5em {a}_4\left({I}_4-1\right) \exp \left[{a}_5{\left({I}_4-1\right)}^2\right]{J}_5, \end{array}\right.\kern1em \begin{array}{c}\hfill \mathrm{if}\kern0.75em {I}_4\le 1\hfill \\ {}\hfill \mathrm{if}\kern0.75em {I}_4>1\hfill \end{array} $$

(11)

where *J*
_{5} is an invariant defined as *J*
_{5} = **Ċ**
^{2} : **N**
_{0} and *a*
_{4} and *a*
_{5} are viscous material parameters. As mentioned previously, only collagen fibers were considered viscoelastic and thus, the anisotropic energy functions are written only for the fiber network. Similar to the elastic part, the functions are zero under compression (*I*
_{4} ≤ 1). This form of short-term viscous function includes the nonlinearities of the short-term viscous response associated with strain (*I*
_{4}) by the nonlinear exponential function of *I*
_{4}.

The second Piola-Kirchhoff stress for elastic and viscous terms are then derived from the energy functions [9,19]:

$$ {\mathbf{S}}_m^e=2{a}_1\mathbf{I} $$

(12)

$$ {\mathbf{S}}_f^e=2{a}_2 \exp \left[{a}_3{\left({I}_4-1\right)}^2\right]\left({I}_4-1\right){\mathbf{N}}_0 $$

(13)

$$ {\mathbf{S}}_s^v={a}_4\left({I}_4-1\right) \exp \left[{a}_5{\left({I}_4-1\right)}^2\right]\left({\mathbf{N}}_0\dot{\mathbf{C}}+\dot{\mathbf{C}}{\mathbf{N}}_0\right) $$

(14)

The total stress can then be written as the summation of the stresses obtained hitherto (Eqs. 12–14) as well as the long-term stress (Eq. 7) as a piecewise function:

$$ \mathbf{S}(t)=\left\{\begin{array}{l}\begin{array}{l}p{\mathbf{C}}^{-1}+2{a}_1\mathbf{I}, \kern14.81em \mathrm{if}\kern0.75em {I}_4\le 1\\ {}\ \end{array}\hfill \\ {}\hfill \begin{array}{l}p{\mathbf{C}}^{-1}+2{a}_1\mathbf{I}+2{a}_2 \exp \left[{a}_3{\left({I}_4-1\right)}^2\right]\left({I}_4-1\right){\mathbf{N}}_0+{a}_4\left({I}_4-1\right) \exp \left[{a}_5{\left({I}_4-1\right)}^2\right]\left({\mathbf{N}}_0\dot{\mathbf{C}}+\dot{\mathbf{C}}{\mathbf{N}}_0\right)\\ {}+{\displaystyle \sum_i{w}_i}{\displaystyle {\int}_{\delta}^t \exp \left[-\left(t-T\right)/{\tau}_i\right]{\mathbf{S}}_s^vdT, \kern6.5em \mathrm{if}\kern0.75em {I}_4>1}\end{array}\hfill \end{array}\right. $$

(15)

where *p* is a Lagrange multiplier that was introduced to enforce the tissue incompressibility.

The material parameters can be obtained from uniaxial tensile tests. Assuming the *x* direction to be the direction of loading and considering the material to be incompressible, the deformation gradient is of the following form:

$$ \mathbf{F}={\left(\lambda,\ \frac{1}{\sqrt{\lambda }},\ \frac{1}{\sqrt{\lambda }}\right)}^T\mathbf{I} $$

(16)

from which the right Cauchy-Green deformation tensor can be obtained:

$$ \mathbf{C}={\left({\lambda}^2,\ \frac{1}{\lambda },\ \frac{1}{\lambda}\right)}^T\mathbf{I} $$

(17)

Using the aforementioned matrix (Eq. 17) along with the structure tensor obtained from the unit vector of fibers \( {n}_0=\overrightarrow{i} \) to obtain *I*
_{4}, the components of the hyperelastic second Piola-Kirchhoff stress will be obtained from Eqs. (12–14):

$$ {S}_{11}^e=2{a}_1+p/{\lambda}^2+2{a}_2\left({\lambda}^2-1\right) \exp \left[{a}_3{\left({\lambda}^2-1\right)}^2\right] $$

(18)

$$ {S}_{22}^e={S}_{33}^e=2{a}_1+p\lambda $$

(19)

As there is no stress in the lateral directions, i.e., \( {S}_{22}^e={S}_{33}^e=0 \), the Lagrange multiplier, *p*, can be calculated as:

$$ p=-2{a}_1/\lambda $$

(20)

The time derivative of the right Cauchy-Green deformation tensor, **Ċ**, is also needed for calculation of the viscous stresses:

$$ \dot{\mathbf{C}}={\left(2\lambda \dot{\lambda},-\dot{\lambda}/{\lambda}^2,-\dot{\lambda}/{\lambda}^2\right)}^T\mathbf{I} $$

(21)

The viscous stress can also be obtained from Eq. (14) similarly by substituting the corresponding matrices (Eqs. 17 and 21):

$$ {S}_{11}^v=4{a}_4\lambda \left({\lambda}^2-1\right) \exp \left({a}_5{\left({\lambda}^2-1\right)}^2\right)\dot{\lambda} $$

(22)

The first Piola-Kirchhoff or nominal stress is commonly used in experimental studies. It can be obtained from the second Piola-Kirchhoff stress according to:

$$ \mathbf{P}={J}^{-1}\ \mathbf{F}\ \mathbf{S} $$

(23)

Therefore, the nominal stress in the direction of loading was derived for data fit as follows:

$$ {P}_{11}=\left\{\begin{array}{l}2{a}_1\left(\lambda -1/{\lambda}^2\right), \kern14em \mathrm{if}\kern0.75em {I}_4\le 1\hfill \\ {}\hfill \begin{array}{l}2{a}_1\left(\lambda -1/{\lambda}^2\right)+2{a}_2\lambda \left({\lambda}^2-1\right) \exp \left[{a}_3{\left({\lambda}^2-1\right)}^2\right]+4{a}_4\left({\lambda}^2-1\right){\lambda}^2 \exp \left[{a}_5{\left({\lambda}^2-1\right)}^2\right]\dot{\lambda}\\ {}+{\displaystyle \sum_{i=1}^N{w}_i{\displaystyle {\int}_{\delta}^t \exp \left[-\left(t-T\right)/{\tau}_i\right]{P}_{11}^vdT, \kern5.5em \mathrm{if}\kern0.75em {I}_4>1}}\end{array}\hfill \end{array}\right. $$

(24)

### Particular formulation for articular cartilage

Articular cartilage can be modeled as a fluid-saturated non-fibrillar matrix reinforced by a collagen network [26,27] with various fiber alignments depending on the depth of the tissue and site where the tissue is located in the joint. A thin and narrow strip of cartilage aligned in the fiber direction, as used in uniaxial tensile testing, is similar to a ligament as the fluid pressurization does not affect the tensile properties significantly [28]. The hyperelastic Helmholtz free energy function for cartilage is also composed of isotropic and anisotropic parts representing the non-fibrillar and fibrillar matrices respectively. Similar to the ligament model (Eq. 10), a piece wise function is used to describe the tension-compression difference of the fibers in articular cartilage:

$$ {\Psi}^e=\left\{\begin{array}{l}{b}_1\left({I}_1-3\right), \\ {}{b}_1\left({I}_1-3\right)+\frac{1}{2}{b}_2{\left({I}_4-1\right)}^2+\frac{1}{3}{b}_3{\left({I}_4-1\right)}^3, \kern0.5em \end{array}\right.\begin{array}{c}\hfill \mathrm{if}\kern0.75em {I}_4\le 1\hfill \\ {}\hfill \mathrm{if}\kern0.75em {I}_4>1\hfill \end{array} $$

(25)

The short-term viscous potential was introduced to resemble Equation (11) as follows:

$$ {\Psi}_s^v=\left\{\begin{array}{c}\hfill 0,\kern6.25em \mathrm{if}\kern0.5em {I}_4\le 0\hfill \\ {}\hfill {b}_4\left({I}_4-1\right){J}_5 \ln {I}_4,\kern0.5em \mathrm{if}\kern0.5em {I}_4>0\hfill \end{array}\right. $$

(26)

By enforcing incompressibility for the tissue in tensile testing, the second Piola-Kirchhoff stress is obtained:

$$ \mathbf{S}(t)=\left\{\begin{array}{l}p{\mathbf{C}}^{-1}+2{b}_1\mathbf{I}, \kern12.4em \mathrm{if}\kern0.75em {I}_4\le 1\hfill \\ {}\hfill \begin{array}{l}p{\mathbf{C}}^{-1}+2{b}_1\mathbf{I}+2{b}_2\left({I}_4-1\right){\mathbf{N}}_0+2{b}_3{\left({I}_4-1\right)}^2{\mathbf{N}}_0+2{b}_4 \ln {I}_4\left({I}_4-1\right)\left({\mathbf{N}}_0\dot{\mathbf{C}}+\dot{\mathbf{C}}{\mathbf{N}}_0\right)\\ {}+{\displaystyle \sum_i{w}_i}{\displaystyle {\int}_{\delta}^t \exp \left[-\left(t-T\right)/{\tau}_i\right]{\mathbf{S}}_s^vdT, \kern4em \mathrm{if}\kern0.75em {I}_4>1}\end{array}\hfill \end{array}\right. $$

(27)

This equation is similar to Eq. (15), and *p* is the Lagrange multiplier.

Similar to Eq. (24), the nominal stress in a uniaxial tensile test can be obtained:

$$ {P}_{11}=\left\{\begin{array}{l}2{b}_1\left(\lambda -1/{\lambda}^2\right), \kern11.4em \mathrm{if}\kern0.75em {I}_4\le 1\hfill \\ {}\hfill \begin{array}{l}2{b}_1\left(\lambda -1/{\lambda}^2\right)+2{b}_2\lambda \left({\lambda}^2-1\right)+2{b}_3\lambda {\left({\lambda}^2-1\right)}^2+8{b}_4{\lambda}^2 \ln {\lambda}^2\left({\lambda}^2-1\right){\lambda}^2\\ {}+{\displaystyle \sum_i{w}_i}{\displaystyle {\int}_{\delta}^t \exp \left[-\left(t-T\right)/{\tau}_i\right]{P}_{11}^vdT, \kern2.75em \mathrm{if}\kern0.75em {I}_4>1}\end{array}\hfill \end{array}\right. $$

(28)

This 1D formulation can then be extended to a 3D formulation of the solid matrix for articular cartilage. Darcy’s law must be incorporated to account for fluid pressure in the tissue under compressive loadings.

### Numerical implementation

The constitutive model proposed here was implemented in the finite element software ABAQUS (Simulia, Providence, RI) by means of a UMAT, a user-defined subroutine in FORTRAN. The constitutive behavior used in ABAQUS for the fluid flow is governed by Forchheimer’s law that includes Darcy’s law as a linear case when the fluid velocities are low. The Newton–Raphson method is used in ABAQUS which requires the stress and the stiffness matrices to be updated in each time step. The Jaumann rate of stress is used in calculations and correspondingly, the Jaumann elasticity tensor along with the Cauchy stress are needed. The Cauchy stress tensor can be obtained by pushing the second Piola-Kirchhoff stress forward [29]:

$$ \boldsymbol{\upsigma} ={J}^{-1}\mathbf{F}\ \mathbf{S}\ {\mathbf{F}}^T $$

(29)

The elasticity and viscosity tensors can also be derived from the free energy functions or stresses:

$$ {\operatorname{C}}^e=2\frac{\partial {\mathbf{S}}^e}{\partial \mathbf{C}}=4\frac{\partial^2{\Psi}^e}{\partial {\mathbf{C}}^2} $$

(30)

$$ {\operatorname{C}}^v=2\frac{\partial {\mathbf{S}}^v}{\partial \dot{\mathbf{C}}}=4\frac{\partial^2{\Psi}^v}{\partial {\dot{\mathbf{C}}}^2} $$

(31)

The elasticity and viscosity tensors were obtained and used to arrive at the Jaumann elasticity tensor following the procedure in [30]. The Jacobian matrix was derived analytically, resulting in a quadratic convergence rate instead of a slower convergence rate associated with the numerical approximation of the Jacobian matrix [31,32].

Reinforcing the incompressibility condition strictly in a numerical procedure can cause problems such as mesh locking. Therefore, a slight compressibility should be introduced into the model which requires the deviatoric-volumetric decomposition of the Helmholtz free energy [33]. For nearly incompressible materials, the decomposition can be applied on both isotropic and anisotropic components of the energy function.

$$ \Psi \left({I}_1,{I}_3,{I}_4,{J}_5\right)={\Psi}_{vol}(J)+{\Psi}_m\left({\overline{I}}_1\right)+{\Psi}_f\left({\overline{I}}_4,{\overline{J}}_5\right) $$

(32)

Here, *Ψ*
_{
m
}(*Ī*
_{1}) is the deviatoric term for the isotropic matrix, and \( {\varPsi}_f\left({\overline{I}}_4,{\overline{J}}_5\right) \) for the anisotropic collagen network. The invariants with bar are the invariants of the deviatoric part of the right Cauchy-Green deformation tensor and its time derivative, which are associated with viscoelastic behavior [34].

For the purpose of comparison, the QLV solution below was also numerically implemented following a previously established procedure [27,35]

$$ \mathbf{S}(t)={\mathbf{S}}^e(t)+{\displaystyle \sum_i{g}_i}{\displaystyle {\int}_0^t \exp \left[-\left(t-T\right)/{\tau}_i\right]{\dot{\mathbf{S}}}^edT} $$

(33)

The data fitting was performed with the least square method using the optimization tool box in MATLAB (MathWorks, MA). Constraints on the parameters were applied whenever applicable. For example, the material properties must be positive to ensure the convexity of the energy functions and positive-definiteness of the elasticity tensors.