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 finitestrain 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 homeostaticdriven 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.
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 patientspecific 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 wellestablished 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.
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 nearincompressibility 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 nearincompressibility 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 indepth 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.
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 inhouse code PLCd [1].
PLCd is an implicit FE code developed in Fortran and capable of solving finitestrain nonlinear threedimensional 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 multiscale 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 researchoriented 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 quasiincompressible hyperelastic and finitestrain 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 finitestrain 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 nonmechanistic 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 trialanderror approach, which is often tedious and generally timeconsuming. 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:
These constitutive formulations are developed in the framework of finitestrain continuum solid mechanics and are built on a wellestablished 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 postprocessor interfaces.
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 quasiincompressibility assumption are also described. Then, a generalized finitestrain 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 reversedamage 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 selfcontained 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.
The work included in this monography resulted in the following scientific publications:
In addition, part of the work was presented at the following conferences and workshops:
Finally, part of the work presented in this document is the result of collaborating with external researchers during the following research stays:
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 mid19th century [11]. The pioneering researchers in biomechanics observed that this behaviour is similar to that of rubberlike materials due to the longchain, crosslinked 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.
Continuumbased 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 finitestrain 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 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, nearlyincompressible and isotropic, which can be oversimplifying 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 CauchyGreen deformation tensor , or the principal stretches, which are the square roots of the eigenvalues of . 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 MooneyRivlin type, based on the first and second invariants of , 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 . Exponential and logarithmic forms are also widely used, for example in the Fung and VerondaWestmann models. The reader is referred to the comprehensive review on hyperelastic models for rubberlike 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 quasiincompressible 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 neoHookean and Ogden formulations will be developed. Their numerical implementation in the inhouse FE code PLCd [1] will be described and validated through numerical examples.
Although the hyperelastic models used to reproduce softtissue 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 , attributable to the behaviour of the extracellular matrix (ECM), and an anisotropic one , 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

(2.1) 
where is the amount of fibre families considered. Here, is defined in terms of the right CauchyGreen deformation tensor invariants and should account for the ECM behaviour. depends on both and the fibre orientation, typically through the pseudoinvariants of the fibre orientation vectors, also known as Spencer invariants [20]. The incorporation of the matrixfibre 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 is defined as an isotropic phenomenological hyperelastic function, while the definition of the anisotropic counterpart 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 rubberlike materials subjected to quasistatic cyclic tensile/compressive loading. The degradation is observed at strain levels below the maximum strain achieved in the history of deformation. This stresssoftening 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 rubberlike 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 finitestrain damage models developed up to date, with particular focus on those developed for biological applications. Then, a general finitestrain 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 distinct components (matrix and fibres), the second PiolaKirchhoff stress in the tissue will be given by

(2.2) 
Here, is the volumetric participation of each tissue component and is the second PiolaKirchhoff 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”.
Hyperelastic models are basically higherorder 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. Rubberlike 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.
In elastic materials, the constitutive behaviour is solely a function of the current strain level or deformation state and does not depend on the loading path followed or the history of strains suffered by the material. This implies that the stressstrain curves for loading and unloading are exactly the same and, thus, the original shape of the material is recovered upon unloading. In other words, the work stored during loading is retrieved during the unloading process, so there is no dissipation of internal energy. Hence, strains are reversible and rateindependent, and there is a biunivocal correspondence between stress and strain.
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]:
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 quasiincompressibility 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 quasiincompressible 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 neoHookean and Ogden models in sections 2.2.5 and 2.2.6, respectively. The details of the numerical implementation in the inhouse 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.
(^{1}) An objective stress rate is a time derivative of stress that is frame indifferent.
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 historydependent models such as damage, which will be addressed in section 2.3.
Then, the static equilibrium condition for an elemental discrete volume in the reference configuration is given by

(2.3) 
which must be satisfied for any displacement increment . Here,
The local equilibrium equations of each elemental volume are assembled in order to obtain the global equilibrium equation,

(2.4) 
where 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 stressstrain 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 and, thus, the nonlinearity is introduced in 2.4 through the stress tensor . Geometric nonlinearities appear when there are changes in the geometry which have a significant effect on the loaddeformation 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 displacementstrain compatibility tensor . 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 . 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, pushforward and pullback operators, denoted as and , respectively, are required. These are defined as

(2.5) 
where is the Kirchhoff stress tensor, is the Jacobian determinant and is the EulerAlmansi strain tensor.
Figure 2: Scheme of the total Lagrangian and partially updated Lagrangian formulations implemented in PLCd [1]. The subindex indicates iteration number in the present load increment. The definition of each term is available in the notation list. 
The nearincompressibility 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 quasiincompressible hyperelastic models, albeit the adequacy of the quasiincompressibility hypothesis of soft tissues has been debated in some cases [41,42].
Assuming soft tissue behaviour can be modelled by means of quasiincompressible 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 quasistatic loading case). Thus, the specific strain energy density function in the reference configuration is split into

(2.6) 
where corresponds to the volumepreserving or volumetric part and , to the isochoric or deviatoric one [43].
In addition, a quasiincompressible material has a quasiinfinite bulk elastic modulus, which results in an illconditioned stiffness matrix and locking problems when computing the solution using standard displacementbased 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:
In the present work, a twofield 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

(2.7) 
where is the displacement vector, is the pressure vector, are the stiffness matrices and and 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 , which results in

(2.8) 

(2.9) 
and

(2.10) 
Here and are the classical linear and nonlinear straindisplacement compatibility or transformation tensors, respectively [2]. is the tangent constitutive tensor, is the vector of pressure shape functions, is the bulk modulus and is a vector which in matrix form is the pressurerelated constitutive tensor . The internal forces are computed as

(2.11) 
and

(2.12) 
where is the pressure obtained from the displacement field and is the total element pressure obtained by independent interpolation, . Then, and are constitutive tensors relating strain and pressure, respectively, to stress, i.e.,

(2.13) 
The pressure constitutive tensor is given by

(2.14) 
where is the right CauchyGreen deformation tensor that is defined in terms of the deformation gradient tensor as . Since the volumetric part of the Cauchy stress tensor is given directly by the hydrostatic pressure ^{1}, , the pressure constitutive tensor is, in fact, a pullback operation 2.5 of this term. So, converts the Cauchy stress to Kirchhoff stress () and, through the inverse of the right CauchyGreen deformation tensor, the pullback operation is completed.
Now, the pressure obtained from the displacement field, , is defined positive in compression and computed as

(2.15) 
where can take different forms according to different authors [39,37,43,47], but is always defined in terms of the Jacobian determinant , which is a measure of the change of volume, and the bulk modulus . Possible functions of and their corresponding function for the pressure are

(2.16) 
In the first case, , the value of the bulk modulus directly dictates the value of the pressure when there is no change in volume, i.e., (see Figure 3). In addition, the tetrahedral elements that will be discussed later are very sensitive to the choice of the function . In this sense, the best choice is the one with the term 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

(2.17) 
where has been taken into account.
Figure 3: Pressure (positive in compression) vs. the Jacobian determinant for the different volumetric specific strain energy density functions given in 2.16 in terms of the bulk modulus . 
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 infsup 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 from the second line of 2.7 yields

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

(2.19) 
The formulation obtained is expressed solely in terms of displacement variables, and the new stiffness tensor and internal force vector differ from the classic displacement FE formulation in that they have the pressurerelated 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 infsup 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 infsup condition. Elements Q2P0 (bottom left) and T2P0 (bottom right) satisfy the infsup 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 infsup 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 infsup 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, , and the pressure shape function vector is reduced to . Then, the stiffness matrices 2.9 and 2.10, and internal force vector 2.12 become

(2.20) 

(2.21) 

(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 indicates iteration number in the present load increment. The definition of each term is available in the notation list. 
(^{1}) Here, is the secondorder identity tensor, which is defined as , being the Kronecker delta function.
(^{2}) The infsup 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].
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 ClausiusDuhem 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, , and temperature, , 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, .
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 . Then, the local form of the second law of thermodynamics in terms of internal specific energy density is expressed through the ClausiusDuhem inequality in reference configuration [54] as

(2.23) 
where is the specific entropy, is the heat flux vector and 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

(2.24) 
Then, the time derivative of the energy function is described in terms of the second PiolaKirchhoff stress tensor and its conjugate, the rate of GreenLagrange strain as . 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

(2.25) 
which shows that, as expected, derives from a strain energy function. It is often more convenient to use the right CauchyGreen deformation tensor instead of the GreenLagrange strain tensor. Introducing the relation , the stress results in

(2.26) 
This relation can be expressed using other stress tensors such as the first PiolaKirchhoff stress tensor, . If one recalls that, like and , the tensors and are work conjugates, then

(2.27) 
In any case, the elastic potential must fulfil:
Additionally, a reversible material (with no internal energy dissipation) must also satisfy:
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

(2.28) 
Here, the relation has been considered in the derivation of . This constitutive equation must now be completed by defining an appropriate deviatoric part of the specific strain energy density function, . 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 frameindifference ^{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 CauchyGreen deformation tensor,

(2.29) 
The invariants of are defined as
and the corresponding partial derivatives, which will prove useful in the derivation of the different hyperelastic constitutive equations, are
Expressions of isotropic specific strain energy density functions given in the present configuration are generally defined in terms of the invariants of the left CauchyGreen deformation tensor . This tensor, also known as Finger deformation tensor, is defined as and is known to have identical invariants to those of the right CauchyGreen deformation tensor: , and .
Finally, let us evaluate the rate of change of the constitutive equation 2.28, which can be expressed as

(2.30) 
Here, is a symmetric fourthorder 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 (), here the tangent constitutive tensor is not constant but a function of strain, 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 pushforward operation is performed on the time derivatives of stress and strain. This yields the Eulerian or spatial elasticity tensor, , which relates the Lie derivative ^{3} of the EulerAlmansi strain to the Lie derivative of the Kirchhoff stress,

(2.31) 
The spatial elasticity tensor is, in fact, the pushforward of the material one which, in index notation, is written as

(2.32) 
The tangent constitutive tensor is used to obtain the stiffness tensor 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

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

(2.34) 
where the operator is used to indicate the open product between the two secondorder terms. Further developing this expression in terms of the right CauchyGreen deformation tensor invariants produces

(2.35) 
Here, the fourthorder tensor is defined as

(2.36) 
Performing the pushforward operation 2.32 on 2.35, yields the volumetric part of the spatial elasticity tensor,

(2.37) 
where the fourthorder identity tensor is defined as

(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 must be derived for each hyperelastic model according to the expression of its particular . 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 ClausiusDuhem 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 frameindifference, 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 pulledback to the reference configuration. The material time derivative is performed and, then, the variable is brought back to the current configuration through a pushforward. 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 relative to the vector field : .
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 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 direction is

(2.39) 
where is the stretch in the loading direction, defined as . Here, is the original length and is the deformed (stretched) length. Then, the second PiolaKirchhoff stress in the loading direction is derived from the specific strain energy density function following

(2.40) 
Note that because the third invariant is due to the nearincompressibility assumption. Here, the expressions and have been used. Finally, the first PiolaKirchhoff and Cauchy stresses in the loading direction are obtained through the transformations and , respectively.
NeoHooke It is the simplest hyperelastic model, given by the specific strain energy density function [56]

(2.41) 
where is a material constant related to the initial shear modulus of the material, . Then,

(2.42) 
is the corresponding second PiolaKirchhoff stress in the loading direction of an uniaxial loading setup.
MooneyRivlin The general form of the specific strain energy density function [57,58] for this model is

(2.43) 
where are the material constants and is a positive integer that determines the order of the model. Then,

(2.44) 
is the corresponding second PiolaKirchhoff stress in the loading direction of an uniaxial loading setup for the firstorder model (). It is usually assumed that and , yet may be positive and still result in a positivedefinite specific strain energy density function [59]. The firstorder model with becomes the neoHookean model. The stress for the secondorder model () is

(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]

(2.46) 
where , and are material constants. The initial shear modulus is given by . Then,

(2.47) 
is the corresponding second PiolaKirchhoff stress in the loading direction of an uniaxial loading setup.
Fung The specific strain energy density function of this model [3] is

(2.48) 
where and are the positive material constants related to the shear modulus and the stiffening, respectively. Then,

(2.49) 
is the corresponding second PiolaKirchhoff stress in the loading direction of an uniaxial loading setup.
VerondaWestmann This model adds a term based on the second invariant of 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

(2.50) 
where , and are material constants. Then,

(2.51) 
is the corresponding second PiolaKirchhoff stress in the loading direction of an uniaxial loading setup.
Ogden The specific strain energy density function [62] is based on the principal stretches, instead of the invariants of ,

(2.52) 
Here, the shear moduli and the stiffening parameters must satisfy the consistency condition

(2.53) 
is a positive integer that determines the order of the model. The thirdorder model is the most widely used, although the firstorder model () with reduces to the NeoHookean expression and the secondorder model () with and produces the MooneyRivlin expression. For the thirdorder model,

(2.54) 
corresponds to the second PiolaKirchhoff stress in the loading direction of an uniaxial loading setup.
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 timeconsuming 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 stressstretch 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 stressstretch 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 
NeoHooke  
Mooney  
Rivlin  
Mooney  
Rivlin  
()  
Yeoh  
Model  Param.  Value 
Fung  
Veronda  
Westmann  
Ogden  
The neoHookean and firstorder MooneyRivlin models are clearly incapable of reproducing the Jshape of the stressstretch 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 VerondaWestmann model produces comparable results but has not received as much attention in literature.
The thirdorder 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 neoHookean and firstorder MooneyRivlin models with adequate values of its material parameters. The neoHookean model has also been implemented in the code since it is the most basic and simple of hyperelastic models.
The neoHookean function 2.41 is introduced as the deviatoric term in the split specific strain energy density function 2.6, resulting in

(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 is related to the initial shear modulus through , is the bulk modulus and is the modified volumepreserving first invariant of .
As expected, if there is no deformation and . Then, the elastic potential vanishes since and .
Applying 2.28, the constitutive equation in the reference configuration

(2.56) 
is obtained. Performing a pushforward operation 2.5 on this expression, the constitutive equation in the present configuration results in

(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

(2.58) 
where has been introduced and the the fourthorder tensor is computed as in 2.36.
Then, the spatial elasticity tensor

(2.59) 
is obtained by applying the pushforward operation 2.32 on 2.58. The fourthorder identity tensor is defined in 2.38.
The tangent constitutive tensor in the reference and present configurations both satisfy the minor and major symmetries and , respectively.
The terms containing the bulk modulus 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 and then performing a pullback to obtain . In the pUL framework, the calculation of the stiffness tensor is performed in the reference configuration and, therefore, requires instead of (see Figures 2 and 6).
The scheme at Gauss point level of the numerical integration in PLCd [1] of the isotropic neoHookean 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 8noded hexahedral element with 1cm length sides and a single pressure point (Q1P0) is subjected to a displacementdriven pure tensile load applied in steps of 0.1mm (see Figure 8). The material constant has been set to 27.2kPa and a penalizer value of Pa has been considered for the bulk modulus . 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
Given: deformation gradient tensor , elemental pressure and the material property .

Figure 8: Prescribed displacements applied on an 8noded hexahedral linear element with a single pressure integration point (Q1P0) used in the homogeneous uniaxial tensile test example. 
Figure 9: Results for an 8noded hexahedral linear element with a single pressure integration point (Q1P0) under homogeneous uniaxial tensile loading. NeoHookean hyperelastcity with and a penalizer value 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 displacementdriven loading . 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 plane shown in Figure 10 have motion restricted in the direction, while nodes belonging to the symmetry plane have motion in the direction restricted. Accumulative incremental displacements are imposed in the direction on the nodes of the top part of the specimen, with the other directions left unrestrained. The material constant has been set to 7.5kPa and a penalizer value of Pa has been considered for the bulk modulus .
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 direction of the quarter of the specimen. Figure 11 (bottom) shows the pressure and the principal second PiolaKirchhoff stress distribution for the final displacement value mm. The convergence curves for each load increment, plotted in Figure 11 (top right), show adequate convergence of the solution. A tolerance of of 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.
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 
The specific strain energy density function of the thirdorder Ogden hyperelastic material 2.52 is introduced as the deviatoric term in the split specific strain energy density function 2.6, resulting in

(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 , and , and the stiffness constants , and are empirically determined material constants that must satisfy the consistency condition 2.53. The volumeinvariant or deviatoric principal stretches , and are related to the principal stretches through with and .
Then, applying 2.28, the constitutive equation in the reference configuration

(2.61) 
is obtained. Here, the tensor is given by , where is the eigenvector of the right CauchyGreen deformation tensor such that . In addition, the square of the principal stretches are the eigenvalues of , i.e., . The scalar is related to the deviatoric principal stretches through

(2.62) 
Performing a pushforward operation 2.5 on 2.61 and considering that yields

(2.63) 
which corresponds to the constitutive equation in the present configuration. Here, is defined as in 2.62 and the tensor is given by , where is the eigenvector of the left CauchyGreen deformation tensor, i.e., . Analogous to the reference configuration, the square of the principal stretches are the eigenvalues of , i.e., .
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

(2.64) 
Here, has been introduced, the fourthorder tensor is computed as in 2.36 and the scalar is related to the deviatoric principal stretches through

(2.65) 
The term is given by

(2.66) 
where the fourthorder identity tensor is defined in 2.38. Here, the scalar is computed as and its derivative is .
Then, the spatial elasticity tensor

(2.67) 
is obtained by applying the pushforward operation 2.32 on 2.64. The term is given by

(2.68) 
where the fourthorder tensor is defined as

(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 neoHookean 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 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 thirdorder 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 8noded hexahedral element with 1cm length sides and a single pressure point (Q1P0) is subjected to a displacementdriven pure tensile load applied in steps of 0.1mm (see Figure 8). The material constants have been set to , , , , and . A penalizer value of Pa has been considered for the bulk modulus . 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
Given: deformation gradient tensor , elemental pressure and material properties and for .

Figure 14: Results for an 8noded hexahedral linear element with a single pressure integration point (Q1P0) under homogeneous uniaxial tensile loading. Ogden hyperelasticity with , , , , , and a penalizer value . 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 displacementdriven loading (see Figure 10) in the previous section 2.2.5 is now modelled with Ogden hyperelasticity. The material constants considered are , , , , and , in addition to a penalizer value . 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 direction of the quarter of the specimen. Figure 15 (bottom) shows the pressure distribution and the principal second PiolaKirchhoff stress distribution for the final displacement value mm. The convergence curves for each load increment, plotted in Figure 15 (top right), show adequate convergence of the solution. A tolerance of 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.
The neoHookean 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 pushforwards and pullbacks 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 pushforward operation on the second PiolaKirchhoff 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 postprocessing 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 mm for the membrane with a hole modelled using neoHookean hyperelasticity in a TL framework and meshed with Q1P0 elements (above, see Figure 11) and Q2P0 elements (below, see Figure 12). Principal second PiolaKirchhoff stress (left), smoothed principal second PiolaKirchhoff stress (centre left), Kirchhoff stress (centre right) and smoothed Kirchhoff stress (right) distributions. Real deformation (1) is plotted. 
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 microvoids 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 macrocracks. The small fissures and microvoids 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 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 , this fourthorder tensor characterises the state of damage in an anisotropic model and transforms the “real” stress tensor, , into the effective stress tensor,

(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 . In this model, the internal damage variable is rewritten as , where is the fourthorder symmetric identity tensor defined in 2.38. Then, 2.70 becomes

(2.71) 
where the damage parameter is a measure of the loss of rigidity in the material and must be within the limits . A value of represents an undamaged material state whilst 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 volumetricdeviatoric 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 finitestrain 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 volumetricdeviatoric hyperelastic constitutive model desired. Thus, the result is a generalpurpose 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.
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

(2.72) 
where is the undamaged isochoric or deviatoric part and is its undamaged volumetric one, both given in the reference configuration. The Jacobian determinant is related to the right CauchyGreen deformation tensor, , through . The tilde in indicates that it is the deviatoric or volumepreserving part of , given by . The functions chosen for and must be such that and hold if and only if and , respectively.
Expression 2.72 introduces an internal scalar damage variable , which defines a reduction factor similar to the one first proposed by Kachanov [70].
For an isothermal case with uniform temperature distribution and other standard arguments [88], the ClausiusDuhem inequality in the reference configuration 2.23 is reduced to

(2.73) 
where is the second PiolaKirchhoff stress tensor. Considering , the expression becomes

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

(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

(2.76) 
Setting the term in brackets in 2.75 to zero yields

(2.77) 
which is the finitestrain 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 elasticdamage tangent constitutive tensor is derived as with

(2.78) 
Now, the onset and evolution of the damage variable in the constitutive equation 2.77 must be defined to complete the finitestrain damage model.
Various models based on this approach have been proposed for characterizing the degradation or damage in biological tissues since the foundations of the finitestrain damage formulation were first developed in the 1980s.
Among the major known causes of physiological damage in biological tissues are the tensile and shearing structural failures due to the relative motions between tissue components [89]. Damage in soft tissues is caused by fibre tearing and rupture caused by excessively high stretches. Collagen fibres are, in fact, a bundle of collagen fibrils assembled together by proteoglycan crosslinks. As stretch on the fibre increases, the proteoglycan crosslinks begin to rupture, weakening the collagen fibre [31]. This results in the stressstretch response shown in Figure 20. Initially, the collagen fibres in the tissue are crimped and most of the loading is borne by the ECM (I: toe region). As loading progresses, the fibres begin to straighten and start to bear load. The fibres are said to be “recruited”. The high nonlinearity observed (II: heel region) and fast increase in stiffness value is due to more and more fibres being recruited. Once all fibres have been recruited, a constant stiffness region ensues (III: linear region) in which fibres are completely straightened and at full loadbearing capacity. At certain loading, fibres reach their maximum possible extension and start to rupture. At first, only a few fibres are disrupted and the tissue as a whole can still continue bearing load, albeit with a lower stiffness (IV: progressive fibre failure). Fibres continue to rupture as loading on the tissue increases, reaching a point where the tissue can no longer bear any loading. The tissue stiffness suddenly decreases just before complete rupture of the tissue (V: complete failure).
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 loadbearing 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 microvoids 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 quasiincompressibility behaviour of soft tissues, attributed to the high volume fraction of water present in most soft tissues [11]. For supraphysiological loadings, the injury produced could potentially introduce changes in the water content of the tissue, which would require revisiting the quasiincompressibility 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, , 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 is

(2.79) 
where is the history variable [83,84].
Miehe [83] introduced a continuous damage variable into the definition of the total damage, . 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 arclength of the undamaged strain energy function as

(2.80) 
where the initial condition is . 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 tensilecompressive 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 microvoids 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 microvoids 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

(2.81) 
where is a nonnegative scalar named damage consistency parameter used to define the loading, unloading and reloading conditions through the KarushKuhnTucker complementary conditions

(2.82) 
The damage surface is

(2.83) 
where is a damage evolution law given in terms of the norm , and is a scalar function of the damage threshold , which is the maximum reached value of in the history of strains. The damage criterion proposed by Simo and Ju [69],

(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 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 finitestrain 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.
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 . 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 is defined as a scalar function with linear arguments

(2.85) 
where and