## Abstract

Living biological tissues are complex structures that have the capacity of evolving in response to external loads and environmental stimuli. The adequate modelling of soft biological tissue behaviour is a key issue in successfully reproducing biomechanical problems through computational analysis.

This study presents a general constitutive formulation capable of representing the behaviour of these tissues through finite element simulation. It is based on phenomenological models that, used in combination with the generalized mixing theory, can numerically reproduce a wide range of material behaviours.

First, the passive behaviour of tissues is characterized by means of hyperelastic and finite-strain damage models. A generalized damage model is proposed, providing a flexible and versatile formulation that can reproduce a wide range of tissue behaviour. It can be particularized to any hyperelastic model and requires identifying only two material parameters. Then, the use of these constitutive models with generalized mixing theory in a finite strain framework is described and tools to account for the anisotropic behaviour of tissues are put forth.

The active behaviour of tissues is characterized through constitutive models capable of reproducing the growth and remodelling phenomena. These are built on the hyperelastic and damage formulations described above and, thus, represent the active extension of the passive tissue behaviour. A growth model considering biological availability is used and extended to include directional growth. In addition, a novel constitutive model for homeostatic-driven turnover remodelling is presented and discussed. This model captures the stiffness recovery that occurs in healing tissues, understood as a recovery or reversal of damage in the tissue, which is driven by both mechanical and biochemical stimuli.

Finally, the issue of correctly identifying the material parameters for computational modelling is addressed. An inverse method using optimization techniques is developed to facilitate the identification of these parameters.

# 1 Introduction

## 1.1 Motivation

Computational biomechanics is an emerging field that embraces a broad range of computational modelling techniques used in the numerical simulation and analysis of biological systems. Research in biomechanics, which seeks to understand the mechanics of living systems [3], is generally aimed at improving our knowledge of the human body to advance in medical science and technology. Unsurprisingly, biomechanics plays an important role in what is now termed in silico medicine [4]. This new discipline aims at capturing aspects of the physiology and pathology of the human body in computer models that aid in the clinical prevention, diagnosis and treatment of injury and disease. In this context, finite element analysis (FEA) [5] is proving to be a powerful tool to perform certain patient-specific biomechanical studies and, in this way, provide additional data for clinical decisions.

Computational biomechanics and in silico medicine overlap in areas such as orthopaedics, rehabilitation, gait analysis, tissue engineering, hemodynamics and mechanobiology, to name but a few. The latter is of particular interest for researchers in continuum mechanics and structural analysis since it poses the challenge of extending well-established constitutive formulations to be able to model biological tissue behaviour. Constitutive models describe the macroscopic behaviour resulting from the internal constitution of a material in order to characterize its response to external stimuli [6].

The constitutive modelling of biological tissues presents the difficulty of having to account for the fact that these materials are alive, so they respond and adapt to external and internal actions of both mechanic and metabolic origins. Furthermore, computational modelling of living tissues tends to involve complex and evolving geometries, large displacements and strains, and often their loads, boundary conditions and interactions are hard to quantify and establish with accuracy [7]. Further challenges include accounting for their hierarchical structure and the multiple biophysical stimuli at the different spatial and time scales involved.

## 1.2 Overview

Biological tissues are constituted by different components arranged in a hierarchical structure whose properties strongly depend on the size, distribution and geometry of these constituents. Biomaterials are typically classified into either hard or soft tissues. The former include mineralized tissues with a rigid intercellular substance, e.g., bone and teeth. They suffer small deformations and behave nearly elastically in the physiological range. The latter are tissues composed of an extracellular matrix (ECM) of collagen and elastin fibres embedded in ground substance. Examples of soft tissues are tendons, ligaments, skin, muscles and blood vessels, amongst others. They are highly deformable and exhibit a nonlinear elastic behaviour. In addition, the ground substance, which is basically a hydrophilic gel, confers near-incompressibility to these tissues.

The properties of both hard and soft tissues are classified into two separate groups, passive and active, according to their predominantly mechanical or biological nature. Properties that are not directly determined by the biochemical and biophysical processes taking place in the tissue are known as passive properties. In addition to the nonlinearity and near-incompressibility mentioned above, soft tissues may exhibit passive properties such as anisotropy, residual stresses, viscoelasticity, plasticity, damage and failure. In contrast, active properties are directly dependent on the metabolism, which keeps the tissues alive and allows them to adapt and evolve in response to their environment. Growth, atrophy, remodelling, healing, regeneration and ageing are all active properties.

Constitutive models to represent these properties can be developed following either a mechanistic or a phenomenological approach. The mechanistic approach tries to understand and model the mechanisms that regulate tissue reactions at cellular level. By incorporating the biophysics behind the tissue behaviour, these mathematical models allow testing different hypothesis and, ultimately, provide a better understanding of the mechanobiological interactions involved in the processes being modelled. However, they require in-depth knowledge of the tissue histology and the biochemical reactions taking place at cellular and molecular levels, which in some cases are not completely understood yet.

Phenomenological models, on the other hand, establish direct relations between the external stimuli and the observed tissue response without trying to explain the mechanisms behind the observations. They describe tissue behaviour using functional relations that closely fit experimental studies. The phenomenological approach avoids quantifying microscopic quantities. Instead, the internal variables of the model are directly associated with the behaviour observed at macroscopic level. These variables are expressed in continuous terms even though they are indirectly related to the mechanisms taking place at cellular and molecular levels. Nonetheless, biological tissues are extremely heterogeneous and, often, cannot be modelled with precision by using this approach alone.

The characterization of material inhomogeneity due to the local variation of the tissue structure is sometimes addressed with microstructural models. Histological data is used to describe the tissue architecture and then, by means of homogenization techniques, the macroscopic behaviour of the biomaterial is obtained. The underlying microstructural relations are mathematically complex and accurate quantification of the interactions between the constituents at the different organizational scales is required. Humphrey and Yin [8] proposed using mixing theory as an alternative, combining in this way desirable features from the phenomenological and microstructural approaches.

Mixing theory idealizes the composite material (the tissue) as composed of several individual compounds (fibres and matrix). Each of these components or simple materials is modelled, in turn, by the constitutive model of choice. Phenomenological models are generally used to describe the behaviour of the simple material, although mechanistic, microstructural or, even, mixing theory could be used for this purpose. The use of phenomenological models provides a compromise between the physical accuracy of the microstructural approach and the mathematical simplicity of the purely phenomenological one.

## 1.3 Goals

The aim of this study is to set the bases for a general constitutive formulation that represents the behaviour of soft biological tissues through numerical simulation, specifically, the finite element method (FEM). This formulation is based on phenomenological models used in combination with the generalized mixing theory [9] and implemented in the in-house code PLCd [1].

PLCd is an implicit FE code developed in Fortran and capable of solving finite-strain nonlinear three-dimensional solid mechanics problems. It uses the direct sparse solver Pardiso [10] and a full Newton algorithm. Originally developed for the modelling of fracture, dynamics behaviour and analysis of composite structures through mixing theory, its scope has extended over the years to topics such as multi-scale homogenization, fatigue analysis, the study of masonry structures and, now, biological tissue simulation. Among its strengths are the thorough implementation of the generalized mixing theory and the straightforward organization and accessibility of the code, which facilitates the introduction of new constitutive models. Nevertheless, since PLCd is research-oriented and in constant development, the implementation of new constitutive formulations often requires additional changes and improvements in the general structure of the code.

The mixing theory is understood as a constitutive equation manager that allows mixing and matching simple materials, whose behaviour is modelled by means of phenomenological constitutive formulations, to obtain the overall behaviour of the composite material. In this study, biological tissues are assimilated to a composite material and, therefore, the appropriate constitutive models must be formulated to represent its simple components. Hence, the idea of a general constitutive formulation for biological tissues: the use of mixing theory in conjunction with the simple constitutive models confers flexibility and versatility and, ultimately, generality to the overall formulation.

Exploiting the aforementioned strong points of the formulation, phenomenological constitutive models capable of reproducing biological tissue behaviour are developed and implemented in order to achieve the goal of this study, a general constitutive formulation. This task is tackled in two distinct parts: the passive and active properties.

First, the passive behaviour of biological tissues is addressed by means of quasi-incompressible hyperelastic and finite-strain damage models. This is, of course, a simplification since tissues are known to exhibit more complex behaviours such as viscoelasticity. However, the aim of the study is to obtain a general formulation and the hyperelastic and damage models are deemed enough for the purpose of representing the basic passive properties of the tissues.

Then, the active behaviour of the tissues is modelled through continuum growth and healing models, which are based on the previous hyperelastic and finite-strain damage formulations. Again, the active properties of living tissues include much more than the growth and remodelling modelled in this study. Yet, the scope of the general constitutive formulation is limited to these under the assumption that they are the most notable of the active properties observed in adult functional tissues.

The use of phenomenological models requires the correct identification of its material parameters to obtain an adequate computational representation of the real tissue behaviour. Due to the non-mechanistic nature of these parameters, it is often difficult to establish their values based on physical or structural information. Identifying the correct material parameters tends to require a trial-and-error approach, which is often tedious and generally time-consuming. Parameter identification through inverse methods by means of optimization techniques is an alternative to the manual identification. Such a method is developed, tailored to the particular applications of this study, but flexible enough to be easily modifiable if required for other purposes.

In view of the above, the aim of this study is divided into the following objectives, each of which corresponds to a chapter in this monography:

1. Development of a set of constitutive formulations capable of representing the basic passive properties of biological tissues.
2. Development of a set of constitutive formulations capable of representing the basic active properties of biological tissues.
3. Identification of the material parameters of the constitutive formulations through an inverse method and adequate optimization techniques.

These constitutive formulations are developed in the framework of finite-strain continuum solid mechanics and are built on a well-established thermodynamic basis. The computational implementation and validation of the models require numerical improvements in areas of the code such as FE formulation, integration scheme, loading configurations and pre- and post-processor interfaces.

## 1.4 Outline

This monography is organized as follows:

Chapter 2 A set of constitutive equations capable of reproducing the basic passive properties of biological tissues is proposed. First, several hyperelastic formulations are reviewed and the implementation of two of these models is described in detail. The finite strain framework in which these formulations are developed and the numerical requirements imposed by the quasi-incompressibility assumption are also described. Then, a generalized finite-strain damage model is developed and implemented based on these hyperelastic models. A few examples serve to illustrate its main characteristics and validate the formulation. The use of these constitutive models together with generalized mixing theory is studied. Improvements to better represent the passive properties of biological tissues are proposed and discussed.

Chapter 3 A set of constitutive equations capable of reproducing the main active properties of biological tissues is proposed. Continuum modelling of growth and remodelling is reviewed and a constitutive model for volumetric growth that takes into account the metabolic contribution is described and implemented. An extension of this model to include directional growth is also presented. Then, a novel reverse-damage model to account for healing in soft tissues is formulated, implemented and validated through numerical examples. Similarly to the growth formulation, healing is driven by mechanical stimuli but subject to biological availability, as allowed by the tissue's metabolism. The integration of the growth and healing models in a common formulation in order to better represent the active properties of biological tissues is proposed and discussed.

Chapter 4 The constitutive formulation described in Chapter 2 is used to model the basic components of the human cervical spine and an attempt to numerically reproduce a discectomy spine surgery is presented. The difficulties encountered in the modelling of the spine are discussed and a solution is proposed to solve the material parameter identification problem. An inverse method using optimization techniques, tailored to this particular application, is developed. Then, the method is modified and used in the material parameter identification of a generic composite specimen, demonstrating the versatility of the method proposed. Finally, the improvements required in the cervical spine modelling are discussed, in addition to the capabilities of the proposed material parameter identification method.

Chapter 5 The achievements of this study are summed up, final conclusions are drawn and future work lines are outlined.

Note that a chapter devoted to reviewing the state of the art is not included at the beginning of the monography. The self-contained nature of each chapter lends itself to including specific literature reviews in each chapter. Likewise, conclusions related to the specific content of each chapter are pointed out at the end of the same.

## 1.5 Research dissemination

The work included in this monography resulted in the following scientific publications:

## Chapter 2

1. E. Comellas, F.J. Bellomo and S. Oller. A generalized finite-strain damage model for quasi-incompressible hyperelasticity using hybrid formulation. International Journal for Numerical Methods in Engineering, 2015. doi:10.1002/nme.5118.

## Chapter 3

1. F.J. Bellomo, E. Comellas, L. Nallim and S. Oller. Numerical simulation of the directioned growth and remodelling of soft biological tissues generated by mechanical stimuli (in Spanish). In: Mecánica Computacional Vol. XXXIII. Number 41. Computational Modeling in Bioengineering Applications (A), pages 2667–2676. Asociación Argentina de Mecánica Computacional, 2014.
2. E. Comellas, T.C. Gasser, F.J. Bellomo and S. Oller. A homeostatic-driven turnover remodelling constitutive model for healing in soft tissues. Submitted to the Journal of the Royal Society Interface. October 2015.

## Chapter 4

1. E. Comellas, S. Oller, J. Poblete, J. Berenguer and A. Prats-Galino. Numerical modelling of a cervical spine discectomy. In: Mecánica Computacional Vol. XXXI. Number 24. Computational Modeling in Bioengineering and Biomedical Systems (A), pages 3811-3826. Asociación Argentina de Mecánica Computacional, 2012.
2. E. Comellas, S.I. Valdez, S. Oller and S. Botello. Optimization method for the determination of material parameters in damaged composite structures. Composite Structures, 2015, 122:417–424. doi:10.1016/j.compstruct.2014.12.014.

In addition, part of the work was presented at the following conferences and workshops:

## Unpublished Conference Presentations

1. E. Comellas, T.C. Gasser, F.J. Bellomo and S. Oller. A reverse damage model for healing in soft biological tissues, VI International Conference on Computational Bioengineering (ICCB 2015), Barcelona, September 14 – 16, 2015.
2. E. Comellas, F.J. Bellomo and S. Oller. Ogden parameter optimization for finite element modelling of cervical ligaments, invited to present in Biomechanics Session of XXXII Annual Conference of the Spanish Society of Biomedical Engineering (CASEIB 2014), Barcelona, November 26 – 28, 2014.
3. E. Comellas, F.J. Bellomo and S. Oller. Ogden parameter optimization for finite element modelling of cervical ligaments using hybrid formulation, 11${\displaystyle ^{th}}$ World Congress on Computational Mechanics (WCCM XI), Barcelona, June 20 – 25, 2014.

## Unpublished Workshop Presentations

1. E. Comellas and S. Oller. Non-linear constitutive modelling of biological tissues, Course on Advanced Structural Analysis Interface, MuMoLaDe Michaelmas School a Marie Curie ITN, Barcelona, October 1 – 2, 2015.
2. E. Comellas, T.C. Gasser, F.J. Bellomo and S. Oller. A damage-driven model for remodelling/healing of soft tissues considering biological availability, International Workshop on Modelling across the Biology-Mechanics Interface, Castro Urdiales, September 1 – 4, 2015.
3. E. Comellas, F.J. Bellomo and S. Oller. Ogden parameter optimization for finite element modelling of cervical ligaments. IV Meeting of the Spanish Chapter of the European Society of Biomechanics (ESB), Valencia, November 20 – 21, 2014. Best valued works, 2${\displaystyle ^{nd}}$ position.
4. E. Comellas. TCAiNMaND Progress in WP5. TCAiNMaND First Annual Workshop, Barcelona, July 18, 2014.
5. E. Comellas, F.J. Bellomo and S. Oller. FE modeling of cervical ligaments using Ogden hyperelasticity and hybrid formulation, III Meeting of the Spanish Chapter of the European Society of Biomechanics (ESB), Barcelona, October 23 – 24, 2013.
6. E. Comellas, S. Oller, J. Poblete, J Berenguer and A. Prats-Galino. Numerical simulation of a cervical discectomy surgery (in Spanish), II Meeting of the Spanish Chapter of the European Society of Biomechanics (ESB), Sevilla, October 25, 2012.

Finally, part of the work presented in this document is the result of collaborating with external researchers during the following research stays:

1. KTH Royal Institute of Technology, 6-month doctoral research stay. Worked under the direct supervision of Prof. Christian Gasser in the Solid Mechanics Department of the KTH Royal Institute of Technology. Stockholm, Sweden. February – July 2015. The article “A homeostatic-driven turnover remodelling constitutive model for healing in soft tissues” is the result of the work developed during the stay.
2. CIMAT Center for Research in Mathematics, 2-month research stay in the framework of the TCAiNMaND project, a Marie Curie International Research Staff Exchange Scheme (IRSES) under the European Union 7${\displaystyle ^{th}}$ Framework Programme with grant agreement n. 612607. Worked under the direct supervision of Dr. Salvador Botello in the Computational Sciences Department of the Center for Research in Mathematics (CIMAT). Guanajuato, México. February – March 2014. The article “Optimization method for the determination of material parameters in damaged composite structures” is the result of the work developed during the stay.

# 2 Constitutive modelling of passive properties

## 2.1 Background

The adequate modelling of a biological tissue's constitutive behaviour is key in successfully reproducing biomechanical problems through computational analysis. The nonlinearity exhibited by soft tissue in response to loading was studied and experimentally quantified as early as mid-19th century [11]. The pioneering researchers in biomechanics observed that this behaviour is similar to that of rubber-like materials due to the long-chain, cross-linked polymeric structures of both types of materials. Hence, the introduction of continuum finite elasticity theory, used in the modelling elastomers, to model soft tissue behaviour.

Continuum-based theories deal with matter at a macroscopic scale, i.e., the behaviour of materials is studied considering matter as a continuous medium rather than as an heterogeneous mass formed by discrete particles. Also, the assumption is made that the solid and its properties are describable by continuous functions which have continuous derivatives. The existence of a continuous displacement field in this type of problems makes the implementation of finite element techniques to solve them much simpler.

Fung, considered by many the father of modern biomechanics, was among the first to characterize soft tissue biomechanics within the framework of finite-strain elasticity [12,3]. So, these first efforts to numerically reproduce soft tissue behaviour were mainly focused at determining the adequate elastic potentials or strain energy functions ${\textstyle \Psi }$ for such purpose.

Most of the models employed in representing soft tissue behaviour were borrowed or adapted, and still are, from elastomer applications. These models assume the material is homogeneous, nearly-incompressible and isotropic, which can be over-simplifying hypotheses for some biological tissues. Yet, these assumptions allow defining the elastic potential in terms of either the first and second invariant values of the right Cauchy-Green deformation tensor ${\textstyle \mathbf {C} }$, or the principal stretches, which are the square roots of the eigenvalues of ${\textstyle \mathbf {C} }$. This is especially useful in the mathematical development and numerical implementation of the hyperelastic formulation, as will be seen in section 2.2.

The main challenge is then to choose an appropriate set of invariants such that the material model includes as few parameters as possible but is still able to adequately reproduce complex deformation states of the material it is representing. In addition, the material parameters should be determined by fitting solely a small number of simple experiments. Among the most popular models of these type are the Mooney-Rivlin type, based on the first and second invariants of ${\textstyle \mathbf {C} }$, and the Ogden models, based on the principal stretches. Another commonly used polynomial form of strain energy function is the one introduced by Yeoh, which uses power terms of the first invariant of ${\textstyle \mathbf {C} }$. Exponential and logarithmic forms are also widely used, for example in the Fung and Veronda-Westmann models. The reader is referred to the comprehensive review on hyperelastic models for rubber-like materials by Steinmann et al. [13] for detailed information on many other constitutive models not mentioned in this study. Chagnon et al. [14] provide a comprehensive review of strain energy density functions for soft biological tissues. Studies of hyperelastic models for particular applications are also available, e.g. Auricchio et al. [15] review hyperelastic models for human aortic valve tissue.

The numerical and thermodynamic bases of isotropic quasi-incompressible hyperelasticity will be reviewed in section 2.2. The hyperelastic models typically used in biomechanical applications will be briefly described and discussed and, then, the complete neo-Hookean and Ogden formulations will be developed. Their numerical implementation in the in-house FE code PLCd [1] will be described and validated through numerical examples.

Although the hyperelastic models used to reproduce soft-tissue behaviour described above are able to adequately reproduce the behaviour of biological tissues in some applications, the anisotropy exhibited in most cases cannot be ignored. In order to address this issue, several authors [16,17,18] developed elastic potentials in terms of the components of the strain tensor. Therefore, the orientation of each material point in the tissue could be specified in terms of a certain coordinate system. Typically, the cylindrical polar coordinates are chosen for this purpose since these models were first applied on arterial tissue. However, this type of formulation is inherently limited to specific kinematics and the applicability of these models is limited, as discussed by Holzapfel [19].

Humphrey and Yin [8] proposed a different approach to modelling anisotropy. The elastic potential is defined as the sum of an isotropic part ${\textstyle \Psi _{iso}}$, attributable to the behaviour of the extracellular matrix (ECM), and an anisotropic one ${\textstyle \Psi _{ani}}$, accountable for the contribution of the collagen fibres. The latter, in turn, is defined as the sum of the contribution of each family of fibres 1. Then, the elastic potential of the tissue is

 ${\displaystyle \Psi =\Psi _{iso}+\Psi _{ani},\quad {\textrm {with}}\quad \Psi _{ani}=\sum \limits _{i=1}^{k}\Psi _{f}^{i},}$
(2.1)

where ${\textstyle k}$ is the amount of fibre families considered. Here, ${\textstyle \Psi _{iso}}$ is defined in terms of the right Cauchy-Green deformation tensor ${\textstyle \mathbf {C} }$ invariants and should account for the ECM behaviour. ${\textstyle \Psi _{f}}$ depends on both ${\textstyle \mathbf {C} }$ and the fibre orientation, typically through the pseudo-invariants of the fibre orientation vectors, also known as Spencer invariants [20]. The incorporation of the matrix-fibre distinction and the fibre orientation confers a certain mechanistic character to this formulation, even if phenomenological models are used as basis to define the basic elastic potentials.

Over the years, many researchers have used this same concept to develop constitutive models to represent soft tissue behaviour. One of the most popular of these models is the HGO model, developed by Holzapfel, Gasser and Ogden to represent arterial tissue behaviour [19]. Generally, the isotropic elastic potential ${\textstyle \Psi _{iso}}$ is defined as an isotropic phenomenological hyperelastic function, while the definition of the anisotropic counterpart ${\textstyle \Psi _{ani}}$ depends on the histological structure of the particular tissue being studied. Numerous constitutive formulations based on this isotropic/anisotropic split have been developed in the past decade and a half, extending the models to account for characteristics such as viscoelasticity [21,22,23,24,25], softening [26,27,28,29,30], or both [31,32].

Biological tissues exhibit several known softening phenomena: the Mullins effect, preconditioning, permanent set and damage softening [29]. The Mullins effect is a phenomenon originally investigated and documented by Mullins et al. [33] in rubbers. This effect is phenomenologically characterized by the degradation of the elastic properties of rubber-like materials subjected to quasi-static cyclic tensile/compressive loading. The degradation is observed at strain levels below the maximum strain achieved in the history of deformation. This stress-softening is associated with internal damage at microscopic level and structural realignment in the material. However, there is no general agreement on its physical source. Diani et al. [34] review both phenomenological and macromolecular models put forth by researchers in the field to model the Mullins effect in rubber-like materials. Preconditioning or hysteresis is characterized by an initial continuous softening in the first cycles of a testing procedure which is identically repeated on the same tissues [11]. Softening in the tissue halts when it reaches a steady or “saturated” state, that is, it suffers no further changes in its internal structure. Permanent set is the residual stretch observed in tissues after unloading [35]. Finally, softening in biological tissues may occur as a result of damage or degradation in the material structure. Rupture, fissures and voids at microscopic level occur in the tissue as the (macroscopic) deformation proceeds, resulting in an overall softening behaviour of the tissue.

The latter type of softening is commonly described through models formulated in a continuum damage mechanics (CDM) framework. These damage formulations are based on the hyperelastic strain energy functions, as will be seen in section 2.3. The CDM bases will be described in this section, together with a brief review of finite-strain damage models developed up to date, with particular focus on those developed for biological applications. Then, a general finite-strain damage model that can be particularized for any hyperelastic strain energy function will be developed and its numerical implementation and main characteristics will be described.

Let us recall the models based on a split elastic potential, initially introduced by Humphrey and Yin [8]. They can be interpreted as an ad hoc mixing theory, tailored to each biomechanical application. Now, since the aim of this study is to set the bases for a constitutive formulation to represent soft tissue behaviour, we propose a more general approach by directly using the generalized mixing theory [9], which works as a constitutive equation manager. In this way, one can build a library of constitutive equations which are known to adequately represent certain biological tissues or their main components and, then, choose the relevant models and quantify their contribution according to each application. So, there should be no need to redefine a new constitutive equation for each new biomechanical problem. The mixing theory is no longer introduced in the definition of the elastic potential, instead, it is applied at stress level. Then, for a tissue composed by ${\textstyle n}$ distinct components (matrix and fibres), the second Piola-Kirchhoff stress in the tissue will be given by

 ${\displaystyle \mathbf {S} =\sum \limits _{i=1}^{n}{\textrm {v}}_{i}\mathbf {S} _{i}.}$
(2.2)

Here, ${\textstyle {\textrm {v}}_{i}}$ is the volumetric participation of each tissue component and ${\textstyle \mathbf {S} _{i}}$ is the second Piola-Kirchhoff stress of each of these components, computed from its own constitutive equation.

The mixing theory formulation used, developed in a finite strain framework, will be reviewed in section 2.4. Its capability in reproducing soft tissue behaviour will be discussed as well as its shortcomings in the present form. In particular, the modelling of anisotropic behaviour in tissues will be addressed and the numerical tools that can be used to explicitly account for this behaviour will be discussed.

(1) Humphrey and Yin [8] define a family of fibres as “a collection of locally parallel fibres with identical material properties”.

## 2.2 Hyperelasticity

Hyperelastic models are basically higher-order forms of linear elastic ones and, therefore, are a particular type of nonlinear elastic models. Yet, hyperelasticity accounts for both nonlinear material behaviour and large geometric changes. Rubber-like materials are the most common example of a hyperelastic material but many elastomers and most soft biological tissues are also modelled using a hyperelastic idealization.

The simplest elastic model is the linear elastic one, in which stress is directly proportional to strain through a linear relationship. However, in general, an elastic solid which undergoes large deformations will not follow linear elasticity. Many different constitutive relations can be developed to represent the behaviour of the elastic materials under the finite strain theory. In addition, the same constitutive behaviour, i.e., the same material model, can be written in several ways because of the different stress and deformation measures used in finite strain. These models are said to follow nonlinear elasticity and can be loosely grouped into [36,37,38]:

1. Kirchhoff material. It is a straightforward generalization of linear elasticity to finite strain, typically used for applications in which the large deformation effects are due to geometric nonlinearities. The most general model is defined through the relation ${\textstyle \mathbf {S} ={\boldsymbol {\boldsymbol {C}}}:\mathbf {E} ,}$ where ${\textstyle \mathbf {S} }$ is the second Piola-Kirchhoff stress tensor, ${\textstyle {\boldsymbol {\boldsymbol {C}}}}$ is the constitutive tensor and ${\textstyle \mathbf {E} }$ is the Green-Lagrange strain tensor.
2. Hypoelasticity. The constitutive equation is defined in terms of increments or rates, for example, ${\textstyle {\dot {\boldsymbol {\sigma }}}=f\left({\boldsymbol {\sigma }},\mathbf {d} \right)}$, where ${\textstyle {\dot {\boldsymbol {\sigma }}}}$ is the rate of the Cauchy stress tensor, ${\textstyle {\boldsymbol {\sigma }}}$ is the Cauchy stress tensor and ${\textstyle \mathbf {d} }$ is the rate of deformation tensor. These type of laws do not strictly reflect the path independence of elasticity. In addition, the derivation of objective stress rates 1 and the corresponding stiffness tensors is not trivial.
3. Cauchy elasticity. The stress is computed by means of a material response function, ${\textstyle {\boldsymbol {\sigma }}=\mathbf {\mathcal {G}} \left(\mathbf {F} \right)}$, where ${\textstyle \mathbf {F} }$ is the deformation gradient tensor. The material response function ${\textstyle \mathbf {\mathcal {G}} }$ has an explicit dependence on position, time and choice of reference configuration.
4. Hyperelasticity or Green elasticity. Stress is obtained from a scalar specific strain energy density function, also known as Helmholtz free energy function or elastic potential, ${\textstyle \Psi }$. This function is a measure of the energy stored in the deforming material such that, upon unloading, the energy is gradually released as the material recovers its original shape. The specific strain energy density function typically takes the form ${\textstyle \mathbf {S} ={\partial \Psi \left(\mathbf {E} \right)}/{\partial \mathbf {E} }}$, although expressions for different stress measures are easily obtained through appropriate transformations.

The finite strain framework in which hyperelasticity is developed assumes that both rotations and strains are large and, thus, the reference or undeformed configuration is significantly different to the present or deformed one (see Figure 1). A clear distinction must be made between them and it is of utmost importance to know in which configuration one is working in at all times.

The implications of developing constitutive models for FEA in a finite strain framework are detailed in section 2.2.1. Then, the constraints imposed by the quasi-incompressibility of soft biological tissues in the numerical modelling are addressed in section 2.2.2. Section 2.2.3 presents the derivation of the basic constitutive formulation of quasi-incompressible hyperelasticity. Several hyperelastic models and their applicability to reproducing soft tissue behaviour are outlined in section 2.2.4. Then, the formulation derived in section 2.2.3 is particularized for the neo-Hookean and Ogden models in sections 2.2.5 and 2.2.6, respectively. The details of the numerical implementation in the in-house FE code PLCd [1] are described for both, and, finally, the numerical and modelling limitations of these two formulations are discussed in section 2.2.7.

 Figure 1: Configurations of the continuous medium, adapted and reproduced with permission from Oliver [Oliver2003]. Point ${\displaystyle P}$ in the reference configuration ${\displaystyle \Omega _{0}}$ (at the initial time ${\displaystyle t_{0}}$) corresponds to point ${\displaystyle P'}$ in the present configuration ${\displaystyle \Omega _{t}}$ (at time ${\displaystyle t}$). ${\displaystyle \mathbf {X} }$ and ${\displaystyle \mathbf {x} }$ are the position vectors of this point for the reference and present configurations, respectively. F is the deformation gradient tensor.

(1) An objective stress rate is a time derivative of stress that is frame indifferent.

### 2.2.1 Finite strain framework

The great majority of solid mechanics FE solvers work with a Lagrangian mesh, and PLCd is no exception. This type of mesh is advantageous in computational solid mechanics because nodes and elements move with the material and, therefore, the problem boundaries and interfaces remain coincident with the element edges, simplifying their computational treatment. Also, the constitutive equations are evaluated always at the same material points because the integration points also move with the material and, thus, always coincide with the exact same point of the material. This is especially useful in history-dependent models such as damage, which will be addressed in section 2.3.

Then, the static equilibrium condition for an elemental discrete volume ${\textstyle e}$ in the reference configuration is given by

 ${\displaystyle {\boldsymbol {\boldsymbol {R}}}^{e}={\underset {\textrm {internalforces}}{\underbrace {\int _{V_{0}^{e}}\mathbf {S} ^{e}\mathbf {B} ^{e}\,d{V}_{0}} }}-{\underset {{\textrm {externalforces,}}{\boldsymbol {\boldsymbol {F}}}^{e}{\textrm {}}}{\underbrace {\left(\int _{V_{0}^{e}}\rho _{0}\mathbf {b} \,\mathbf {N} ^{e}\,d{V}_{0}+\int _{S_{0}^{e}}\mathbf {t} \,\mathbf {N} ^{e}\,d{S}_{0}\right)} }}=\mathbf {0} }$
(2.3)

which must be satisfied for any displacement increment ${\textstyle \delta \mathbf {U} ^{^{e}}}$. Here,

1. ${\textstyle {\boldsymbol {\boldsymbol {R}}}}$ is the vector of residual forces, which must be null;
2. ${\textstyle V_{0}}$, ${\textstyle S_{0}}$ and ${\textstyle \rho _{0}}$ are the reference volume, reference surface and reference density of the discrete volume, respectively;
3. ${\textstyle \mathbf {S} }$ is the second Piola-Kirchhoff stress tensor;
4. ${\textstyle \mathbf {B} }$ is the strain-displacement compatibility or transformation tensor, given by ${\textstyle {\mathbf {B} =\nabla ^{s}\mathbf {N} }}$;
5. ${\textstyle \mathbf {N} }$ is the shape function;
6. ${\textstyle \mathbf {b} }$ is the vector of body forces acting on the volume per unit of mass; and
7. ${\textstyle \mathbf {t} }$ is the vector of surface forces acting on the surface of the volume per unit of surface.

The local equilibrium equations of each elemental volume are assembled in order to obtain the global equilibrium equation,

 ${\displaystyle {\boldsymbol {\boldsymbol {R}}}={\underset {e}{\mathrm {\mathbf {A} } }}\left\{\int _{V_{0}^{e}}\mathbf {S} ^{e}\mathbf {B} ^{e}\,dV_{0}\right\}-{\boldsymbol {\boldsymbol {F}}}=\mathbf {0} ,}$
(2.4)

where ${\textstyle \mathrm {\mathbf {A} } }$ is a lineal operator that represents the sum of the different force components, according to the position and direction of the local contributions.

This equation may include nonlinearities, which can be of material or geometric origin. Material nonlinearities arise when the stress-strain behaviour is nonlinear due to the particular constitutive equation of the material. Constitutive nonlinearity is strictly due to the changes in the material properties during its mechanical behaviour. It is directly reflected in the constitutive tensor ${\textstyle {\boldsymbol {\boldsymbol {C}}}}$ and, thus, the nonlinearity is introduced in 2.4 through the stress tensor ${\textstyle \mathbf {S} }$. Geometric nonlinearities appear when there are changes in the geometry which have a significant effect on the load-deformation behaviour. In large displacement cases, the large translations and rotations taking place induce changes in the local reference system of the different solid points. This introduces nonlinearities in the displacement-strain compatibility tensor ${\textstyle \mathbf {B} }$. Large strains, in addition to the large displacement effects, also introduce nonlinearities in the strain field due to the change in configuration of the solid. This change also modifies the constitutive tensor, which results in a nonlinear dependency of stresses on strains.

The reader is referred to reference textbooks in the FEM field such as Zienkiewicz and Taylor [5], Bathe [39], Crisfield [37] and Belytschko et al. [40] for detailed information on the numerical bases of nonlinear FE analysis.

FE solutions for problems which use Lagrangian meshes are typically divided into total Lagrangian (TL) and updated Lagrangian (UL) formulations. In both cases the dependent variables are functions of the material coordinates and time, i.e., they use Lagrangian descriptions. However, in the TL scheme all variables are referred to the reference (initial) configuration at time 0 whilst in the UL one they are referred to the current (deformed) configuration at time ${\textstyle t}$. Thus, in a TL formulation integrals are taken over the reference configuration and derivatives are calculated with respect to the material coordinates. In contrast, variables are integrated over the current configuration and derived with respect to spatial coordinates in an UL formulation. Nonetheless, both formulations include nonlinear effects due to large displacements, large rotations and large strains. Whether the large strain behaviour is modelled correctly or not will depend on the constitutive relations used and the results will be the same for both formulations if the appropriate constitutive models are defined in each formulation.

A decisive factor in choosing one formulation over the other is the computational efficiency. In the present study, the use of mixing theory might require the definition of angles to indicate the fibre alignment in the matrix. An UL approach would require updating the angle information, as well as the integration volumes, each time the equilibrium of forces were to be enforced. This was regarded as a considerable computational drawback and the decision was made to use a TL framework. However, because one might need to define certain constitutive equations in the reference configuration, a partially updated Lagrangian (pUL) framework was implemented. In the pUL framework, all steps of the resolution scheme are performed in the reference configuration, except for the computation of the stress tensor through the constitutive equation, which is done in the present configuration (see Figure 2). To switch between the two configurations, push-forward and pull-back operators, denoted as ${\textstyle {\overrightarrow {\phi }}}$ and ${\textstyle {\overleftarrow {\phi }}}$, respectively, are required. These are defined as

 ${\displaystyle {\begin{array}{l}{\boldsymbol {\tau }}=J{\boldsymbol {\sigma }}={\overrightarrow {\phi }}\left(\mathbf {S} \right)=\mathbf {F} \cdot \mathbf {S} \cdot \mathbf {F} ^{T},\\\mathbf {e} ={\overrightarrow {\phi }}\left(\mathbf {E} \right)=\mathbf {\mathbf {F} } ^{-T}\cdot \mathbf {E} \cdot \mathbf {F} ^{-1},\\\\\mathbf {S} ={\overleftarrow {\phi }}\left({\boldsymbol {\tau }}\right)=\mathbf {F} ^{-1}\cdot {\boldsymbol {\tau }}\cdot \mathbf {F} ^{-T}=J\,\mathbf {F} ^{-1}\cdot {\boldsymbol {\sigma }}\cdot \mathbf {F} ^{-T},\\\mathbf {E} ={\overleftarrow {\phi }}\left(\mathbf {e} \right)=\mathbf {\mathbf {F} } ^{T}\cdot \mathbf {e} \cdot \mathbf {F} ,\end{array}}}$
(2.5)

where ${\textstyle {\boldsymbol {\tau }}}$ is the Kirchhoff stress tensor, ${\textstyle J=\det {\mathbf {F} }}$ is the Jacobian determinant and ${\textstyle \mathbf {e} }$ is the Euler-Almansi strain tensor.

 Figure 2: Scheme of the total Lagrangian and partially updated Lagrangian formulations implemented in PLCd [1]. The subindex ${\displaystyle k}$ indicates iteration number in the present load increment. The definition of each term is available in the notation list.

### 2.2.2 Quasi-incompressibility and hybrid elements

The near-incompressibility characteristic of many tissues is attributed to the high volume fraction of water present in most soft tissues [11]. Hence, soft tissue behaviour is typically represented through quasi-incompressible hyperelastic models, albeit the adequacy of the quasi-incompressibility hypothesis of soft tissues has been debated in some cases [41,42].

Assuming soft tissue behaviour can be modelled by means of quasi-incompressible hyperelasticity implies that any hydrostatic pressure can be applied on said material without changing its shape and, thus, stress cannot be uniquely determined from strains. In this case, the hyperelastic constitutive law only specifies the deviatoric stresses and the hydrostatic stresses must be calculated by solving the equilibrium equations with the appropriate boundary conditions (assuming a quasi-static loading case). Thus, the specific strain energy density function in the reference configuration is split into

 ${\displaystyle \Psi =\Psi _{vol}+{\widetilde {\Psi }},}$
(2.6)

where ${\textstyle \Psi _{vol}}$ corresponds to the volume-preserving or volumetric part and ${\textstyle {\widetilde {\Psi }}}$, to the isochoric or deviatoric one [43].

In addition, a quasi-incompressible material has a quasi-infinite bulk elastic modulus, which results in an ill-conditioned stiffness matrix and locking problems when computing the solution using standard displacement-based FE formulations.

Gadala [44] reviews in his work the formulations proposed to overcome these difficulties, detailing the advantages and disadvantages of one approach over another. All of them are based on the decomposition of the deformation gradient into a deviatoric or isochoric part and a volumetric or dilatational one. The most popular approaches used in literature can be grouped into:

1. Multifield or mixed principles. The variational principle behind the FE formulation not only includes the displacement field, but also the volumetric deformation and/or pressure fields [45,43]. In their most general form, they have the disadvantage of introducing a considerable number of unknowns into the problem, which is computationally costly. However, static condensation at element level can eliminate the pressure degrees of freedom and, thus, this problem is overcome [39].
2. Reduced selective integration penalty approach. The displacement field corresponding to the deviatoric deformation is integrated normally but the volumetric terms are integrated using low order rules [46]. This reduces the over-stiffening caused by volumetric locking. Unfortunately, this approach does not work well with elements subjected to large strains or highly distorted elements [44].

In the present work, a two-field approach has been used. The u/p formulation described by Crisfield [37] and Bathe [39], consisting in a mixed interpolation of the displacements and pressure has been implemented into PLCd [1]. This formulation is based on the classical displacement FEM but includes an additional unknown variable, the hydrostatic stress distribution or pressure, which is interpolated separately from the displacement variable [46]. Then, the equations of motion for an elemental discrete volume are

 ${\displaystyle \left[{\begin{array}{cc}\mathbf {K} _{UU}&\mathbf {K} _{UP}\\\mathbf {K} _{PU}&\mathbf {K} _{PP}\end{array}}\right]\left[{\begin{array}{c}\mathbf {U} \\\mathbf {P} \end{array}}\right]=\left[{\begin{array}{c}{\boldsymbol {\mathbb {F} }}_{0}^{ext}\\0_{\,}\end{array}}\right]-\left[{\begin{array}{c}{\boldsymbol {\mathbb {F} }}_{0,\,U}^{int}\\{\boldsymbol {\mathbb {F} }}_{0,\,P}^{int}\end{array}}\right],}$
(2.7)

where ${\textstyle \mathbf {U} }$ is the displacement vector, ${\textstyle \mathbf {P} }$ is the pressure vector, ${\textstyle \mathbf {K} _{\left\{\bullet \right\}\left\{\circ \right\}}}$ are the stiffness matrices and ${\textstyle {{\boldsymbol {\mathbb {F} }}_{0}^{ext}}}$ and ${\textstyle {\boldsymbol {\mathbb {F} }}_{0,\,\left\{\bullet \right\}}^{int}}$ are the vector of forces corresponding to the external and internal loads, respectively. All variables are nodal and given in the reference configuration. The stiffness matrices are defined as ${\textstyle {\mathbf {K} _{\left\{\bullet \right\}\left\{\circ \right\}}=\partial \,{\boldsymbol {\mathbb {F} }}_{0,\,\left\{\bullet \right\}}^{int}/\partial \left\{\circ \right\}}}$, which results in

 ${\displaystyle \mathbf {K} _{UU}=\int _{V_{\mathit {_{0}}}}\mathbf {B} _{0}^{T}:{\boldsymbol {\boldsymbol {C}}}^{tan}:\mathbf {B} _{0}\,dV_{\mathit {_{0}}}+\int _{V_{\mathit {_{0}}}}\mathbf {B} _{0,\;NL}^{T}\cdot \mathbf {S} \cdot \mathbf {B} _{0,\;NL}\;dV_{\mathit {_{0}}},}$
(2.8)

 ${\displaystyle \mathbf {K} _{UP}=\int _{V_{\mathit {_{0}}}}\mathbf {B} _{0,\;NL}^{T}\cdot \mathbf {g} \,\mathbf {h} _{p}^{T}\;dV_{\mathit {_{0}}}=\mathbf {K} _{PU}^{T}}$
(2.9)

and

 ${\displaystyle \mathbf {K} _{PP}=-\int _{V_{\mathit {_{0}}}}{\frac {1}{\kappa }}\,\mathbf {h} _{P}\,\mathbf {h} _{P}^{T}\;dV_{\mathit {_{0}}}.}$
(2.10)

Here ${\textstyle \mathbf {B} _{0}}$ and ${\textstyle \mathbf {B} _{0,\;NL}}$ are the classical linear and nonlinear strain-displacement compatibility or transformation tensors, respectively [2]. ${\textstyle {\boldsymbol {\boldsymbol {C}}}^{tan}}$ is the tangent constitutive tensor, ${\textstyle \mathbf {h} _{P}}$ is the vector of pressure shape functions, ${\textstyle \kappa }$ is the bulk modulus and ${\textstyle \mathbf {g} }$ is a vector which in matrix form is the pressure-related constitutive tensor ${\textstyle \mathbf {G} }$. The internal forces ${\textstyle {\boldsymbol {\mathbb {F} }}_{\left\{\bullet \right\}}^{int}}$ are computed as

 ${\displaystyle {\boldsymbol {\mathbb {F} }}_{U}^{int}=\int _{V_{\mathit {_{0}}}}\mathbf {B} _{0,\;NL}^{T}\cdot \mathbf {S} \;dV_{\mathit {_{0}}}}$
(2.11)

and

 ${\displaystyle {\boldsymbol {\mathbb {F} }}_{P}^{int}=\int _{V_{\mathit {_{0}}}}{\frac {1}{\kappa }}\left({\overline {p}}-{\widetilde {p}}\right){\frac {\partial {\widetilde {p}}}{\partial \mathbf {P} }}\;dV_{\mathit {_{0}}},}$
(2.12)

where ${\textstyle {\overline {p}}}$ is the pressure obtained from the displacement field and ${\textstyle {\widetilde {p}}}$ is the total element pressure obtained by independent interpolation, ${\textstyle {\widetilde {p}}=\mathbf {h} _{P}\cdot \mathbf {P} }$. Then, ${\textstyle {\boldsymbol {\boldsymbol {C}}}^{tan}}$ and ${\textstyle \mathbf {G} }$ are constitutive tensors relating strain and pressure, respectively, to stress, i.e.,

 ${\displaystyle \mathbf {\dot {S}} ={\boldsymbol {\boldsymbol {C}}}^{tan}:\mathbf {\dot {E}} +\mathbf {G} \cdot {\dot {\mathbf {P} }}.}$
(2.13)

The pressure constitutive tensor ${\textstyle \mathbf {G} }$ is given by

 ${\displaystyle \mathbf {G} =-J\,\mathbf {C} ^{-1},}$
(2.14)

where ${\textstyle \mathbf {C} }$ is the right Cauchy-Green deformation tensor that is defined in terms of the deformation gradient tensor as ${\textstyle \mathbf {C} =\mathbf {F} ^{T}\cdot \mathbf {F} }$. Since the volumetric part of the Cauchy stress tensor is given directly by the hydrostatic pressure 1, ${\textstyle {\boldsymbol {\sigma }}_{vol}=-p\,\mathbf {I} }$, the pressure constitutive tensor is, in fact, a pull-back operation 2.5 of this term. So, ${\textstyle J}$ converts the Cauchy stress to Kirchhoff stress (${\textstyle {\boldsymbol {\tau }}=J{\boldsymbol {\sigma }}}$) and, through the inverse of the right Cauchy-Green deformation tensor, the pull-back operation is completed.

Now, the pressure obtained from the displacement field, ${\textstyle {\overline {p}}}$, is defined positive in compression and computed as

 ${\displaystyle {\overline {p}}=-{\frac {\partial \Psi _{vol}}{\partial J}},}$
(2.15)

where ${\textstyle \Psi _{vol}}$ can take different forms according to different authors [39,37,43,47], but is always defined in terms of the Jacobian determinant ${\textstyle J}$, which is a measure of the change of volume, and the bulk modulus ${\textstyle \kappa }$. Possible functions of ${\textstyle \Psi _{vol}}$ and their corresponding function for the pressure ${\textstyle {\overline {p}}}$ are

 ${\displaystyle {\begin{array}{l}\Psi _{vol}^{(a)}={\frac {1}{2}}\kappa \left(J-1\right)^{2}\;\;\Rightarrow \quad {\overline {p}}=-\kappa \left(J-1\right),\\[4mm]\Psi _{vol}^{(b)}={\frac {1}{2}}\kappa \left(\ln J\right)^{2}\;\;\Rightarrow \quad {\overline {p}}=-\kappa {\dfrac {\ln J}{J}},\quad {\textrm {and}}\\[4mm]\Psi _{vol}^{(c)}={\frac {1}{2}}\kappa \left(J^{2}-2\ln J\right)\;\;\Rightarrow \quad {\overline {p}}=-\kappa {\dfrac {\left(J^{2}-1\right)}{J}}.\end{array}}}$
(2.16)

In the first case, ${\textstyle \Psi _{vol}^{(a)}}$, the value of the bulk modulus ${\textstyle \kappa }$ directly dictates the value of the pressure ${\textstyle {\overline {p}}}$ when there is no change in volume, i.e., ${\textstyle J=0}$ (see Figure 3). In addition, the tetrahedral elements that will be discussed later are very sensitive to the choice of the function ${\textstyle \Psi _{vol}}$. In this sense, the best choice is the one with the term ${\textstyle \left(J-1\right)}$ because it usually has a stabilizing effect [48]. For these reasons, the first expression in 2.16 of the volumetric specific strain energy density function has been used in PLCd [1]. Then, the internal force corresponding to the volumetric part is reduced to

 ${\displaystyle {\boldsymbol {\mathbb {F} }}_{P}^{int}=-\int _{V_{\mathit {_{0}}}}\left(\left(J-1\right)+{\frac {\widetilde {p}}{\kappa }}\right)\mathbf {h} _{P}\,dV_{\mathit {_{0}}},}$
(2.17)

where ${\textstyle {\partial {\widetilde {p}}}/{\partial \mathbf {P} }=\mathbf {h} _{P}}$ has been taken into account.

 Figure 3: Pressure ${\displaystyle {\overline {p}}}$ (positive in compression) vs. the Jacobian determinant ${\displaystyle J}$ for the different volumetric specific strain energy density functions ${\displaystyle \Psi _{vol}}$ given in 2.16 in terms of the bulk modulus ${\displaystyle \kappa }$.

Figure 4 shows two possible configurations of mixed elements based on the serendipitous hexahedral and the tetrahedral families of elements already implemented in PLCd [1]. Both element types satisfy the inf-sup condition 2 and, therefore, will produce robust and reliable solutions. Figure 5 shows alternative configurations of the distribution of the pressure nodes in which these nodes do not coincide with the displacement nodes. Such elements are often referred to as hybrid elements and have the advantage of allowing for static condensation of the pressure variables in the equation of motion. Isolating ${\textstyle \mathbf {P} }$ from the second line of 2.7 yields

 ${\displaystyle \mathbf {P} =-\mathbf {K} _{PP}^{-1}\cdot \left(\mathbb {F} _{P}^{int}+\mathbf {K} _{UP}^{T}\cdot \mathbf {U} \right),}$
(2.18)

which, replaced into the first line of 2.7, results in the condensed equation of motion to be solved,

 ${\displaystyle {\underset {\mathbf {\displaystyle {\overline {K}}} }{\underbrace {\left(\mathbf {K} _{UU}-\mathbf {K} _{UP}\cdot \mathbf {K} _{PP}^{-1}\cdot \mathbf {K} _{UP}^{T}\right)} }}\cdot \mathbf {U} ={\boldsymbol {\mathbb {F} }}_{0}^{ext}-{\underset {{\overline {\displaystyle {\boldsymbol {\mathbb {F} }}}}^{int}}{\underbrace {\left({\boldsymbol {\mathbb {F} }}_{U}^{int}-\mathbf {K} _{UP}\cdot \mathbf {K} _{PP}^{-1}\cdot {\boldsymbol {\mathbb {F} }}_{P}^{int}\right)} }}.}$
(2.19)

The formulation obtained is expressed solely in terms of displacement variables, and the new stiffness tensor ${\textstyle {\overline {K}}}$ and internal force vector ${\textstyle {\overline {\boldsymbol {\mathbb {F} }}}^{int}}$ differ from the classic displacement FE formulation in that they have the pressure-related subtracting terms. Therefore, the same FE solver can be used because new variables have not been added, only the stiffness tensor and internal force vector require modification.

 Figure 4: Possible configurations for mixed u/p elements that satisfy the inf-sup condition [39]: elements Q2P1c (left) and P2P1c (right). The blue circles indicate displacement nodes and the orange crosses indicate pressure nodes. Nodes are numbered according to the local numbering criterion used in PLCd [1].
 Figure 5: Possible configurations for mixed u/p or hybrid elements that allow for static condensation of the pressure variables in the equation of motion. Elements Q1P0 (top left) and T1P0 (top right) do not satisfy the inf-sup condition. Elements Q2P0 (bottom left) and T2P0 (bottom right) satisfy the inf-sup condition but the constant pressure assumption may require fine discretization [39]. The blue circles indicate displacement nodes and the orange crosses indicate pressure nodes. Nodes are numbered according to the local numbering criterion used in PLCd [1].

Now, not all possible configurations of hybrid elements shown in Figure 5 perform equally. Element Q1P0 does not satisfy the inf-sup condition but is capable of predicting reasonably good displacements and is commonly used because of its simplicity. However, due to the constant pressure assumption, stress predictions may be inaccurate [39]. Element T1P0 is known to exhibit volumetric locking for nearly incompressible deformations and most authors recommend directly avoiding its use in problems of this type [37,49,38]. Nonetheless, tetrahedral linear elements present considerable advantages in terms of computational cost and straightforward meshing of complex geometries. This motivated some authors [50,51,52] to develop the averaged nodal pressure tetrahedral element, which overcomes the volumetric problems of the standard T1P0 element. This option was studied for PLCd [1] but eventually discarded since it was seen that quadratic tetrahedrons do not entail an excessively high computational cost for the applications we could consider modelling with this code. Elements Q2P0 and T2P0 both satisfy the inf-sup condition but, again, fine discretization may be required in order to accurately predict stresses due to the constant pressure assumption. To overcome this problem, linear or quadratic pressure variable distributions can be introduced but, to obtain adequate results, switching to higher order displacement distributions is also required [39]. Thus, Q1P0, Q2P0 and T2P0 elements are deemed to be sufficient for the purposes of this study.

When a constant pressure distribution is considered, the condensed equation of motion 2.19 is further simplified. There is a single pressure variable, therefore, ${\textstyle \mathbf {P} ={\widetilde {p}}=p}$, and the pressure shape function vector is reduced to ${\textstyle \mathbf {h} _{P}=1}$. Then, the stiffness matrices 2.9 and 2.10, and internal force vector 2.12 become

 ${\displaystyle \mathbf {K} _{UP}=-\int _{V_{\mathit {_{0}}}}\mathbf {B} _{0,\,NL}^{T}\cdot J\,\mathbf {C} ^{-1}\,dV_{\mathit {_{0}}}=\mathbf {K} _{PU}^{T},}$
(2.20)

 ${\displaystyle \mathbf {K} _{PP}=-\int _{V_{\mathit {_{0}}}}{\frac {1}{\kappa }}\,dV_{\mathit {_{0}}},\quad {\textrm {and}}}$
(2.21)

 ${\displaystyle {\boldsymbol {\mathbb {F} }}_{P}^{int}=-\int _{V_{\mathit {_{0}}}}\left(\left(J-1\right)+{\frac {p}{\kappa }}\right)\,dV_{\mathit {_{0}}}.}$
(2.22)

The definition of the pressure constitutive tensor 2.14 has been taken into account in 2.20.

Figure 6 summarizes the condensed mixed u/p formulation implemented in PLCd [1] for the total Lagrangian framework. Extension to the partially updated Lagrangian framework is straightforward and requires no additional changes (see Figure 2).

 Figure 6: Scheme of the condensed mixed u/p or hybrid formulation implemented in PLCd [1] for a total Lagrangian framework (reference configuration). The subindex ${\displaystyle k}$ indicates iteration number in the present load increment. The definition of each term is available in the notation list.

(1) Here, ${\displaystyle \mathbf {I} }$ is the second-order identity tensor, which is defined as ${\displaystyle \left[\mathbf {I} \right]_{ij}=\delta _{ij}}$, being ${\displaystyle \delta _{ij}}$ the Kronecker delta function.

(2) The inf-sup condition is the basic mathematical criterion that determines whether a mixed finite element discretization is stable and convergent and, hence, will yield a reliable solution [39].

### 2.2.3 Thermodynamic basis of hyperelastic formulations

The Axiom of Local State postulates that, in continuous media, the thermodynamic state at a given point and time instant is completely defined by a certain number of variables at that time instant. The physical phenomena can be described with precision depending on the nature and number of state variables chosen for that. The process defined in this way will be thermodynamically admissible if, at any time instant, the Clausius-Duhem inequality 1 is satisfied. The state variables, also named thermodynamic or independent variables, are the observable and the internal variables. Internal variables describe the internal structure of the medium, which is hidden to the eye of the external observer who simply sees a “black box”. It must be pointed out, though, that the notion of internal depends upon the considered level of observation: a variable considered as internal from the point of view of macroscopic observation will probably be an observable variable from a mesoscopic point of view. Thus, the tensor and physical nature of the internal variables should be known. Furthermore, these variables are, in practice, measurable but not controllable.

It follows that the dependent values in the continuous medium such as stress are, at any given time instant, a function of both the values of the independent observable variables and the internal variables. In the particular case under consideration, the observable variables are strain, ${\textstyle \mathbf {E} }$, and temperature, ${\textstyle \theta }$, whilst no internal variables are considered. We anticipate that, in the following section, when damage is introduced, there will be an internal variable associated with the damage in the material, ${\textstyle D}$.

Once the state variables have been defined, one postulates the existence of a thermodynamic potential from which the state laws can be derived. This potential must be a continuous scalar function, concave with respect to the temperature, convex with respect to the other state variables and containing the origin. Taking the Helmholtz free energy, one has ${\textstyle \Psi =\Psi \left(\mathbf {E} ,\theta \right)}$. Then, the local form of the second law of thermodynamics in terms of internal specific energy density is expressed through the Clausius-Duhem inequality in reference configuration [54] as

 ${\displaystyle \Xi =-{\dot {\Psi }}-\eta \,{\dot {\theta }}+\mathbf {S} :\mathbf {\dot {E}} -{\frac {1}{\theta }}\;{\boldsymbol {q}}_{0}\cdot \nabla \theta \geq {0},}$
(2.23)

where ${\textstyle \eta }$ is the specific entropy, ${\textstyle {\boldsymbol {q}}_{0}}$ is the heat flux vector and ${\textstyle \Xi }$ is the dissipation. Considering solely purely deformation processes where there does not intervene any changes in entropy or temperature, and taking into account that hyperelastic materials have null internal dissipation due to the reversible character of their loading processes, this equation is reduced to

 ${\displaystyle \Xi =-{\dot {\Psi }}+\mathbf {S} :\mathbf {\dot {E}} =0.}$
(2.24)

Then, the time derivative of the energy function is described in terms of the second Piola-Kirchhoff stress tensor and its conjugate, the rate of Green-Lagrange strain as ${\textstyle {\dot {\Psi }}=\mathbf {S} :\mathbf {\dot {E}} }$. The rate of energy, which is totally due to deformation, is equal to the stress power. Expanding the derivative of the energy function leads to

 ${\displaystyle {\frac {\partial \Psi \left(\mathbf {E} \right)}{\partial \mathbf {E} }}:{\frac {\partial \mathbf {E} }{\partial t}}=\mathbf {S} :\mathbf {\dot {E}} ,}$
(2.25)

which shows that, as expected, ${\textstyle \mathbf {S} }$ derives from a strain energy function. It is often more convenient to use the right Cauchy-Green deformation tensor instead of the Green-Lagrange strain tensor. Introducing the relation ${\textstyle \mathbf {E} =\left(\mathbf {C} -\mathbf {I} \right)/2}$, the stress results in

 ${\displaystyle \mathbf {S} ={\frac {\partial \Psi \left(\mathbf {E} \right)}{\partial \mathbf {E} }}=2{\frac {\partial \Psi \left(\mathbf {C} \right)}{\partial \mathbf {C} }}.}$
(2.26)

This relation can be expressed using other stress tensors such as the first Piola-Kirchhoff stress tensor, ${\textstyle \mathbf {T} }$. If one recalls that, like ${\textstyle \mathbf {S} }$ and ${\textstyle \mathbf {\dot {E}} }$, the tensors ${\textstyle \mathbf {T} }$ and ${\textstyle \mathbf {\dot {F}} }$ are work conjugates, then

 ${\displaystyle \mathbf {T} ={\frac {\partial \Psi \left(\mathbf {F} \right)}{\partial \mathbf {F} }}.}$
(2.27)

In any case, the elastic potential ${\textstyle \Psi }$ must fulfil:

1. Normalization condition: the function must be zero when the material is completely unloaded, ${\textstyle \Psi \left(\mathbf {F} =\mathbf {I} \right)=0}$ .
2. The energy function must grow monotonously with deformation, ${\textstyle \Psi \left(\mathbf {F} \right)\geq {0}}$.

Additionally, a reversible material (with no internal energy dissipation) must also satisfy:

1. The work performed by the forces must be independent of the path followed, ${\textstyle \int _{\Gamma _{1}}\mathbf {S} :d\mathbf {E} =\int _{\Gamma _{2}}\mathbf {S} :d\mathbf {E} }$.
2. For any closed deformation cycle, the deformation work must be zero, ${\textstyle {\oint \mathbf {S} :d\mathbf {E} =0}}$.

Now, introducing the additive split of the specific strain energy density function 2.6 into 2.26, and considering the first volumetric specific strain energy density function in 2.16, yields the constitutive equation

 ${\displaystyle \mathbf {S} =2{\frac {\partial {\widetilde {\Psi }}\left(\mathbf {C} \right)}{\partial \mathbf {C} }}+2{\frac {\partial \Psi _{vol}}{\partial J}}{\frac {\partial J}{\partial \mathbf {C} }}=2{\frac {\partial {\widetilde {\Psi }}\left(\mathbf {C} \right)}{\partial \mathbf {C} }}-p\,J\mathbf {C} ^{-1}.}$
(2.28)

Here, the relation ${\textstyle J=\left[\det \mathbf {C} \,\right]^{1/2}}$ has been considered in the derivation of ${\textstyle {\partial J/\partial \mathbf {C} =-J\,\mathbf {C} ^{-1}/\,2}}$. This constitutive equation must now be completed by defining an appropriate deviatoric part of the specific strain energy density function, ${\textstyle {\widetilde {\Psi }}}$. Many different functions have been proposed over the years in order to adequately represent the behaviour of the material being modelled. These must satisfy the aforementioned conditions of the elastic potential in addition to the principle of material frame-indifference 2. The isotropic functions typically used in biomechanical applications will be reviewed in the following section.

An isotropic material exhibits identical behaviour in all material directions. In mathematical terms, isotropy requires 2.26 to be independent of the chosen material axes. Thus, the specific strain energy density function must be defined solely in terms of the invariants of the right Cauchy-Green deformation tensor,

 ${\displaystyle \Psi \left(\mathbf {C} \right)=\Psi \left(I_{\mathbf {C} }^{(1)},I_{\mathbf {C} }^{(2)},I_{\mathbf {C} }^{(3)}\right).}$
(2.29)

The invariants of ${\textstyle \mathbf {C} }$ are defined as

1. First invariant: ${\textstyle I_{\mathbf {C} }^{(1)}=Tr\left(\mathbf {C} \right)=C_{ii}}$,
2. Second invariant: ${\textstyle I_{\mathbf {C} }^{(2)}=\left(\mathbf {C} :\mathbf {C} -I_{\mathbf {C} }^{2}\right)/\,2=\left(C_{ij}C_{ji}-C_{ii}^{2}\right)/\,2}$,
3. Third invariant: ${\textstyle I_{\mathbf {C} }^{(3)}=\det \left(\mathbf {C} \right)=J^{2}}$,

and the corresponding partial derivatives, which will prove useful in the derivation of the different hyperelastic constitutive equations, are

1. ${\textstyle {\partial I_{\mathbf {C} }^{(1)}}/\,{\partial \mathbf {C} }=\mathbf {I} }$,
2. ${\textstyle {\partial I_{\mathbf {C} }^{(2)}}/\,{\partial \mathbf {C} }=2\,\mathbf {C} }$   and
3. ${\textstyle {\partial I_{\mathbf {C} }^{(3)}}/\,{\partial \mathbf {C} }=J^{2}\,\mathbf {C} ^{-1}}$.

Expressions of isotropic specific strain energy density functions given in the present configuration are generally defined in terms of the invariants of the left Cauchy-Green deformation tensor ${\textstyle \mathbf {B} }$. This tensor, also known as Finger deformation tensor, is defined as ${\textstyle \mathbf {B} =\mathbf {F} \cdot \mathbf {F} ^{T}}$ and is known to have identical invariants to those of the right Cauchy-Green deformation tensor: ${\textstyle I_{\mathbf {C} }^{(1)}=I_{\mathbf {B} }^{(1)}}$, ${\textstyle I_{\mathbf {C} }^{(2)}=I_{\mathbf {B} }^{(2)}}$ and ${\textstyle I_{\mathbf {C} }^{(3)}=I_{\mathbf {B} }^{(3)}}$.

Finally, let us evaluate the rate of change of the constitutive equation 2.28, which can be expressed as

 ${\displaystyle \mathbf {\dot {S}} =\underbrace {4{\frac {\partial ^{2}\Psi \left(\mathbf {C} \right)}{\partial \mathbf {C} \;\partial \mathbf {C} }}} _{\displaystyle {\boldsymbol {\boldsymbol {C}}}^{tan}}:\mathbf {\dot {C}} .}$
(2.30)

Here, ${\textstyle {\boldsymbol {\boldsymbol {C}}}^{tan}}$ is a symmetric fourth-order tensor named tangent constitutive tensor, tangent stiffness or elasticity tensor. It can be noted that, even though this expression is similar to Hooke's Law (${\textstyle {{\boldsymbol {\sigma }}={\boldsymbol {\boldsymbol {C}}}:{\boldsymbol {\varepsilon }}}}$), here the tangent constitutive tensor is not constant but a function of strain, ${\textstyle {{\boldsymbol {\boldsymbol {C}}}^{tan}\left(\mathbf {C} \right)}}$ and it relates the material time derivatives of strain and stress. Because the stress and strain variables belong to the reference configuration, this tensor is known as the Lagrangian or material elasticity tensor. The equivalent tensor in current configuration is obtained if a push-forward operation is performed on the time derivatives of stress and strain. This yields the Eulerian or spatial elasticity tensor, ${\textstyle {\boldsymbol {\boldsymbol {c}}}^{tan}}$, which relates the Lie derivative 3 of the Euler-Almansi strain to the Lie derivative of the Kirchhoff stress,

 ${\displaystyle {\mathcal {L_{\phi }}}\left({\boldsymbol {\tau }}\right)=J{\boldsymbol {\boldsymbol {c}}}^{tan}:{\mathcal {L_{\phi }}}\left(\mathbf {e} \right).}$
(2.31)

The spatial elasticity tensor is, in fact, the push-forward of the material one which, in index notation, is written as

 ${\displaystyle {\boldsymbol {c}}_{_{ijkl}}^{tan}=F_{iI}\cdot F_{jJ}\cdot F_{kK}\cdot F_{lL}\cdot {\boldsymbol {C}}_{_{IJKL}}^{tan}.}$
(2.32)

The tangent constitutive tensor is used to obtain the stiffness tensor ${\textstyle \mathbf {K} _{UU}}$ by means of 2.8. The correct implementation of the tangent constitutive tensor is of particular importance in highly nonlinear models because an error in its numerical formulation can lead to lack of convergence in the calculation, especially in problems involving triaxial stress states.

Introducing the additive split of the specific strain energy density function 2.6 into 2.30 yields the material tangent constitutive tensor

 ${\displaystyle {\boldsymbol {\boldsymbol {C}}}^{tan}={\widetilde {\boldsymbol {\boldsymbol {C}}}}^{tan}+{\boldsymbol {\boldsymbol {C}}}_{vol}^{tan}=4{\frac {\partial ^{2}{\widetilde {\Psi }}}{\partial \mathbf {C} \;\partial \mathbf {C} }}+4{\frac {\partial ^{2}{\Psi }_{vol}}{\partial \mathbf {C} \;\partial \mathbf {C} }}.}$
(2.33)

Considering the first volumetric specific strain energy density function in 2.16 for the volumetric tangent constitutive tensor results in

 ${\displaystyle {\boldsymbol {\boldsymbol {C}}}_{vol}^{tan}=2p{\frac {\partial \left(J\mathbf {C} ^{-1}\right)}{\partial \mathbf {C} }}+2J\mathbf {C} ^{-1}\otimes p{\frac {\partial p}{\partial \mathbf {C} }},}$
(2.34)

where the operator ${\textstyle \otimes }$ is used to indicate the open product between the two second-order terms. Further developing this expression in terms of the right Cauchy-Green deformation tensor invariants produces

 ${\displaystyle {\boldsymbol {\boldsymbol {C}}}_{vol}^{tan}=-p\,{\big (}I_{\mathbf {C} }^{(3)}{\big )}^{1/2}\left(\mathbf {C} ^{-1}\otimes \mathbf {C} ^{-1}-2{\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}\right)+\kappa \;I_{\mathbf {C} }^{(3)}\mathbf {C} ^{-1}\otimes \mathbf {C} ^{-1}.}$
(2.35)

Here, the fourth-order tensor ${\textstyle {\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}}$ is defined as

 ${\displaystyle \left[{\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}\right]_{IJKL}={\frac {1}{2}}\left(\left[\mathbf {C} ^{-1}\right]_{IK}\left[\mathbf {C} ^{-1}\right]_{JL}+\left[\mathbf {C} ^{-1}\right]_{IL}\left[\mathbf {C} ^{-1}\right]_{JK}\right).}$
(2.36)

Performing the push-forward operation 2.32 on 2.35, yields the volumetric part of the spatial elasticity tensor,

 ${\displaystyle {\boldsymbol {\boldsymbol {c}}}_{vol}^{tan}=-p\,{\big (}I_{\mathbf {B} }^{(3)}{\big )}^{1/2}\left(\,\mathbf {I} \otimes \mathbf {I} -2{\boldsymbol {\boldsymbol {I}}}\,\right)+\kappa \;I_{\mathbf {B} }^{(3)}\mathbf {I} \otimes \mathbf {I} ,}$
(2.37)

where the fourth-order identity tensor ${\textstyle {\boldsymbol {\boldsymbol {I}}}}$ is defined as

 ${\displaystyle \left[{\boldsymbol {\boldsymbol {I}}}\right]_{\;ijkl}={\frac {1}{2}}\left(\delta _{ik}\;\delta _{jl}+\delta _{il}\;\delta _{jk}\right).}$
(2.38)

The last terms in 2.35 and 2.37 are not included in the definition of the tangent tensor at constitutive level. This term corresponds to the purely volumetric component of the tangent constitutive tensor and is already accounted for separately at element level in the implementation of the hybrid element (see section 2.2.2).

The deviatoric part of the tangent constitutive tensor ${\textstyle {\widetilde {\boldsymbol {\boldsymbol {C}}}}^{tan}}$ must be derived for each hyperelastic model according to the expression of its particular ${\textstyle {\widetilde {\Psi }}}$. Unfortunately, this is not always an easy task. Alternatively, the complete tangent constitutive tensor, or only its deviatoric part, can be computed numerically through the perturbation technique outlined by Miehe [55].

(1) The Clausius-Duhem inequality is a way of expressing the second law of thermodynamics in local form which is extensively used in continuum mechanics. This inequality is particularly useful in determining whether the constitutive relation of a material is thermodynamically allowable [53].

(2) The principle of material frame-indifference, also referred to in literature as material objectivity, requires that the constitutive equations not depend on the external frame of reference used to describe them or, in other words, the constitutive equations should be independent of the observer.

(3) To perform a time derivative of a variable in the current configuration, this variable must be first pulled-back to the reference configuration. The material time derivative is performed and, then, the variable is brought back to the current configuration through a push-forward. This is so because the derivative must be performed in an invariant frame. The Lie derivative is defined [19] as the change of a spatial field ${\displaystyle f\left(\mathbf {x} ,t\right)}$ relative to the vector field ${\displaystyle {\mathbf {v} ={d\mathbf {x} }/{dt}}}$: ${\displaystyle {{\mathcal {L_{\phi }}}\left(f\right)={\overrightarrow {\phi }}\left[{\frac {D}{Dt}}\left({\overleftarrow {\phi }}\left[f\right]\right)\right]}}$.

### 2.2.4 Characterization of soft biological tissue behaviour using hyperelasticity

The specific strain energy density functions most commonly used in modelling soft biological tissue behaviour are presented in this section. These correspond to the deviatoric part of the elastic potential ${\textstyle {\widetilde {\Psi }}}$ in 2.6, however, the tilde has been dropped in this section for simplicity. The expression for uniaxial stress is also given for each model and is then employed in a simple example to illustrate the characteristics of each model.

The deformation gradient tensor obtained when loading uniaxially in the ${\displaystyle x}$-direction is

 ${\displaystyle \left[{F}\right]=\left[{\begin{array}{ccc}\lambda \;\;&0&0\\0&\lambda ^{-1/2}&0\\0&0&\lambda ^{-1/2}\\\end{array}}\right],}$
(2.39)

where ${\textstyle \lambda }$ is the stretch in the loading direction, defined as ${\textstyle {\lambda =l/L_{0}}}$. Here, ${\textstyle L_{0}}$ is the original length and ${\textstyle l}$ is the deformed (stretched) length. Then, the second Piola-Kirchhoff stress in the loading direction is derived from the specific strain energy density function following

 ${\displaystyle S_{x}=2{\frac {\partial \Psi }{\partial I_{\mathbf {C} }^{(1)}}}{\Big (}1-\lambda ^{-3}{\Big )}+2{\frac {\partial \Psi }{\partial I_{\mathbf {C} }^{(2)}}}{\Big (}\lambda ^{-1}-\lambda ^{-4}{\Big )}}$
(2.40)

Note that ${\textstyle \Psi =\Psi {\big (}I_{\mathbf {C} }^{(1)},I_{\mathbf {C} }^{(2)}{\big )}}$ because the third invariant is ${\textstyle {I_{\mathbf {C} }^{(3)}=\det \left(\mathbf {C} \right)\approx {1}}}$ due to the near-incompressibility assumption. Here, the expressions ${\textstyle {I_{\mathbf {C} }^{(1)}=\lambda ^{2}+2\lambda ^{-1}}}$ and ${\textstyle I_{\mathbf {C} }^{(2)}=2\lambda +\lambda ^{-2}}$ have been used. Finally, the first Piola-Kirchhoff and Cauchy stresses in the loading direction are obtained through the transformations ${\textstyle {T=\lambda \,S}}$ and ${\textstyle {\sigma =\lambda ^{2}S}}$, respectively.

Neo-Hooke It is the simplest hyperelastic model, given by the specific strain energy density function [56]

 ${\displaystyle \Psi =C_{1}\left(I_{\mathbf {C} }^{(1)}-3\right),}$
(2.41)

where ${\textstyle C_{1}=\mu /2>0}$ is a material constant related to the initial shear modulus of the material, ${\textstyle \mu }$. Then,

 ${\displaystyle S=2\,C_{1}\left(1-\lambda ^{-3}\right)}$
(2.42)

Mooney-Rivlin The general form of the specific strain energy density function [57,58] for this model is

 ${\displaystyle \Psi =\sum \limits _{i+j=1}^{N}C_{ij}\left(I_{\mathbf {C} }^{(1)}-3\right)^{i}\left(I_{\mathbf {C} }^{(2)}-3\right)^{j},}$
(2.43)

where ${\textstyle C_{ij}}$ are the material constants and ${\textstyle N}$ is a positive integer that determines the order of the model. Then,

 ${\displaystyle S=2\left(C_{10}+C_{01}\lambda ^{-1}\right){\Big (}1-\lambda ^{-3}{\Big )}}$
(2.44)

is the corresponding second Piola-Kirchhoff stress in the loading direction of an uniaxial loading set-up for the first-order model (${\textstyle N=1}$). It is usually assumed that ${\textstyle C_{01}>0}$ and ${\textstyle C_{10}\leq {0}}$, yet ${\textstyle C_{10}}$ may be positive and still result in a positive-definite specific strain energy density function [59]. The first-order model with ${\textstyle C_{01}=0}$ becomes the neo-Hookean model. The stress for the second-order model (${\textstyle N=2}$) is

 ${\displaystyle {\begin{array}{l}S=2{\Big (}\,C_{10}+C_{01}\lambda ^{-1}+3\,C_{11}\left(\lambda ^{-2}-1\right)\left(1-\lambda \right)+\\+2\,C_{20}\left(\lambda ^{2}+2\lambda ^{-1}-3\right)2\lambda ^{-1}C_{02}\left(\lambda ^{-2}+2\lambda -3\right){\Big )}{\Big (}1-\lambda ^{-3}{\Big )}.\end{array}}}$
(2.45)

In general, higher order models are not as used since they require fitting a considerable amount of material constants to the experimental data.

Yeoh This model is defined by the specific strain energy density function [60]

 ${\displaystyle \Psi =\sum \limits _{i=1}^{3}C_{i}\left(I_{\mathbf {C} }^{(1)}-3\right)^{i},}$
(2.46)

where ${\textstyle C_{1}}$, ${\textstyle C_{2}}$ and ${\textstyle C_{3}}$ are material constants. The initial shear modulus is given by ${\textstyle \mu =2\,C_{1}>0}$. Then,

 ${\displaystyle S=2\left(C_{1}+2\,C_{2}\left(\lambda ^{2}+2\lambda ^{-1}-3\right)+3\,C_{3}\left(\lambda ^{2}+2\lambda ^{-1}-3\right)^{2}\right){\Big (}1-\lambda ^{-3}{\Big )}}$
(2.47)

Fung The specific strain energy density function of this model [3] is

 ${\displaystyle \Psi ={\frac {C_{1}}{2\,C_{2}}}{\Big (}\exp \left[{C_{2}\left(I_{\mathbf {C} }^{(1)}-3\right)}\right]-1{\Big )},}$
(2.48)

where ${\textstyle C_{1}}$ and ${\textstyle C_{2}}$ are the positive material constants related to the shear modulus and the stiffening, respectively. Then,

 ${\displaystyle S=C_{1}\exp \left[{C_{2}\left(\lambda ^{2}+2\lambda ^{-1}-3\right)}\right]{\Big (}1-\lambda ^{-3}{\Big )}}$
(2.49)

Veronda-Westmann This model adds a term based on the second invariant of ${\textstyle \mathbf {C} }$ to the exponential term of the Fung model that depends solely on the first invariant of this tensor. Thus, the specific strain energy density function [61] is given by

 ${\displaystyle \Psi =C_{1}{\Big (}\exp \left[{C_{3}\left(I_{\mathbf {C} }^{(1)}-3\right)}\right]-1{\Big )}+C_{2}\left(I_{\mathbf {C} }^{(2)}-3\right),}$
(2.50)

where ${\textstyle C_{1}}$, ${\textstyle C_{2}}$ and ${\textstyle C_{3}}$ are material constants. Then,

 ${\displaystyle S=2{\Big (}C_{1}C_{3}\exp \left[{C_{3}\left(\lambda ^{2}+2\lambda ^{-1}-3\right)}\right]+C_{2}\lambda ^{-1}{\Big )}{\Big (}1-\lambda ^{-3}{\Big )}}$
(2.51)

Ogden The specific strain energy density function [62] is based on the principal stretches, instead of the invariants of ${\textstyle \mathbf {C} }$,

 ${\displaystyle \Psi =\sum \limits _{i=1}^{N}{\frac {\mu _{i}}{\alpha _{i}}}{\big (}{\lambda }_{1}^{\alpha _{i}}+{\lambda }_{2}^{\alpha _{i}}+{\lambda }_{3}^{\alpha _{i}}-3{\big )}.}$
(2.52)

Here, the shear moduli ${\textstyle \mu _{i}}$ and the stiffening parameters ${\textstyle \alpha _{i}}$ must satisfy the consistency condition

 ${\displaystyle 2\mu =\sum \limits _{i=1}^{N}\mu _{i}\alpha _{i}\qquad \mathrm {with} \quad \mu _{i}\alpha _{i}>0\quad \mathrm {for\quad } i=\left\{1,2,...,N\right\}.}$
(2.53)

${\textstyle N}$ is a positive integer that determines the order of the model. The third-order model is the most widely used, although the first-order model (${\textstyle N=1}$) with ${\textstyle \alpha =2}$ reduces to the Neo-Hookean expression and the second-order model (${\textstyle N=2}$) with ${\textstyle \alpha _{1}=2}$ and ${\textstyle \alpha _{2}=-2}$ produces the Mooney-Rivlin expression. For the third-order model,

 ${\displaystyle S=\sum \limits _{i=1}^{N}{\frac {\mu _{i}}{\alpha _{i}}}\left({\lambda }^{\left(\alpha _{i}-2\right)}-{\lambda }^{-\left(\alpha _{i}/2+2\right)}\right)}$
(2.54)

As a general rule, a model with more parameters has a higher chance of producing a good fit to experimental data. The drawback is, of course, identifying the adequate parameter values. A model with a large amount of material parameters may be burdensome and time-consuming to fit. Therefore, the simplest model which gives a reasonable fit to the desired experimental data is usually selected for a given application. Experimental data that includes several types of test (uniaxial, biaxial, torsion, etc.) tends to require more complex models in order to correctly reproduce the tissue behaviour.

Figure 7 shows the stress-stretch responses for the hyperelastic models described above, whose material parameters have been fitted to experimental data of liver tissue under tensile and compressive loadings [63]. The values used for these parameters are given in Table 1. The fitting was performed manually through trial and error and does not pretend to be an exhaustive study but merely an illustrative example. The reader is referred to the comprehensive studies of hyperelastic material models for modelling biological tissue behaviour by Martins et al. [64] and Wex et al. [65] for more information on this subject.

 Figure 7: Uniaxial stress-stretch responses for the different hyperelastic models used to fit the experimental data [63] of liver tissue under tensile and compressive loadings. The material parameters in Table 1 have been used.
 Model Param. Value Neo-Hooke ${\displaystyle C_{1}}$ ${\displaystyle 100\,{\textrm {kPa}}}$ Mooney- ${\displaystyle C_{10}}$ ${\displaystyle 100\,{\textrm {kPa}}}$ Rivlin ${\displaystyle C_{01}}$ ${\displaystyle 30\,{\textrm {kPa}}}$ Mooney- ${\displaystyle C_{10}}$ ${\displaystyle 1\,{\textrm {kPa}}}$ Rivlin ${\displaystyle C_{01}}$ ${\displaystyle 0.03\,{\textrm {kPa}}}$ (${\displaystyle N=2}$) ${\displaystyle C_{11}}$ ${\displaystyle 2\,{\textrm {kPa}}}$ ${\displaystyle C_{20}}$ ${\displaystyle 32\,{\textrm {kPa}}}$ ${\displaystyle C_{02}}$ ${\displaystyle 20\,{\textrm {kPa}}}$ Yeoh ${\displaystyle C_{1}}$ ${\displaystyle 15\,{\textrm {kPa}}}$ ${\displaystyle C_{2}}$ ${\displaystyle 2.1\,{\textrm {kPa}}}$ ${\displaystyle C_{3}}$ ${\displaystyle 13\,{\textrm {kPa}}}$
 Model Param. Value Fung ${\displaystyle C_{1}}$ ${\displaystyle 22\,{\textrm {kPa}}}$ ${\displaystyle C_{2}}$ ${\displaystyle 1.36\,}$ Veronda- ${\displaystyle C_{1}}$ ${\displaystyle 1.2\,{\textrm {kPa}}}$ Westmann ${\displaystyle C_{2}}$ ${\displaystyle 32\,{\textrm {kPa}}}$ ${\displaystyle C_{3}}$ ${\displaystyle 1.98}$ Ogden ${\displaystyle \mu _{1}}$ ${\displaystyle 30\,{\textrm {kPa}}}$ ${\displaystyle \alpha _{1}}$ ${\displaystyle 1.4}$ ${\displaystyle \mu _{2}}$ ${\displaystyle 0.005\,{\textrm {kPa}}}$ ${\displaystyle \alpha _{2}}$ ${\displaystyle 16.2}$ ${\displaystyle \mu _{3}}$ ${\displaystyle -0.8\,{\textrm {kPa}}}$ ${\displaystyle \alpha _{3}}$ ${\displaystyle -12}$

The neo-Hookean and first-order Mooney-Rivlin models are clearly incapable of reproducing the J-shape of the stress-stretch curves characteristic of soft tissue behaviour. Nonetheless, they are typically used when the stretches of interest are in the initial linear interval of the tissue response due to their simplicity and ease in fitting. In addition, it is often used to represent the ECM behaviour in the context of mixing theory constitutive formulations.

The most popular hyperelastic models for representing soft tissue behaviour are the Fung and Yeoh models because they predict responses with reasonable accuracy and require a small amount of parameters. The Veronda-Westmann model produces comparable results but has not received as much attention in literature.

The third-order Ogden model generally produces accurate fits but requires identifying six material parameter values. Yet, its use in soft tissue modelling is also considerably extended owing to the fact that it is based on the principal stretches, which are directly measurable quantities.

Precisely for this reason, the Ogden hyperelastic model has been chosen for implementation in PLCd [1]. In addition, it can be easily reduced to the neo-Hookean and first-order Mooney-Rivlin models with adequate values of its material parameters. The neo-Hookean model has also been implemented in the code since it is the most basic and simple of hyperelastic models.

### 2.2.5 Neo-Hookean hyperelasticity

The neo-Hookean function 2.41 is introduced as the deviatoric term in the split specific strain energy density function 2.6, resulting in

 ${\displaystyle \Psi ={\widetilde {\Psi }}+\Psi _{vol}=C_{1}\left({\widetilde {I}}_{\mathbf {C} }^{(1)}-3\right)+{\frac {1}{2}}\kappa \left(J-1\right)^{2},}$
(2.55)

where the first expression in 2.16 has been taken into account for the volumetric part of the elastic potential. Here, the material constant ${\textstyle C_{1}}$ is related to the initial shear modulus through ${\textstyle C_{1}={\mu }/{2}}$, ${\textstyle \kappa }$ is the bulk modulus and ${\textstyle {\widetilde {I}}_{\mathbf {C} }^{(1)}=J^{(-2/3)}I_{\mathbf {C} }^{(1)}}$ is the modified volume-preserving first invariant of ${\textstyle \mathbf {C} }$.

As expected, if there is no deformation ${\textstyle \mathbf {C} =\mathbf {I} }$ and ${\textstyle J=1}$. Then, the elastic potential vanishes since ${\textstyle {\widetilde {I}}_{\mathbf {C} }^{(1)}=J^{(-2/3)}Tr\left(\mathbf {C} \right)=3}$ and ${\textstyle {\left(J-1\right)=0}}$.

Applying 2.28, the constitutive equation in the reference configuration

 ${\displaystyle \mathbf {S} =2\,C_{1}J^{-2/3}\left(\mathbf {I} -{\frac {1}{3}}I_{\mathbf {C} }^{(1)}\mathbf {C} ^{-1}\right)-p\,J\mathbf {C} ^{-1}}$
(2.56)

is obtained. Performing a push-forward operation 2.5 on this expression, the constitutive equation in the present configuration results in

 ${\displaystyle {\boldsymbol {\tau }}=2\,C_{1}J^{-2/3}\left(\mathbf {B} -{\frac {1}{3}}I_{\mathbf {B} }^{(1)}\mathbf {I} \right)-p\,J\mathbf {I} .}$
(2.57)

Considering the split definition 2.33 of the material elasticity tensor and the expression of its volumetric part 2.35 corresponding to the selected volumetric elastic potential, the material elasticity tensor is computed as

 ${\displaystyle {\begin{array}{l}{\boldsymbol {\boldsymbol {C}}}^{tan}={\frac {4}{3}}C_{1}J^{(-2/3)}{\Big (}\,{\frac {1}{3}}I_{\mathbf {C} }^{(1)}\mathbf {C} ^{-1}\otimes \mathbf {C} ^{-1}-\mathbf {I} \otimes \mathbf {C} ^{-1}-\mathbf {C} ^{-1}\otimes \mathbf {I} +I_{\mathbf {C} }^{(1)}{\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}{\Big )}\\[5mm]-p\,J\left(\mathbf {C} ^{-1}\otimes \mathbf {C} ^{-1}-2{\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}\right)+\kappa \,J^{2}\mathbf {C} ^{-1}\otimes \mathbf {C} ^{-1},\end{array}}}$
(2.58)

where ${\textstyle I_{\mathbf {C} }^{(3)}=J^{2}}$ has been introduced and the the fourth-order tensor ${\textstyle {\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}}$ is computed as in 2.36.

Then, the spatial elasticity tensor

 ${\displaystyle {\begin{array}{l}{\boldsymbol {\boldsymbol {c}}}^{tan}={\frac {4}{3}}C_{1}J^{(-2/3)}{\Big (}\,{\frac {1}{3}}I_{\mathbf {B} }^{(1)}\mathbf {I} \otimes \mathbf {I} -\mathbf {B} \otimes \mathbf {I} -\mathbf {I} \otimes \mathbf {B} +I_{\mathbf {B} }^{(1)}{\boldsymbol {\boldsymbol {I}}}{\Big )}\\[5mm]-p\,J\left(\,\mathbf {I} \otimes \mathbf {I} -2{\boldsymbol {\boldsymbol {I}}}\,\right)+\kappa \,J^{2}\mathbf {I} \otimes \mathbf {I} \end{array}}}$
(2.59)

is obtained by applying the push-forward operation 2.32 on 2.58. The fourth-order identity tensor ${\textstyle {\boldsymbol {\boldsymbol {I}}}}$ is defined in 2.38.

The tangent constitutive tensor in the reference and present configurations both satisfy the minor and major symmetries ${\textstyle {{\boldsymbol {C}}_{IJKL}^{tan}={\boldsymbol {C}}_{KLIJ}^{tan}={\boldsymbol {C}}_{JIKL}^{tan}={\boldsymbol {C}}_{IJLK}^{tan}}}$ and ${\textstyle {\boldsymbol {c}}_{ijkl}^{tan}={\boldsymbol {c}}_{klij}^{tan}={\boldsymbol {c}}_{jikl}^{tan}={\boldsymbol {c}}_{ijlk}^{tan}}$, respectively.

The terms containing the bulk modulus ${\textstyle \kappa }$ in 2.58 and 2.59 are not included in the definition of the tangent tensor at constitutive level when hybrid elements are used. This term corresponds to the purely volumetric component of the tangent constitutive tensor and is already accounted for separately at element level in the implementation of the hybrid element (see section 2.2.2). In addition, the tangent constitutive tensor in a pUL framework could be directly computed in the reference configuration instead of computing ${\textstyle {\boldsymbol {\boldsymbol {c}}}^{tan}}$ and then performing a pull-back to obtain ${\textstyle {\boldsymbol {\boldsymbol {C}}}^{tan}}$. In the pUL framework, the calculation of the stiffness tensor ${\textstyle \mathbf {K} _{UU}}$ is performed in the reference configuration and, therefore, requires ${\textstyle {\boldsymbol {\boldsymbol {C}}}^{tan}}$ instead of ${\textstyle {\boldsymbol {\boldsymbol {c}}}^{tan}}$ (see Figures 2 and 6).

The scheme at Gauss point level of the numerical integration in PLCd [1] of the isotropic neo-Hookean hyperelastic model is outlined in Table 2. The implementation has been validated by means of a simple example consisting in a single element subjected to homogeneous uniaxial tensile loading. An 8-noded hexahedral element with 1cm length sides and a single pressure point (Q1P0) is subjected to a displacement-driven pure tensile load applied in steps of 0.1mm (see Figure 8). The material constant ${\textstyle C_{1}}$ has been set to 27.2kPa and a penalizer value of ${\textstyle 10^{12}}$Pa has been considered for the bulk modulus ${\textstyle \kappa }$. The results for both TL and uPL formulations coincide with the reference curve from [64] and are shown in Figure 9.

 Algorithm at each load increment ${\displaystyle n}$ Given: deformation gradient tensor ${\displaystyle \mathbf {F} }$, elemental pressure ${\displaystyle p}$ and the material property ${\displaystyle C_{1}}$. If (TL framework) then Compute the right Cauchy-Green deformation tensor ${\textstyle \mathbf {C} =\mathbf {F} ^{T}\cdot \mathbf {F} }$ and its inverse ${\textstyle \mathbf {C} ^{-1}}$. Calculate the stress from the constitutive equation 2.56, ${\textstyle \mathbf {S} =2\,C_{1}J^{-2/3}{\Big (}\,\mathbf {I} -{\frac {1}{3}}Tr\left(\mathbf {C} \right)\mathbf {C} ^{-1}{\Big )}-pJ\mathbf {C} ^{-1}}$. Calculate the corresponding material elasticity tensor 2.58, ${\textstyle {\begin{array}{l}{\boldsymbol {\boldsymbol {C}}}^{tan}={\frac {4}{3}}C_{1}J^{(-2/3)}{\Big (}\,{\frac {1}{3}}Tr\left(\mathbf {C} \right)\mathbf {C} ^{-1}\otimes \mathbf {C} ^{-1}-\mathbf {I} \otimes \mathbf {C} ^{-1}-\mathbf {C} ^{-1}\otimes \mathbf {I} \\[1mm]+Tr\left(\mathbf {C} \right){\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}{\Big )}-p\,J\left(\mathbf {C} ^{-1}\otimes \mathbf {C} ^{-1}-2{\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}\right).\end{array}}}$ else (pUL framework) then Compute the left Cauchy-Green deformation tensor ${\textstyle \mathbf {B} =\mathbf {F} \cdot \mathbf {F} ^{T}}$. Calculate the stress from the constitutive equation 2.57, ${\textstyle {\boldsymbol {\tau }}=2\,C_{1}J^{-2/3}\left(\mathbf {B} -{\frac {1}{3}}I_{\mathbf {B} }^{(1)}\mathbf {I} \right)-p\,J\mathbf {I} }$. Calculate the corresponding spatial elasticity tensor 2.59, ${\textstyle {\begin{array}{l}{\boldsymbol {\boldsymbol {c}}}^{tan}={\frac {4}{3}}C_{1}J^{(-2/3)}{\Big (}\,{\frac {1}{3}}I_{\mathbf {B} }^{(1)}\mathbf {I} \otimes \mathbf {I} -\mathbf {B} \otimes \mathbf {I} -\mathbf {I} \otimes \mathbf {B} +I_{\mathbf {B} }^{(1)}{\boldsymbol {\boldsymbol {I}}}{\Big )}\\[1mm]-p\,J\left(\,\mathbf {I} \otimes \mathbf {I} -2{\boldsymbol {\boldsymbol {I}}}\,\right).\end{array}}}$ end
 Figure 8: Prescribed displacements applied on an 8-noded hexahedral linear element with a single pressure integration point (Q1P0) used in the homogeneous uniaxial tensile test example.
 Figure 9: Results for an 8-noded hexahedral linear element with a single pressure integration point (Q1P0) under homogeneous uniaxial tensile loading. Neo-Hookean hyperelastcity with ${\displaystyle C_{1}=27.2\,{\textrm {kPa}}}$ and a penalizer value ${\displaystyle {\kappa =10^{12}\,{\textrm {Pa}}}}$ has been considered. PLCd [1] results for both TL (blue crosses) and pUL (green squares) formulations coincide with the reference results from [64] (black line).

A second example reproducing a membrane with a hole at its centre subjected to tensile loading is performed to ensure the correct numerical implementation of the formulation. The membrane depicted in Figure 10 is subjected to the indicated displacement-driven loading ${\textstyle u}$. Due to the symmetry in the specimen, only a quarter of the membrane has been discretized using 360 hexahedral elements. Symmetry conditions are imposed, thus, nodes belonging to the symmetry ${\textstyle y-z}$ plane shown in Figure 10 have motion restricted in the ${\textstyle x-}$direction, while nodes belonging to the symmetry ${\textstyle x-z}$ plane have motion in the ${\textstyle y-}$direction restricted. Accumulative incremental displacements are imposed in the ${\textstyle y-}$direction on the nodes of the top part of the specimen, with the other directions left unrestrained. The material constant ${\textstyle C_{1}}$ has been set to 7.5kPa and a penalizer value of ${\textstyle 10^{12}}$Pa has been considered for the bulk modulus ${\textstyle \kappa }$.

The mechanical response for the TL formulation is illustrated in Figure 11 (top left) by means of the vertical reaction vs. stretch curve. The vertical reaction plotted is the total resultant reaction force in the ${\textstyle y-}$direction of the quarter of the specimen. Figure 11 (bottom) shows the pressure and the principal second Piola-Kirchhoff stress distribution for the final displacement value ${\textstyle u=75}$mm. The convergence curves for each load increment, plotted in Figure 11 (top right), show adequate convergence of the solution. A tolerance of of ${\textstyle 10^{-5}}$ has been used and only two iterations per load increment are required to reduce the maximum residue below this value.

The same example is repeated with quadratic hexahedral (Q2P0) and tetrahedral (T2P0) meshes (see Figures 12 and 13). The tetrahedral mesh contains elements whose side lengths are, on average, half the side length of the hexahedral element. Table 3 summarizes the characteristics of each mesh.

 Figure 10: Geometry (${\displaystyle r=100}$ mm, ${\displaystyle h=w=200}$ mm, ${\displaystyle t=20}$ mm) and loading of the membrane with a hole as described in [Waffenschmidt2014] (left); and hexahedral mesh and boundary conditions imposed on the quarter of the membrane that has been discretized (right).

 Element type Total number of elements Gauss points per element Total number of nodes Q1P0 360 8 630 Q2P0 360 27 2217 T2P0 39458 4 14739
 Figure 11: Membrane with a hole meshed with Q1P0 elements and subjected to tensile displacement-driven loading. Neo-Hookean hyperelasticity in a TL framework with ${\displaystyle C_{1}=7.5\,{\textrm {kPa}}}$ and ${\displaystyle \kappa =10^{12}\,{\textrm {Pa}}}$. Vertical reaction vs. stretch response (top left) and convergence curves of each load step (top right). Pressure ${\displaystyle p}$ (bottom left) and principal second Piola-Kirchhoff stress ${\displaystyle S_{1}}$ (bottom right) distributions at an imposed displacement value of ${\displaystyle u=75}$mm. Real deformation (${\displaystyle \times }$1) is plotted.
 Figure 12: Membrane with a hole meshed with Q2P0 elements and subjected to tensile displacement-driven loading. Neo-Hookean hyperelasticity in a TL framework with ${\displaystyle C_{1}=7.5\,{\textrm {kPa}}}$ and ${\displaystyle \kappa =10^{12}\,{\textrm {Pa}}}$. Vertical reaction vs. stretch response (top left) and convergence curves of each load step (top right). Pressure ${\displaystyle p}$ (bottom left) and principal second Piola-Kirchhoff stress ${\displaystyle S_{1}}$ (bottom right) distributions at an imposed displacement value of ${\displaystyle u=75}$mm. Real deformation (${\displaystyle \times }$1) is plotted.
 Figure 13: Membrane with a hole meshed with T2P0 elements and subjected to tensile displacement-driven loading. Neo-Hookean hyperelasticity in a TL framework with ${\displaystyle C_{1}=7.5\,{\textrm {kPa}}}$ and ${\displaystyle \kappa =10^{12}\,{\textrm {Pa}}}$. Vertical reaction vs. stretch response (top left) and convergence curves of each load step (top right). Pressure ${\displaystyle p}$ (bottom left) and principal second Piola-Kirchhoff stress ${\displaystyle S_{1}}$ (bottom right) distributions at an imposed displacement value of ${\displaystyle u=75}$mm. Real deformation (${\displaystyle \times }$1) is plotted.

### 2.2.6 Ogden hyperelasticity

The specific strain energy density function of the third-order Ogden hyperelastic material 2.52 is introduced as the deviatoric term in the split specific strain energy density function 2.6, resulting in

 ${\displaystyle \Psi ={\widetilde {\Psi }}+\Psi _{vol}=\sum \limits _{i=1}^{3}{\frac {\mu _{i}}{\alpha _{i}}}\left({\widetilde {\lambda }}_{1}^{\alpha _{i}}+{\widetilde {\lambda }}_{2}^{\alpha _{i}}+{\widetilde {\lambda }}_{3}^{\alpha _{i}}-3\right)+{\frac {1}{2}}\kappa \left(J-1\right)^{2}.}$
(2.60)

Here, the first expression in 2.16 has been taken into account for the volumetric part of the elastic potential. The shear moduli ${\textstyle \mu _{1}}$, ${\textstyle \mu _{2}}$ and ${\textstyle \mu _{3}}$, and the stiffness constants ${\textstyle \alpha _{1}}$, ${\textstyle \alpha _{2}}$ and ${\textstyle \alpha _{3}}$ are empirically determined material constants that must satisfy the consistency condition 2.53. The volume-invariant or deviatoric principal stretches ${\textstyle {\widetilde {\lambda }}_{1}}$, ${\textstyle {\widetilde {\lambda }}_{2}}$ and ${\textstyle {\widetilde {\lambda }}_{3}}$ are related to the principal stretches through ${\textstyle {{\widetilde {\lambda }}_{i}=J^{(-1/3)}\lambda _{i}}}$ with ${\textstyle {J=\lambda _{1}\lambda _{2}\lambda _{3}}}$ and ${\textstyle {{\widetilde {\lambda }}_{1}{\widetilde {\lambda }}_{2}{\widetilde {\lambda }}_{3}=1}}$.

Then, applying 2.28, the constitutive equation in the reference configuration

 ${\displaystyle \mathbf {S} =\sum \limits _{A=1}^{3}\beta _{A}\mathbf {M} _{A}-p\,J\mathbf {C} ^{-1}}$
(2.61)

is obtained. Here, the tensor ${\textstyle \mathbf {M} _{A}}$ is given by ${\textstyle \mathbf {M} _{A}=\lambda _{A}^{-2}\mathbf {N} _{A}\otimes \mathbf {N} _{A}}$, where ${\textstyle \mathbf {N} _{A}}$ is the eigenvector of the right Cauchy-Green deformation tensor such that ${\textstyle {\mathbf {C} =\sum _{A=1}^{3}\lambda _{A}^{2}\mathbf {N} _{A}\otimes \mathbf {N} _{A}}}$. In addition, the square of the principal stretches are the eigenvalues of ${\textstyle \mathbf {C} }$, i.e., ${\textstyle \mathbf {C} \cdot \mathbf {N} _{A}=\lambda _{i}^{2}\mathbf {N} _{A}}$. The scalar ${\textstyle \beta _{A}}$ is related to the deviatoric principal stretches through

 ${\displaystyle \beta _{A}=\sum \limits _{A=1}^{3}\mu _{i}\left({\widetilde {\lambda }}_{A}^{\alpha _{i}}-{\frac {1}{3}}\sum \limits _{p=1}^{3}{\widetilde {\lambda }}_{p}^{\alpha _{i}}\right)}$
(2.62)

Performing a push-forward operation 2.5 on 2.61 and considering that ${\textstyle {\boldsymbol {\tau }}=J\mathbf {\boldsymbol {\sigma }} }$ yields

 ${\displaystyle {\boldsymbol {\tau }}=\sum \limits _{A=1}^{3}\beta _{A}\mathbf {m} _{A}-p\,J\mathbf {I} ,}$
(2.63)

which corresponds to the constitutive equation in the present configuration. Here, ${\textstyle \beta _{A}}$ is defined as in 2.62 and the tensor ${\textstyle \mathbf {m} _{A}}$ is given by ${\textstyle \mathbf {m} _{A}=\mathbf {n} _{A}\otimes \mathbf {n} _{A}}$, where ${\textstyle \mathbf {n} _{A}}$ is the eigenvector of the left Cauchy-Green deformation tensor, i.e., ${\textstyle \mathbf {B} =\sum _{A=1}^{3}\lambda _{A}^{2}\mathbf {n} _{A}\otimes \mathbf {n} _{A}}$. Analogous to the reference configuration, the square of the principal stretches are the eigenvalues of ${\textstyle \mathbf {B} }$, i.e., ${\textstyle \mathbf {B} \cdot \mathbf {n} _{A}=\lambda _{i}^{2}\mathbf {n} _{A}}$.

Considering the split definition 2.33 of the material elasticity tensor and the expression of its volumetric part 2.35 corresponding to the selected volumetric elastic potential, the material elasticity tensor is computed as

 ${\displaystyle {\begin{array}{l}{\boldsymbol {\boldsymbol {C}}}^{tan}=\sum \limits _{A=1}^{3}\sum \limits _{B=1}^{3}\gamma _{AB}\,\mathbf {M} _{A}\otimes \mathbf {M} _{B}+2\sum \limits _{A=1}^{3}\beta _{A}{\dfrac {\partial \mathbf {M} _{A}}{\partial \mathbf {C} }}\\[5mm]-p\,J\left(\mathbf {C} ^{-1}\otimes \mathbf {C} ^{-1}-2{\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}\right)+\kappa \,J^{2}\mathbf {C} ^{-1}\otimes \mathbf {C} ^{-1}.\end{array}}}$
(2.64)

Here, ${\textstyle I_{\mathbf {C} }^{(3)}=J^{2}}$ has been introduced, the fourth-order tensor ${\textstyle {\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}}$ is computed as in 2.36 and the scalar ${\textstyle \gamma _{AB}}$ is related to the deviatoric principal stretches through

 ${\displaystyle \gamma _{AB}=\left\{{\begin{array}{ll}\sum \limits _{A=1}^{3}\mu _{i}\alpha _{i}\left({\frac {1}{3}}{\widetilde {\lambda }}_{A}^{\alpha _{i}}+{\frac {1}{9}}\sum \limits _{p=1}^{3}{\widetilde {\lambda }}_{p}^{\alpha _{i}}\right)&\mathrm {if} \quad A=B,\\[6mm]\sum \limits _{A=1}^{3}\mu _{i}\alpha _{i}\left(-{\frac {1}{3}}{\widetilde {\lambda }}_{A}^{\alpha _{i}}-{\frac {1}{3}}{\widetilde {\lambda }}_{B}^{\alpha _{i}}+{\frac {1}{9}}\sum \limits _{p=1}^{3}{\widetilde {\lambda }}_{p}^{\alpha _{i}}\right)&\mathrm {if} \quad A\neq B.\end{array}}\right.}$
(2.65)

The term ${\textstyle {\partial \mathbf {M} _{A}}/{\partial \mathbf {C} }}$ is given by

 ${\displaystyle {\begin{array}{cl}{\dfrac {\partial \mathbf {M} _{A}}{\partial \mathbf {C} }}&={\Big (}{\boldsymbol {\boldsymbol {I}}}-\mathbf {I} \otimes \mathbf {I} +J^{2}\lambda _{A}^{-2}\left(\mathbf {C} ^{-1}\otimes \mathbf {C} ^{-1}-{\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}\right){\Big )}{\big /}D_{A}\\[3mm]&+{\Big (}\lambda _{A}^{2}\left(\mathbf {I} \otimes \mathbf {M} _{A}+\mathbf {M} _{A}\otimes \mathbf {I} \right)-{\frac {1}{2}}{\dot {D}}_{A}\lambda _{A}\left(\mathbf {M} _{A}\otimes \mathbf {M} _{A}\right){\Big )}{\big /}D_{A}\\[3mm]&-J^{2}\lambda _{A}^{-2}\left(\mathbf {C} ^{-1}\otimes \mathbf {M} _{A}+\mathbf {M} _{A}\otimes \mathbf {C} ^{-1}\right){\big /}D_{A}.\end{array}}}$
(2.66)

where the fourth-order identity tensor ${\textstyle {\boldsymbol {\boldsymbol {I}}}}$ is defined in 2.38. Here, the scalar ${\textstyle D_{A}}$ is computed as ${\textstyle D_{A}=2\lambda _{A}^{4}-Tr\left(\mathbf {C} \right)\lambda _{A}^{2}+J^{2}\lambda _{A}^{2}}$ and its derivative is ${\textstyle {{\dot {D}}_{A}=8\lambda _{A}^{3}-2Tr\left(\mathbf {C} \right)\lambda _{A}-2J^{2}\lambda _{A}^{-3}}}$.

Then, the spatial elasticity tensor

 ${\displaystyle {\begin{array}{l}{\boldsymbol {\boldsymbol {c}}}^{tan}={\frac {1}{J}}\sum \limits _{A=1}^{3}\sum \limits _{B=1}^{3}\gamma _{AB}\,\mathbf {m} _{A}\otimes \mathbf {m} _{B}+{\frac {2}{J}}\sum \limits _{A=1}^{3}\beta _{A}{\dfrac {\partial \mathbf {m} _{A}}{\partial \mathbf {g} }}\\[5mm]-p\,J\left(\,\mathbf {I} \otimes \mathbf {I} -2{\boldsymbol {\boldsymbol {I}}}\,\right)+\kappa \,J^{2}\mathbf {I} \otimes \mathbf {I} \end{array}}}$
(2.67)

is obtained by applying the push-forward operation 2.32 on 2.64. The term ${\textstyle {\partial \mathbf {m} _{A}}/{\partial \mathbf {g} }}$ is given by

 ${\displaystyle {\begin{array}{cl}{\dfrac {\partial \mathbf {m} _{A}}{\partial \mathbf {g} }}&={\Big (}{\boldsymbol {\boldsymbol {I}}}_{\mathbf {B} }-\mathbf {B} \otimes \mathbf {B} +J^{2}\lambda _{A}^{-2}\left(\mathbf {I} \otimes \mathbf {I} -{\boldsymbol {\boldsymbol {I}}}\right){\Big )}{\big /}D_{A}\\[3mm]&+{\Big (}\lambda _{A}^{2}\left(\mathbf {B} \otimes \mathbf {m} _{A}+\mathbf {m} _{A}\otimes \mathbf {B} \right)-{\frac {1}{2}}{\dot {D}}_{A}\lambda _{A}\left(\mathbf {m} _{A}\otimes \mathbf {m} _{A}\right){\Big )}{\big /}D_{A}\\[3mm]&-J^{2}\lambda _{A}^{-2}\left(\mathbf {I} \otimes \mathbf {m} _{A}+\mathbf {m} _{A}\otimes \mathbf {I} \right){\big /}D_{A},\end{array}}}$
(2.68)

where the fourth-order tensor ${\textstyle {\boldsymbol {\boldsymbol {I}}}_{\mathbf {B} }}$ is defined as

 ${\displaystyle \left[{\boldsymbol {\boldsymbol {I}}}_{\mathbf {B} }\right]_{ijkl}={\frac {1}{2}}\left(\left[\mathbf {B} \right]_{ik}\left[\mathbf {B} \right]_{jl}+\left[\mathbf {B} \right]_{il}\left[\mathbf {B} \right]_{jk}\right).}$
(2.69)

The reader is referred to the work by Simo and Taylor [43] for further details on the derivation of the constitutive equation and tangent constitutive tensor from the Ogden hyperelastic specific strain energy density function.

As in the neo-Hookean model, the tangent constitutive tensor in the reference and present configurations both satisfy the minor and major symmetries. Again, the terms containing the bulk modulus ${\textstyle \kappa }$ in 2.64 and 2.67 are not included in the definition of the tangent tensor at constitutive level since they correspond to the purely volumetric component that is already accounted for separately at element level in the implementation of the hybrid element (see section 2.2.2).

The scheme at Gauss point level of the numerical integration in PLCd [1] of the isotropic third-order Ogden hyperelastic model is outlined in Table 4. The implementation has been validated by means of a simple example consisting in a single element subjected to homogeneous uniaxial tensile loading. An 8-noded hexahedral element with 1cm length sides and a single pressure point (Q1P0) is subjected to a displacement-driven pure tensile load applied in steps of 0.1mm (see Figure 8). The material constants have been set to ${\textstyle \mu _{1}=12.069\,{\textrm {Pa}}}$, ${\textstyle \mu _{2}=3.773\,{\textrm {MPa}}}$, ${\textstyle \mu _{3}=-52.171\,{\textrm {kPa}}}$, ${\textstyle \alpha _{1}=8.395}$, ${\textstyle \alpha _{2}=1.882}$ and ${\textstyle \alpha _{3}=-2.2453}$. A penalizer value of ${\textstyle 10^{12}}$ Pa has been considered for the bulk modulus ${\textstyle \kappa }$. The results for both TL and uPL formulations coincide with the reference curve from [66] and are shown in Figure 14.

 Algorithm at each load increment ${\displaystyle n}$ Given: deformation gradient tensor ${\displaystyle \mathbf {F} }$, elemental pressure ${\displaystyle p}$ and material properties ${\displaystyle \mu _{i}}$ and ${\displaystyle \alpha _{i}}$ for ${\displaystyle i=\left\{1,2,3\right\}}$. If (TL framework) then Compute the right Cauchy-Green deformation tensor ${\textstyle \mathbf {C} =\mathbf {F} ^{T}\cdot \mathbf {F} }$ and its inverse ${\textstyle \mathbf {C} ^{-1}}$. Calculate the stress from the constitutive equation 2.61, ${\textstyle \mathbf {S} =\sum \limits _{A=1}^{3}\beta _{A}\mathbf {M} _{A}-p\,J\mathbf {C} ^{-1}}$. Calculate the corresponding material elasticity tensor 2.64, ${\textstyle {\begin{array}{l}{\boldsymbol {\boldsymbol {C}}}^{tan}=\sum \limits _{A=1}^{3}\sum \limits _{B=1}^{3}\gamma _{AB}\,\mathbf {M} _{A}\otimes \mathbf {M} _{B}+2\sum \limits _{A=1}^{3}\beta _{A}{\frac {\partial \mathbf {M} _{A}}{\partial \mathbf {C} }}\\[1mm]+Tr\left(\mathbf {C} \right){\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}{\Big )}-p\,J\left(\mathbf {C} ^{-1}\otimes \mathbf {C} ^{-1}-2{\boldsymbol {\boldsymbol {I}}}_{\mathbf {C} ^{-1}}\right).\end{array}}}$ else (pUL framework)then Compute the left Cauchy-Green deformation tensor ${\textstyle \mathbf {B} =\mathbf {F} \cdot \mathbf {F} ^{T}}$. Calculate the stress from the constitutive equation 2.63, ${\textstyle {\boldsymbol {\tau }}=\sum \limits _{A=1}^{3}\beta _{A}\mathbf {m} _{A}-p\,J\mathbf {I} }$. Calculate the corresponding spatial elasticity tensor 2.67, ${\textstyle {\begin{array}{l}{\boldsymbol {\boldsymbol {c}}}^{tan}={\frac {1}{J}}\sum \limits _{A=1}^{3}\sum \limits _{B=1}^{3}\gamma _{AB}\,\mathbf {m} _{A}\otimes \mathbf {m} _{B}+{\frac {2}{J}}\sum \limits _{A=1}^{3}\beta _{A}{\frac {\partial \mathbf {m} _{A}}{\partial \mathbf {g} }}\\[1mm]-p\,J\left(\,\mathbf {I} \otimes \mathbf {I} -2{\boldsymbol {\boldsymbol {I}}}\,\right).\end{array}}}$ end
 Figure 14: Results for an 8-noded hexahedral linear element with a single pressure integration point (Q1P0) under homogeneous uniaxial tensile loading. Ogden hyperelasticity with ${\displaystyle \mu _{1}=12.069\,{\textrm {Pa}}}$, ${\displaystyle \mu _{2}=3.773\,{\textrm {MPa}}}$, ${\displaystyle \mu _{3}=-52.171\,{\textrm {kPa}}}$, ${\displaystyle \alpha _{1}=8.395}$, ${\displaystyle \alpha _{2}=1.882}$, ${\displaystyle \alpha _{3}=-2.2453}$ and a penalizer value ${\displaystyle \kappa =10^{12}}$. PLCd [1] results for both TL (blue crosses) and pUL (green squares) formulations coincide with the reference results from [66] (black line).

A second set of examples are performed to ensure the correct numerical implementation of the formulations. The membrane with a hole subjected to displacement-driven loading ${\textstyle u}$ (see Figure 10) in the previous section 2.2.5 is now modelled with Ogden hyperelasticity. The material constants considered are ${\textstyle \mu _{1}=0.04\,{\textrm {kPa}}}$, ${\textstyle \mu _{2}=3.7\,{\textrm {kPa}}}$, ${\textstyle \mu _{3}=-0.05\,{\textrm {kPa}}}$, ${\textstyle \alpha _{1}=6.4}$, ${\textstyle \alpha _{2}=1.9}$ and ${\textstyle \alpha _{3}=-4.2}$, in addition to a penalizer value ${\textstyle \kappa =10^{12}\,{\textrm {Pa}}}$. Again, the example has been repeated with three different meshes, each of which uses a different type of hybrid element (see Table 3).

The mechanical response of the TL formulation for the Q1P0 mesh is illustrated in Figure 15 (top left) by means of the vertical reaction vs. stretch curve. The vertical reaction plotted is the total resultant reaction force in the ${\textstyle y-}$direction of the quarter of the specimen. Figure 15 (bottom) shows the pressure distribution and the principal second Piola-Kirchhoff stress distribution for the final displacement value ${\textstyle u=200}$mm. The convergence curves for each load increment, plotted in Figure 15 (top right), show adequate convergence of the solution. A tolerance of ${\textstyle 10^{-5}}$ has been used and, again, only two iterations per load increment are required to reduce the maximum residue below this value. Figures 16 and 17 show these results for the Q2P0 and T2P0 meshes, respectively.

 Figure 15: Membrane with a hole meshed with Q1P0 elements and subjected to tensile displacement-driven loading. Ogden hyperelasticity in a TL framework with ${\displaystyle \mu _{1}=0.04\,{\textrm {kPa}}}$, ${\displaystyle \mu _{2}=3.7\,{\textrm {kPa}}}$, ${\displaystyle \mu _{3}=-0.05\,{\textrm {kPa}}}$, ${\displaystyle \alpha _{1}=6.4}$, ${\displaystyle \alpha _{2}=1.9}$ and ${\displaystyle \alpha _{3}=-4.2}$. A penalizer value ${\displaystyle \kappa =10^{12}\,{\textrm {Pa}}}$ has been considered. Vertical reaction vs. stretch response (top left) and convergence curves of each load step (top right). Pressure ${\displaystyle p}$ (bottom left) and principal second Piola-Kirchhoff stress ${\displaystyle S_{1}}$ (bottom right) distributions at an imposed displacement value of ${\displaystyle u=200}$ mm. Real deformation (${\displaystyle \times }$1) is plotted.
 Figure 16: Membrane with a hole meshed with Q2P0 elements and subjected to tensile displacement-driven loading. Ogden hyperelasticity in a TL framework with ${\displaystyle \mu _{1}=0.04\,{\textrm {kPa}}}$, ${\displaystyle \mu _{2}=3.7\,{\textrm {kPa}}}$, ${\displaystyle \mu _{3}=-0.05\,{\textrm {kPa}}}$, ${\displaystyle \alpha _{1}=6.4}$, ${\displaystyle \alpha _{2}=1.9}$ and ${\displaystyle \alpha _{3}=-4.2}$. A penalizer value ${\displaystyle \kappa =10^{12}\,{\textrm {Pa}}}$ has been considered. Vertical reaction vs. stretch response (top left) and convergence curves of each load step (top right). Pressure ${\displaystyle p}$ (bottom left) and principal second Piola-Kirchhoff stress ${\displaystyle S_{1}}$ (bottom right) distributions at an imposed displacement value of ${\displaystyle u=200}$ mm. Real deformation (${\displaystyle \times }$1) is plotted.
 Figure 17: Membrane with a hole meshed with T2P0 elements and subjected to tensile displacement-driven loading. Ogden hyperelasticity in a TL framework with ${\displaystyle \mu _{1}=0.04\,{\textrm {kPa}}}$, ${\displaystyle \mu _{2}=3.7\,{\textrm {kPa}}}$, ${\displaystyle \mu _{3}=-0.05\,{\textrm {kPa}}}$, ${\displaystyle \alpha _{1}=6.4}$, ${\displaystyle \alpha _{2}=1.9}$ and ${\displaystyle \alpha _{3}=-4.2}$. A penalizer value ${\displaystyle \kappa =10^{12}\,{\textrm {Pa}}}$ has been considered. Vertical reaction vs. stretch response (top left) and convergence curves of each load step (top right). Pressure ${\displaystyle p}$ (bottom left) and principal second Piola-Kirchhoff stress ${\displaystyle S_{1}}$ (bottom right) distributions at an imposed displacement value of ${\displaystyle u=200}$ mm. Real deformation (${\displaystyle \times }$1) is plotted.

### 2.2.7 Discussion

The neo-Hookean and Ogden hyperelastic models have been implemented successfully in PLCd [1] for both TL and uPL formulations, as confirmed by the validation examples provided. Both models produce practically the same results regardless of the configuration (reference or present one) chosen for the definition of the constitutive equation. These two definitions of a same model are, strictly speaking, different constitutive equations that reproduce a same behaviour.

The pUL formulation used requires a series of push-forwards and pull-backs that the TL one does not require. This obviously increases the overall calculation time and might induce numerical errors which could explain the slightly delayed convergence in the pUL examples reproducing complex stress states. For this reason, when possible, the TL formulation will be favoured over the pUL one. On occasions, a constitutive model, or the parameters defining a constitutive model, could be only available in the present configuration. In such case, the use of the pUL would be mandatory.

The analytical constitutive tensors have been implemented correctly since a low number of iterations per load step is required in all examples. In future models that might not allow for the analytical derivation of their constitutive tensor, the calculation by perturbations [55] will be required although this will certainly increase the computational cost.

Regarding the element order, lineal elements are known to be inadequate for the analysis of problems involving bending [49]. However, refining the mesh can improve substantially the quality of the results. Quadratic elements will always yield better results than their linear counterparts but, of course, are computationally more expensive. On the other hand, when the geometry allows it, hexahedral elements should be favoured over tetrahedral ones as they produce much better results using fewer elements and, hence, less computation time.

The fact that the hybrid elements used all have a single pressure point might translate into unrealistic stress distributions. A constant pressure per element results in an approximately constant volumetric stress per element, which means the stress is no longer continuous across elements. In addition, the incompressibility condition is enforced at elemental level and, thus, the individual Gauss points of a same element might have slight volume changes which even out at elemental level. The Cauchy and Kirchhoff stresses are obtained through a push-forward operation on the second Piola-Kirchhoff stress, which uses the deformation gradient tensor and its determinant. Then, these volume differences are reflected in the stress distribution at the present configuration. Applying a smoothing technique on the stress distribution in the post-processing with GiD [67] of the results reduces the visibility of this effect (see Figure 18), although the best solution would undoubtedly be to implement mixed elements with higher number of pressure integration nodes.

Nonetheless, the single pressure elements Q1P0, Q2P0 and T2P0 already implemented in the code will be used in the remainder of this study. The linear element Q1P0 will be favoured when possible owing to its computational efficiency. The aim of this study is to develop a constitutive formulation for representing the behaviour of biological tissues. Thus, the shortcomings of these elemental formulations are acknowledged, but it is not the purpose of this study to address them.

 Figure 18: Stress distributions at an imposed displacement value of ${\displaystyle u=75}$mm for the membrane with a hole modelled using neo-Hookean hyperelasticity in a TL framework and meshed with Q1P0 elements (above, see Figure 11) and Q2P0 elements (below, see Figure 12). Principal second Piola-Kirchhoff stress ${\displaystyle S_{1}}$ (left), smoothed principal second Piola-Kirchhoff stress ${\displaystyle S_{1}}$ (centre left), Kirchhoff stress ${\displaystyle \tau _{1}}$ (centre right) and smoothed Kirchhoff stress ${\displaystyle \tau _{1}}$ (right) distributions. Real deformation (${\displaystyle \times }$1) is plotted.

## 2.3 Finite-strain damage

A damaged material is characterised by a certain loss of load carrying capacity with respect to the original, undamaged material. This damage is a result of micro-voids and small fissures, present inside the matter, starting to grow. As they increase in size, there is a progressive material deterioration which can be measured through a loss in strength and stiffness [68]. This deterioration culminates in crack initiation, growth and final fracture. The latter are, obviously, occurring at larger scales than the material imperfections overlooked when dealing with continuum mechanics. Therefore, continuum concepts are somewhat difficult to introduce here. However, for the former phase, the process of progressive material deterioration leading up to the initial crack formation, a continuum approach can be followed to model the behaviour of the damaged material.

Continuum damage mechanics (CDM) is used to describe the progressive degradation experienced by the mechanical properties of materials prior to the initiation of macro-cracks. The small fissures and micro-voids occurring before are modelled as continuous, disregarding the discontinuity they introduce into the material properties but taking into account how they globally affect the value of these properties. The theories developed in the continuum damage framework are based on the thermodynamics of irreversible processes and use internal state variables [69]. They offer complementary possibilities to fracture mechanics, which deals with the actual fracture phenomena and requires the modelling of the cracks and voids present in the material.

The extent of damage in a given material is not directly measurable as strain is in elasticity and an alternative must be found in order to quantify the damage. The effective stress tensor ${\textstyle {\boldsymbol {\sigma }}_{0}}$ is introduced, which interprets the change in mechanical behaviour of a damaged material as a loss of effective area. Thus, a damage variable that is, in general, a tensor quantity is defined. Denoted by ${\textstyle {\boldsymbol {\mathbb {M} }}}$, this fourth-order tensor characterises the state of damage in an anisotropic model and transforms the “real” stress tensor, ${\textstyle {\boldsymbol {\sigma }}}$, into the effective stress tensor,

 ${\displaystyle {\boldsymbol {\sigma }}_{0}={\boldsymbol {\mathbb {M} }}^{-1}:{\boldsymbol {\sigma }}.}$
(2.70)

The concept of effective stress was first introduced in 1958 by Kachanov [70] through the use of a reduction factor associated with the amount of damage in the material. This factor would be equal to unity at the initial moment when no damage was present and reduce to zero, either at the moment of fracture localization or at the moment of rupture. In addition, a hypothesis of strain equivalence can be introduced. It states that “the strain associated with a damaged state under the applied stress is equivalent to the strain associated with its undamaged state under the effective stressâ [69], as illustrated in Figure 19.

 Figure 19: Schematic illustration of the effective stress concept, adapted and reproduced with permission from Oller [104].

In the case of isotropic damage, the mechanical behaviour of the small cracks and voids is completely independent of their orientation and affects the material properties in the exact same way, whichever material direction is considered. This, however, does not necessarily imply that the original undamaged constitutive tensor of the material is isotropic. The isotropic damage simply preserves the directional characteristics of the initial elastic tensor by degrading it equally in all directions and, therefore, is describable by a single scalar variable ${\textstyle d}$. In this model, the internal damage variable ${\textstyle {\boldsymbol {\mathbb {M} }}}$ is rewritten as ${\textstyle {\boldsymbol {\mathbb {M} }}=(1-d){\boldsymbol {\boldsymbol {I}}}}$, where ${\textstyle {\boldsymbol {\boldsymbol {I}}}}$ is the fourth-order symmetric identity tensor defined in 2.38. Then, 2.70 becomes

 ${\displaystyle {\bar {\boldsymbol {\sigma }}}={\frac {\boldsymbol {\sigma }}{\left(1-d\right)}},}$
(2.71)

where the damage parameter ${\textstyle d}$ is a measure of the loss of rigidity in the material and must be within the limits ${\textstyle d\in \left[0,1\right]}$. A value of ${\textstyle d=0}$ represents an undamaged material state whilst ${\textstyle d=1}$ represents a material state such that the material is completely degraded. This state can represent anything ranging from initial (macro)crack formation to local rupture, depending on what one defines as full material damage.

Since Kachanov first introduced the concept of effective stress, many authors have developed formulations based on this concept of elastic degradation to model damage in materials. Over the years, these formulations have consolidated and are now regarded as indisputable knowledge in the context of CDM [71,69,68,72]. This phenomenological approach is based on a rigorous mathematical and thermodynamic basis that will be reviewed in section 2.3.1. This formulation has proved to be a simple and effective tool in numerical modelling. Although initially formulated in an infinitesimal strain framework and as isotropic, it has been extended to include anisotropy [73,74], has been combined with plasticity [75,76,77] and viscoelasticity [78], and has been formulated for application to specific materials such as concrete [79,80], composites [81] or biological tissues [82], among others.

The first damage models developed in a finite strain context were proposed more than two decades ago, being the work of Simo [69] one of the most renowned. These are generally based on the additive decomposition of the specific strain energy density function (introduced for hyperelasticity in section 2.2.2) with damage affecting only the deviatoric part. Like the formulations by Miehe [83], de Souza [84] and other authors [34], these models were motivated by the Mullins effect (see section 2.1). More recently, damage models based on the decoupled volumetric-deviatoric response have been formulated to model the behaviour of fibred soft biological tissues [28,85,30,86,87]. The main characteristics of some of these models will be reviewed in section 2.3.2.

All these formulations use damage criteria and evolution laws which are defined to particularly suit the specific material behaviour being modelled. In this work, a generalized finite-strain damage softening model is proposed, which includes linear and exponential damage evolution laws that have been translated from an infinitesimal strain framework [9] into the present finite strain one. The novelty of this formulation is that, on the one hand, both proposed evolutions of the damage variable are based on solely two measurable material properties and, on the other hand, the formulation can be particularized for any decoupled volumetric-deviatoric hyperelastic constitutive model desired. Thus, the result is a general-purpose formulation which is versatile enough to model disparate material behaviours without requiring reformulation of the damage model or complex material parameter adjustments. The details of the proposed model are presented and discussed in sections 2.3.3 through 2.3.5.

### 2.3.1 Thermodynamic basis of finite-strain damage formulations

Damage models based on the split elastic potential definition with damage affecting solely the deviatoric part of the elastic potential associate the damage mechanism with the maximum distortional energy of the material and assume damage is independent of hydrostatic pressure [88]. Following this hypothesis [19], the specific strain energy density function is

 ${\displaystyle \Psi \left(\mathbf {C} ,\,D\,\right)=\Psi _{vol}\left(J\right)+\left(1-D\right){\widetilde {\Psi }}_{0}\left({\widetilde {\mathbf {C} }}\right),}$
(2.72)

where ${\textstyle {\widetilde {\Psi }}_{0}}$ is the undamaged isochoric or deviatoric part and ${\textstyle \Psi _{vol}}$ is its undamaged volumetric one, both given in the reference configuration. The Jacobian determinant ${\textstyle J}$ is related to the right Cauchy-Green deformation tensor, ${\textstyle \mathbf {C} }$, through ${\textstyle J=\left[\det \mathbf {C} \,\right]^{1/2}}$. The tilde in ${\textstyle {\widetilde {\mathbf {C} }}}$ indicates that it is the deviatoric or volume-preserving part of ${\textstyle \mathbf {C} }$, given by ${\textstyle {\widetilde {\mathbf {C} }}=J^{(-2/3)}\mathbf {C} }$. The functions chosen for ${\textstyle \Psi _{vol}}$ and ${\textstyle {\widetilde {\Psi }}_{0}}$ must be such that ${\textstyle \Psi _{vol}\left(J\right)=0}$ and ${\textstyle {\widetilde {\Psi }}_{0}({\widetilde {\mathbf {C} }})=0}$ hold if and only if ${\textstyle J=1}$ and ${\textstyle {\widetilde {\mathbf {C} }}={I}}$, respectively.

Expression 2.72 introduces an internal scalar damage variable ${\textstyle D\in \left[0,1\right]}$, which defines a reduction factor ${\textstyle \left(1-D\right)}$ similar to the one first proposed by Kachanov [70].

For an isothermal case with uniform temperature distribution and other standard arguments [88], the Clausius-Duhem inequality in the reference configuration 2.23 is reduced to

 ${\displaystyle -{\dot {\Psi }}+\mathbf {S} :{\frac {\dot {\mathbf {C} }}{2}}\geq 0,}$
(2.73)

where ${\textstyle \mathbf {S} }$ is the second Piola-Kirchhoff stress tensor. Considering ${\textstyle \Psi =\Psi \left(\mathbf {C} ,\,D\,\right)}$, the expression becomes

 ${\displaystyle -\left({\frac {\partial \Psi }{\partial D}}{\dot {D}}+2{\frac {\partial \Psi }{\partial \mathbf {C} }}:{\frac {\dot {\mathbf {C} }}{2}}\right)+\mathbf {S} :{\frac {\dot {\mathbf {C} }}{2}}\geq 0.}$
(2.74)

Then, introducing the specific strain energy density function defined in 2.72 and rearranging, the internal dissipation in the reference configuration

 ${\displaystyle \Xi ={\widetilde {\Psi }}_{0}{\dot {D}}+\left[\mathbf {S} -2\left({\frac {\partial {\Psi }_{vol}}{\partial \mathbf {C} }}+\left(1-D\,\right){\frac {\partial {\widetilde {\Psi }}_{0}}{\partial \mathbf {C} }}\right)\right]:{\frac {\dot {\mathbf {C} }}{2}}\geq {0}}$
(2.75)

is obtained. This inequality must hold true for any strain increment, therefore, the term in brackets must be null and the expression of the dissipation is reduced to

 ${\displaystyle \Xi ={\widetilde {\Psi }}_{0}{\dot {D}}\geq {0.}}$
(2.76)

Setting the term in brackets in 2.75 to zero yields

 ${\displaystyle \mathbf {S} =\mathbf {S} _{vol}+\left(1-D\,\right){\widetilde {\mathbf {S} }}_{0}\qquad {\textrm {with}}\quad \mathbf {S} _{vol}=-p\,J\mathbf {C} ^{-1}\quad {\textrm {and}}\quad {\widetilde {\mathbf {S} }}_{0}=2{\frac {\partial {\widetilde {\Psi }}_{0}}{\partial \mathbf {C} }},}$
(2.77)

which is the finite-strain version of the Kachanov effective stress concept. Here, the same volumetric specific strain energy density function used in hyperelasticity 2.28 has been introduced.

Consider now the material tangent constitutive tensor 2.33 and the expression for its volumetric part derived in 2.34. Then, a general expression for the material elastic-damage tangent constitutive tensor is derived as ${\textstyle {{\boldsymbol {\boldsymbol {C}}}^{tan}={\boldsymbol {\boldsymbol {C}}}_{vol}^{tan}+{\boldsymbol {\widetilde {\boldsymbol {C}}}}^{tan}}}$ with

 ${\displaystyle {\begin{array}{l}{\boldsymbol {\boldsymbol {C}}}_{vol}^{tan}=2{\dfrac {\partial \mathbf {S} _{vol}}{\partial \mathbf {C} }}=2p{\dfrac {\partial \left(J\mathbf {C} ^{-1}\right)}{\partial \mathbf {C} }}+2J\mathbf {C} ^{-1}\otimes p{\dfrac {\partial p}{\partial \mathbf {C} }}\quad {\textrm {and}}\\{\boldsymbol {\widetilde {\boldsymbol {C}}}}^{tan}=2{\dfrac {\partial }{\partial \mathbf {C} }}\left[\left(1-D\right){\widetilde {\mathbf {S} }}_{0}\right]=\left(1-D\right){\boldsymbol {\widetilde {\boldsymbol {C}}}}_{0}^{tan}-2{\dfrac {\partial D}{\partial \mathbf {C} }}\otimes {\widetilde {\mathbf {S} }}_{0}.\end{array}}}$
(2.78)

Now, the onset and evolution of the damage variable ${\textstyle D}$ in the constitutive equation 2.77 must be defined to complete the finite-strain damage model.

### 2.3.2 Characterization of soft biological tissue behaviour using finite-strain damage

Various models based on this approach have been proposed for characterizing the degradation or damage in biological tissues since the foundations of the finite-strain damage formulation were first developed in the 1980s.

 Figure 20: Schematic illustration of fibred soft tissue response to loading, inspired by [94]. (I) Toe region: ECM bears most of the loading, the tissue exhibits low value of practically constant stiffness. (II) Heel region: Fibres are progressively recruited, the stiffness of the tissue increases in a nonlinear manner. (III) Linear region: The majority of fibres have been recruited, the tissue has a high value of nearly constant stiffness. (IV) Progressive fibre failure: Fibres weaken and start to rupture, slightly decreasing the overall stiffness of the tissue. (V) Complete failure: The majority of fibres have ruptured, the tissue decreases drastically its load-bearing capacity until it completely fails.

Therefore, damage in soft tissue is related to the microstructural response of the tissue to loading. From a CDM viewpoint, the fibre rupture and ECM disruption can be assimilated to the small fissures and micro-voids observed in damaging inert materials and, hence, the justification for using CDM theory to model damage in soft tissue. Constitutive models that take into account the mechanisms causing damage at fibrilar level have been proposed [31,87]. Yet, a common approach from a phenomenological perspective is to describe fibre rupture and ECM disruption by means of adequate mathematical expressions of damage onset and evolution. These models mainly vary in how damage is defined to begin and evolve.

It must be noted that the hypothesis assumed here, namely, that damage affects solely the deviatoric part of the elastic potential, has been questioned by some researchers in the field. On the one hand, there is the debated issue of the quasi-incompressibility behaviour of soft tissues, attributed to the high volume fraction of water present in most soft tissues [11]. For supra-physiological loadings, the injury produced could potentially introduce changes in the water content of the tissue, which would require revisiting the quasi-incompressibility assumption and, hence, the use of a split elastic potential. On the other hand, accounting for the possibility of cavitational damage arising in soft tissues [90,91] would mean defining a damage model that exclusively affects the volumetric part of the material definition. For the purposes of the present study, damage will be considered to affect only the deviatoric part of the elastic potential, albeit the aforementioned limitations of this modelling approach.

Continuous and discontinuous damage variables Most authors define damage in terms of a discontinuous variable, ${\textstyle D=D\left(\alpha \right)}$, such that the reloading and the unloading response of the model coincide. Thus, damage accumulation never occurs for strain values below the previous maximum attained strain. A typical form of discontinuous damage variable ${\textstyle \alpha }$ is

 ${\displaystyle \alpha \left(t\right)={\underset {s\in \left[0,t\right]}{\max }}f[\,{\widetilde {\Psi }}_{0}\left(s\right)\,],}$
(2.79)

where ${\textstyle s\in \left[0,t\right]}$ is the history variable [83,84].

Miehe [83] introduced a continuous damage variable ${\textstyle \beta }$ into the definition of the total damage, ${\textstyle d=d\left(\alpha \right)+d\left(\beta \right)}$. Then, part of the damage accumulates continuously within the deformation process, even during unloading and reloading at strain values below the previous maximum. The continuous damage variable is defined in terms of the arc-length of the undamaged strain energy function as

 ${\displaystyle \beta \left(t\right)=\int _{_{0}}^{^{t}}\left|{\dot {\psi }}_{0}\left(s\right)\right|ds,}$
(2.80)

where the initial condition is ${\textstyle \beta \left(0\right)=0}$. Miehe worked in an exclusively spatial setting, but the formulation can be extrapolated to a TL framework. The continuous damage variable has been used for modelling the Mullins effect [92,93] and preconditioning [30,94] in biological tissues. Conversely, to account for the softening exclusively caused by damage (rupture) in soft tissues, researchers usually use a discontinuous damage variable.

Decoupled tensile-compressive models The damage model described in the previous section 2.3.1 has been developed based on the effective stress concept. In other words, the model is based on the hypothesis that damage is directly linked to the history of total strains and, in addition, damage is accumulative in nature. To illustrate this concept, consider the fact that, once the material starts to fissure and micro-voids appear, it is impossible that they later reduce their size and disappear. Thus, damage can grow but will never diminish. Nonetheless, it could happen that, if these small fissures and micro-voids are closed (due to a posterior compressive load state for example) they will not affect the macroscopic behaviour. Then, the damaged state is still present but it is considered to be passive, as opposed to the active state when the damage in the material is growing [68]. If the loading state were subsequently reversed and a tensile load applied again, the damage state would be the previously existing one before it would continue degrading. To account for this phenomena, models in which the specific strain energy density function [95] or the strain tensor [82] is decomposed into its tensile and compressive parts are used and damage is made to affect only the tensile one.

Onset and evolution of damage The evolution of the discontinuous damage variable is given by

 ${\displaystyle {\dot {D}}={\dot {\mu }}{\frac {\partial {\mathcal {F}}}{\partial \tau }}}$
(2.81)

where ${\textstyle {\dot {\mu }}}$ is a non-negative scalar named damage consistency parameter used to define the loading, unloading and reloading conditions through the Karush-Kuhn-Tucker complementary conditions

 ${\displaystyle {\dot {\mu }}\geq 0\quad ;\qquad {\mathcal {F}}\leq 0\quad {\textrm {and}}\qquad {\dot {\mu }}\;{\mathcal {F}}=0.}$
(2.82)

The damage surface is

 ${\displaystyle {\mathcal {F}}=G\left(\tau \right)-G\left(\tau ^{\max }\right)=0,}$
(2.83)

where ${\textstyle G\left(\tau \right)}$ is a damage evolution law given in terms of the norm ${\textstyle \tau }$, and ${\textstyle G\left(\tau ^{\max }\right)}$ is a scalar function of the damage threshold ${\textstyle \tau ^{\max }}$, which is the maximum reached value of ${\textstyle \tau }$ in the history of strains. The damage criterion proposed by Simo and Ju [69],

 ${\displaystyle \tau ={\sqrt {2{\widetilde {\Psi }}_{0}}},}$
(2.84)

is quite extended [28,96,35,85], although other functions of the specific strain energy density function are also used [97,98].

Numerous expressions have been proposed for the damage evolution law ${\textstyle G\left(\tau \right)}$ to characterize damage in soft tissue. Among the type of functions used are exponential [82,97,98,28,93,29,99], root [95], polynomial [96,100] and sigmoidal [101,29] ones. The reader is referred to the article by Peña [85] for a detailed review on the topic. Some of these functions are tailored to particular applications and their parameters can be assigned a physical meaning [27,99]. These functions tend to have a considerable amount of material parameters which are not always easily obtainable from experimental data. Other functions are more phenomenological in nature and require fitting fewer parameters [82,95], although they may not have a direct physical meaning. The functions proposed by Balzani and coworkers [98,30,102] and Calvo, Peña and coworkers [28,100,32,85] in general rely on only two or three material parameters which have been proven to show excellent fit with experimental data. In addition, some of the functions include viscosity, directional damage, separate fibre and matrix contributions and/or other softening effects such as preconditioning, permanent set and the Mullins effect.

In this work, for the purpose of setting the bases of a general constitutive formulation that represents the passive behaviour of soft tissues, the finite-strain extension of two evolution laws first described in an infinitesimal strain context are proposed. These are developed with the aim of reproducing a wide range of softening behaviours and are to be used in conjunction with mixing theory. Thus, they must be simple in their formulation, easy to fit to experiments and versatile enough to reproduce disparate soft tissue behaviour. The damage laws proposed use a discontinuous damage variable since they only aim to represent softening due to damage. Also, damage will be isotropic and affect equally the tensile and compressive states. Possible ways to account for anisotropy and different tensile/compressive behaviours in tissues will be addressed at composite level in section sec:RoM.

### 2.3.3 Proposed damage evolution laws

The linear and exponential explicit scalar functions described in [103,104] as damage evolution laws in an infinitesimal strain context have been translated to a finite strain framework to define ${\textstyle G\left(\tau \right)}$. A notable advantage of these laws is that they are based on only two material parameters with direct physical sense that can be experimentally determined.

Linear softening The damage variable ${\textstyle D}$ is defined as a scalar function with linear arguments

 ${\displaystyle D=G\left(\tau \right)={\frac {1-{\displaystyle {\tau _{0}^{d}}/{\tau }}}{1+H}},}$
(2.85)

where ${\textstyle \tau _{0}^{d}}$ and