Locking treatment of penalty-based gradient-enhanced damage formulation for failure of compressible and nearly incompressible hyperelastic materials

Soft materials are of major interest for biomechanics applications due to their high deformability and susceptibility to experience damage events under different loading scenarios. The present study is concerned with modelling damage evolution processes in these nonlinear materials whose structural responses are prone to locking when low-order kinematic interpolation is employed in the context of nonlinear Finite Element schemes. For this reason, a pair of gradient-enhanced continuum damage schemes are proposed with the aim of tackling mechanical failure problems in applications that exhibit shear and volumetric locking. In particular, we present the consistent formulation and the assessment of the corresponding performance of (i) a mixed displacement-enhanced assumed strain Q1Q1E24 employing a total Lagrangian formulation, and (ii) a three-field mixed displacement-pressure-Jacobian Q1Q1P0 formulation. The novel Q1Q1E24 and Q1Q1P0 formulations are consistently derived and numerically implemented, providing a satisfactory agreement with respect to ABAQUS built-in elements handling the treatment of shear and volumetric locking, respectively, in conjunction to the modelling damage phenomena via the use of a penalty-based gradient-enhanced formulation. This performance is examined via several numerical applications. Furthermore, the final example justifies the need for a formulation combining both mixed FE approaches to simulate problems encompassing both locking issues (shear and volumetric locking), which can be performed using a combination of the Q1Q1E24 and Q1Q1P0 herein proposed.

Among the yet-to-fully-tackle damage applications for gradient-enhanced CDM, deformable nearly incompressible hyperelastic materials constitute a recurrent topic due to their increasing range of applications, from rubber-like polymers to biomedical materials [47,48]. Being amply reported in the related literature, numerical modelling of nearly incompressible hyperelastic materials using displacement-based (single field) FE schemes suffer from the well-known volumetric locking, which overestimates the stiffness of materials when the Poisson ratio is high (ν ∼ 0.5) [49]. In engineering practice, the problem can be solved by the use of meshless techniques [50][51][52][53], smothered FE approaches (sFEM) [54][55][56], discontinuous Galerkin methods [57][58][59][60] or by lower integration algorithms of the divergence terms, among other alternatives. In particular, a mixed three-field displacementpressure-Jacobian mixed FE formulation which uses a lowerorder approximation for the Lagrange multipliers associated with the volumetric pressure [61][62][63][64][65] has been employed, being named Q1P0 or Simo-Taylor-Pister elements. This formulation, pioneered by the aforementioned authors [66] and enhanced by Simo and Taylor [67], and Miehe [68], is based on a local multiplicative split of the deformation gradient into deviatoric and volumetric parts, considering both the pressure and the Jacobian, which appear in the volumetric contribution, as primary unknowns, besides the displacement field. A disadvantage of the Q1P0 formulation is that it might provide unphysical solutions for compressible problems [69]. Being recently applied to model ductile damage [70] and fracture by means of a PF approach [71][72][73], we aim for the development of a Q1Q1P0 formulation with its combination with the gradient-enhanced CDM technique with the target of modelling damage in nearly incompressible hyperelastic materials avoiding volumetric locking issues. It should be mentioned that the proposed formulation for Q1Q1P0 element differs from the one utilized by ABAQUS hybrid elements, as the latter only employs a two-field displacementpressure formulation, being it employed to model CDM by Ostwald et al. [74]. Importantly, shear locking can also alter the compliance of the specimen with bending effects in large deformation analysis [75]. In order to overcome this shear locking pathology, several techniques have been proposed in the last decades: use of higher order Finite Element (FE) formulations [76][77][78], isogeometric analysis [79][80][81][82], reduced integration [83,84] and the application of mixed methods and formulations, such as the approach for incompatible modes [85] or the enhanced assumed strain technique [86][87][88], among others. The latter will be the one implemented within this approach.
The enhanced assumed strain (EAS), also named Q1E[•], being • the number of incompatible modes considered, is a method developed by Simo and Armero [87] for nonlinear continuum elements which, relying on the Hu-Washizu variational principle, enrich the deformation modes stemming from the single field displacement solution by means of several incompatible deformation modes at the element level. This method relies crucially on an additive decomposition of the deformation gradient into a conforming and an enhanced part which accounts for these incompatible modes. The resulting formulation corrects the over-stiffening of the structure by avoiding the phenomenon of shear locking. This technique has been widely applied to solid shells in order to suppress shear locking associated with bending modes [89][90][91][92][93][94]. This technique's major drawback lies in the instability associated with rank deficiency in the stiffness matrix which appears under compressive states [95], which remains an open question in the Computational Mechanics field. Having also been applied to non-local damage approaches such as PF frameworks [96][97][98] and very recently to CDM with reduced integration schemes [99], in this research, we aim at developing a full integration formulation combining the EAS method considering 24 incompatible deformation modes, Q1E24, with a gradient-enhanced CDM approach to analyze damage in samples under bending loads which are prone to display shear locking phenomena, i.e., Q1Q1E24.
In light of the previous discussion and after addressing potential locking events in CDM using non-local formulations, we present two non-local gradient-enhanced CDM approaches to solve locking issues. They are tailored within the geometrically nonlinear setting based on the work of Dimitrijević and Hackl [40], but by introducing an internal damage variable based on the model of Liebe et al. [100], in the line of the work by Waffenschmidt et al. [42]. Within the cited approach, the concept of enhancing the energy function via a gradient term of the independent damage variable is combined with a penalty parameter that simulates the equivalence between the local and non-local damage parameter, being this approach in line with the micromorphic gradienttype dissipative framework proposed by Forest [101]. This coupled two-equation (linear momentum and non-local damage balance) framework for large deformations is formulated in a weak form. The hyperelastic constitutive response is affected by the non-local damage scalar, which is approximated via an exponential law that triggers the deterioration of the structure when it overpasses a threshold.
Two displacement-continuum damage approaches will be built over this primary framework capturing damage in hyperelastic materials prone to shear and volumetric locking. For the former application, the EAS technique is included to encompass 24 incompatible deformation modes, implementing a formulation Q1Q1E24 that is suitable to model compressible samples subjected to bending. Then, separately, the mixed three-field Q1Q1P0 approach is formulated to tackle volumetric locking in nearly incompressible samples. The resulting coupled, highly nonlinear system of equations is solved via two Newton-Raphson type solution schemes: one local and one global employing a user-element subroutine UEL in the FE commercial software ABAQUS. In summary, we are presenting in this manuscript the first full-integration enhanced assumed strain (Q1Q1E24) and a novel mixed displacement-pressure-Jacobian (Q1Q1P0) schemes for gradient-enhanced CDM modelling and we have tested their performance by comparing them with the alreadyformulated standard CDM damage approach, i.e., Q1Q1.
This paper is structured as follows. The basic theory, which includes the constitutive behavior, the non-local gradient-enhanced damage formulation, and the thermodynamical postulates for the standard Q1Q1 element, is developed in Sect. 2. Section 3 displays the variational theorems for this Q1Q1 formulation and the two proposed ones, i.e., Q1Q1E24 and Q1Q1P0, providing a further insight on the numerical implementation in Sect. 4. To validate and test the potential of the proposed frameworks, a wide variety of numerical examples which consists of compressible and nearly incompressible large deformation problems prone to exhibit volumetric and shear locking have been addressed in Sect. 5. Some final remarks and conclusions are provided in Sect. 6.

Theoretical formulation
This section outlines the fundamental concepts and definitions of the current numerical framework addressing the use of gradient-enhanced for a standard CDM scheme, being specialized later for EAS and mixed u-p-J formulations. The proposed numerical methodology is specialized for hyperelastic material models.

Basic definitions and constitutive formulation at local level
Complying with standard nonlinear Continuum Mechanics, let an arbitrary spatial point defined in the current configuration be defined as x := ϕ(X, t), being ϕ(X, t) the nonlinear deformation map which projects the material points X from the initial configuration 0 ⊂ R n to the current one ⊂ R n . The transformation of the differential line elements throughout the deformation process is characterized by the deformation gradient F whose definition renders being 1 the second-order identity tensor and H(X, t), the material displacement gradient tensor. The Jacobian, i.e., the ratio of the deformed to the undeformed volume, being the determinant of F, shall fulfill the condition of J = det[F] > 0. In order to track the motion of the body from the material to the spatial configuration at time t, the displacement vector is defined as: Accordingly, the right and left Cauchy-Green tensors are obtained as, respectively: We postulate the existence of local free energy function . Without loss of generality, we consider a nonlinear compressible neo-Hookean constitutive law. This expression is plotted in Eq. (4) where μ and λ are the shear constant and λ = K − 2 3 μ, respectively, with K as the volumetric constant; and I 1 is the first invariant of the right Cauchy-Green tensor that is defined as I 1 := tr[C]. The particular form given in Eq. (4) also holds for the spatial configuration taking the left Cauchy-Green strain tensor as the main argument, i.e, loc (b) and therefore, From Eq. (4), the second Piola-Kirchhoff tensor S can be computed as follows Jacobians, can be computed directly from the derivation from the material description: where I sym C is a fourth-order tensor that has the following expression: , which employs the non-standard dyadic products. To obtain the spatial counterpart e, we perform the push-forward operation on Eq. (7) where I sym = [1⊗1 + 1⊗1]/2 denotes the fourth-order symmetric identity tensor. Based on Liebe et al. [100], we define a scalar damage function f d (κ), which recalling [42], f d (κ) should be at least twice differentiable, and tracks the material degradation relying on the evolution of a local variable κ ∈ [0, ∞], and whose evolution is ruled by the achievement of a threshold value κ > κ 0 in order to cause a loss in the stiffness of the structure. Therefore, we can state: These conditions guarantee two clearly differentiated states: f d = 1 identifies an intact stiffness at the spatial point level, whereas f d = 0 denotes a fully deteriorated stiffness state. This belongs to the formulation of the local damage model, whose linking procedure with the non-local damage framework will be envisaged in the subsequent sections.

Gradient-enhanced non-local formulation
In line with the approach proposed by Dimitrijević and Hackl [40], a regularized damage material response is achieved by the definition of a gradient-enhanced non-local function term nloc (φ, ∇ X φ, κ) in the reference configuration: This non-local contribution can be split into two separate terms: nloc grd (∇ X φ) containing the material gradient of the non-local damage field variable φ, which stands for the first term of a Taylor series expansion of φ at the material point; and nloc plty (φ, κ) is a penalty term which correlates the local damage variable κ with the non-local damage variable φ. The energy terms are specified as follows: where c d consists in a parameter that characterises the nonlocal character of the formulation; β d , a penalty parameter that enforces the local damage κ and non-local damage φ variables to be equivalent; and γ d , a switch parameter that is introduced to range between a local and non-local gradientenhanced model, respectively and the corresponding value is ranged like γ d ∈ {0, 1}. Consequently, the expression for the internal free energy function considering the previous non-local terms is given by With these expressions at hand, we define the material expressions for the damage-like vector field Y and the scalar damage-like variable Y : whose spatial values y and y are obtained via pushforwarding Eqs. (14)- (15)

Thermodynamic consistency
The thermodynamic consistency of the constitutive framework outlined above is examined through the exploitation of the Clausius-Plank inequality (local internal dissipation (D int ) inequality) [40], which under isothermal conditions is given by Following [42], we focus our attention on the corresponding damage-related terms leading to the definition of a reduced local dissipation D red : where we have introduced the thermodynamic force g as the derivative with respect to the local damage variable κ: We now define a thermodynamic force q ≥ 0, which is conjugated to the classical scalar damage variable d, as follows: Accordingly, the Clausius-Plank inequality holds when the reduced dissipation condition satisfies D red ≥ 0, if ∂ d κ > 0. In fact, q takes the interpretation of the energy release rate, consisting of the addition of a local and a non-local contribution that reads as Complying with these equations, we can now define the damage condition: where d < 0 stands for the purely elastic behavior and d = 0 notes a damaged state. According to Simo and Hughes [102], an optimization problem regarding a Lagrange multiplier λ can be proposed to represent the evolution of the damage variablė where κ d concerns the initial damage threshold. This last equation gives rise to the Karush-Kuhn-Tucker (KKT) conditions to model both the initiation and termination of damage.
which, can be expressed in a more detailed way as The continuous formulation for damage is completed with the definition of the damage function itself f d (κ), which follows an exponential-type law.
with η d > 0 standing for the exponential saturation parameter. It is worth mentioning that we have introduced a damage threshold parameter κ d that differs from the proposed approach by Dimitrijević and Hackl [40] and that was introduced in Waffenschmidt et al. [42], in order to avoid over-compensation of the damage curve due to the logarithmic expression in loc [Eq. (4)] that may lead to enhance the stress-strain curve, rather than weakening it.

Variational formulation of standard gradient-enhanced damage models: material and spatial formulations
The total potential energy of a system, , is obtained from the combination of an internal contribution int , which considers the action of internal forces, and an external contribution ext due to the addition of volume and surface forces, i.e., = int − ext .
Restricting the analysis of conservative loading cases, we can express the total potential of the system in the reference position of the arbitrary body under consideration as follows Since the problem is governed by the principle of minimum potential energy, the expression for the equation concerning the mechanical problem in the material configuration is obtained as that, can be particularized for each independent field as: Let V u = δu ∈ [H 1 ( 0 )] : δu = 0 on ∂ 0,u be the space of admissible displacement variations, and V φ = δφ ∈ [H 1 ( 0 )] : ∇ X δφ · N = 0 on ∂ 0 , the space of admissible test functions for non-local damage function. Furthermore, in the previous expression, the second Piola-Kirchhoff stress tensor renders S := 2 f d (κ)∂ C loc , and F V and T denote the body force and the traction vectors in the reference volume 0 and surface ∂ 0 , respectively.
The previous system of equations can be expressed in the spatial configuration by applying a standard push-forward operation: expressed in the current volume and surface ∂ and with σ identifying the Cauchy stress tensor that is accordingly affected by the degradation function f d (κ). Upon the use of the product rule and the divergence theory on Eqs.
Replicating the previous procedure for Eqs. (32)- (33) in the case of the current configuration, the Euler-Lagrange equations in the spatial form are given by being N and n both the normal vectors in the material and spatial configuration, respectively.

Variational formulation of gradient-enhanced damage models for enhanced assumed strain formulations
This section tailors the already established standard CDM model by combining it with the EAS method. Regarding this novel application, it is performed to alleviate shear locking pathologies in damage using low-order displacement interpolation in the subsequent finite element discretization scheme. We focus our development on the additive decomposition of the Green-Lagrange strain tensor into a displacement derived (E u ) and an enhancing counterpartẼ as follows [103]: E = E u +Ẽ. This differs from the alternative EAS scheme proposed by Simo and Armero [87] that accounts for the additive decomposition of the deformation gradient F = F u +F. Moreover, note that in the following derivation, the free-energy function is expressed in terms of the Green-Lagrange strain tensor. However, there exists a direct relation concerning the right Cauchy-Green tensor.
The point of departure of the formulation is based on the construction of the multi-field Hu-Washizu functional, where the displacement, the enhancing strain, the stress, and the non-local damage variable are the independent fields. This functional is given by where ext (u) identifies the external contribution due to the prescribed domain and boundary actions.
the space of the admissible enhancing strain; and V φ = δφ ∈ [H 1 ( 0 )] : ∇ X δφ · N = 0 on ∂ 0 , the space of admissible test functions for non-local damage function. The first variation of the total potential energy with respect to independent fields gives the following general expression: The previous expression can be particularized as follows from Eq. (42): Exploiting the orthogonality condition between the stress field S and the enhancing strain fieldẼ [87], the weak form of the coupled IBVP (Initial Boundary Value Problem) can be reduced to three independent fields, namely the displacement, the enhancing strain, and the non-local damage fields.
The weak form given in Eq. (44) (recalling S : In Eqs. (45)- (47), δ * int and δ * ext stand for the internal and external contributions of the generic field ( * ). In what follows, we turn our interest to the internal contribution of each independent field.

Variational formulation of gradient-enhanced damage models for penalty-based mixed formulations for nearly incompressible materials
The second presented methodology concerned in this investigation is the mixed Jacobian-pressure formulation originally proposed by Simo et al. [66]. In line with this work, we first perform a modification in the local strain energy density loc to a nearly incompressible neo-Hookean approach: loc vol (48) whereĪ 1 is the first invariant of the isochoric left Cauchy-Green tensorb and reads asĪ For the current three-field variational problem, we consider three fields as primary unknowns of the system {u,p,J }, where: •p is the Lagrange multiplier that accounts for the pressure •J is the dilatation, a constraint variable for the Jacobian of the material J (u).
For convenience, we express the total potential of the system in the current configuration: Again we denote: the space of admissible test functions for non-local damage function.
The first variation of the functional with respect to independent fields renders The previous expression can be expanded as follows.
The weak form of the coupled IBVP (Initial Boundary Value Problem) can be reduced to four independent fields, namely the displacement, the pressure, the dilatation and the non-local damage fields. It is given by with the isochoric contribution of the Cauchy stress σ iso being easily obtained from the derivation of the local strain energy density as In addition to this, it is worth highlighting the expressions for the Jacobians obtained by the volumetric and isochoric contributions: It is observed how the damage function term f d only multiplies the isochoric term. Therefore, for the forthcoming expressions in Sect. 4.1, when we refer to damaged stress and stiffness, they refer to the isochoric terms. The volumetric contribution is left unchanged.

Algorithmic treatment and finite element implementation details
This section outlines the description of the algorithmic description for the general gradient-enhanced damage formulation (Sect. 4.1), and subsequently, in Sect. 4.2, we describe the finite element implementation details describing the resulting operators and the interpolation schemes for each of the formulations given in Sects. 3.1-3.3.

Gradient-enhanced damage framework-algorithmic setting
This section details the algorithmic scheme within the context of an iterative and sequential solution scheme using nonlinear FE. In the sequel, we provide condensed information concerning the material and spatial formulations given in Sect. 3.1 in line with the salient results of Waffenschmidt et al. [42]. Recalling Eq. (24), this expression represents a nonlinear differential equation that should be numerically integrated within the time step interval t = t n+1 − t n ≥ 0 where t n+1 is the current time step and t n is the previous equilibrium solution of the system relying on a Newton-Rahpson-based solution of the corresponding nonlinear FE formulation. The backward Euler integration scheme of the damage variable κ at the current time step n + 1 renders where γ n+1 = tλ n+1 is the Lagrange multiplier at time t n+1 . Therefore, the incremental Karush-Kuhn-Tucker conditions take the form: The flux and source equations can be updated for a material description: where the superscript [•] und refers to undamaged variables. The spatial approach is obtained by just push-forwarding the magnitudes: where κ n+1 = κ n for an elastic incremental step. Complying with Eq. (23), for our incremental scheme, the incremental Lagrange multiplier γ n+1 is obtained by fulfilling the consistency equation: This nonlinear equation is solved by means of a Newton-Raphson (N-R) iterative scheme at the material point level: expanding Eq. (67) in a first-order Taylor series at γ k n+1 for a k-th N-R iteration, we obtain is the residual in the k-th iteration step and d γ r k n+1 is the Jacobian of this residual, which reads Now that we can calculate γ n+1 , the N-R scheme checks if the residual is below a pre-defined tolerance, and accordingly, the internal damage variable (Eq. (59)) and the flux and source terms for stresses and non-local damage magnitudes (Eqs. (61)-(66)) are updated. We now have all the ingredients for the global algorithm setting that is thoroughly summarised in Algorithm 1. For this, we compute the derivatives for all the source values required to complete the Jacobian formulation. For the sake of brevity, we consider f d (κ n+1 ) = f d in this series of equations. Complying with a material configuration approach: and applying a push-forward for the spatial approach: with

Finite element formulation and implementation details
This section addresses the FE derivation and the main implementation details of the proposed coupled system of nonlinear equations for each of the variational formulations proposed in Sect. 3. The baseline kinematic description for the displacement approximation complies with standard first-order 3-D 8-node hexahedral elements. The parametric space is defined as: The related literature has deeply reported the poor performance of this fundamental displacement formulation for bending-dominated applications and nearly incompressible elasticity, motivating the development of several mixed FE formulations. Algorithm 1 Algorithmic box for the local N-R scheme for gradient-enhanced damage constitutive model. MS = material scheme, SS = spatial scheme. 1: (4)), S und n+1 (MS, Eq. (5)), σ und n+1 (SS, Eq. (6)), E und n+1 (MS, Eq. (7)) and e und n+1 (SS, Eq. (8)) 3: Calculate driving force (79) 4: Comprobate damage function: Set κ n+1 = κ n . If Eq. (80) ≤ 0, go to Step 6. Otherwise, go to Step 5. 5: Compute local N-R to obtain incremental Lagrange multiplier γ n+1 iteratively a: Compute residual r k e: Update internal damage variable (83) and go back to Step 5a. For the sake of clarity, we specify the main differences between the three approaches herewith proposed: • A nonlinear CDM approach for the spatial configuration using the formulation of Sect

Discretisation scheme for the displacement and the non-local damage variable
Complying with standard isoparametric FEM, the reference and the current geometries can be interpolated using standard trilinear shape functions N I (N(ξ ) in matrix notation) as where X I and x I stands the nodal positions in the reference and the current configurations, respectively, and setting n n = 8 is the number of nodes. These nodal locations can be expressed in the corresponding vectors: X and x.
The interpolation of the displacements u and the non-local damage variable φ renders where d represents the nodal displacement vector, and φ represents the nodal values of the non-local damage variable; both defined at the element level. The material and spatial gradients of the shape functions N can be read as with ξ referring to the parametric coordinate system with coordinates ξ = {ξ, η, ζ }; and J e and j e as the material and spatial Jacobians of the isoparametric transformation, which allow the computation of the deformation gradient F as follows: With the previous definitions at hand, the corresponding material gradient quantities can be discretized as, for a material description, whereas the spatial gradients render

FE formulation of the gradient-enhanced damage model for spatial configuration-Q1Q1
The point of departure for the finite element formulation of the displacement-based gradient-enhanced damage model recalls the variational formalism defined in Eqs. (32)-(33), defining a coupled problem. The insertion of the interpolation scheme for u and φ leads to a discrete version of the residual forms denoted by R d and R φ that are defined as: For the application of Newton-type solution algorithms for the iterative solution of the boundary value problem, the linearization of the weak form is computed as follows: where * [•] denotes the directional derivative operator with respect to the field * .
Computing the derivatives of the residuals, we reach the Jacobian expressions that are required to solve the global N-R scheme: which leads to the following linearised system of equations that is solved by the global iterative N-R monolithic scheme

FE formulation of the gradient-enhanced damage model EAS-based elements-Q1Q1E24
For its numerical implementation and in line with the previous investigations for Enhanced Assumed Strain (EAS) mixed FE formulations, we herewith recall a material formulation defined in the reference configuration of the body. This formulation has been exploited for its usage in the modelling for solid shells [87,93,94,[96][97][98][104][105][106], as it is proven to block the appearance of shear locking in structures under bending configurations. Recalling from Sect. 3.2, the Cauchy-Green right tensor is computed as follows: The current definition of the enhancing part of the Green-Lagrange tensor relies on the formulation proposed by Andelfinger and Ramm [89] and Bischoff and Ramm [103], whose specific details are omitted for brevity reasons.
The enhanced strain field is defined at the element level via the matrix operator M(ξ ) where T -T 0 is the transpose of the inverse of a matrix that accounts for the terms of the Jacobian from the initial configuration J e that reads: Therefore, the incompatible strains ς are added into the FE implementation as an extra degree of freedom. The consistent linearization of this system is obtained through the Gateaux directional derivative concept, resulting in The linearised system of equations solved for the global N-R monolithic scheme reads ⎡ where the term for the residual of the incompatible strains R ς is given by and the tangent terms that form part of the Jacobian matrix read as where ∇ X N accounts for the nonlinear term that is added to the expression defined in Eq. (86), specific for material configuration schemes [75]. This is expressed for every node by: where N i, j is the derivative of the nodal shape function N i with respect to the j-coordinate and F i j accounts for the terms of the deformation gradient F. This expression is also included in his respective terms for the internal residuals and their associated Jacobians for the reference configuration.
Following the approach proposed by [93], since interelement continuity is not required for the enhanced strains, they can be removed as a DOF through a standard static condensation process, thus reaching the system of equations proposed in Eq. (98), having the element stiffness contributions likẽ and the newly defined residuals as

FE formulation of the gradient-enhanced damage model for spatial configuration employing a mixed FE formulation-Q1Q1P0
For the second novel approach that we have developed, the present formulation relies on the fundamental derivations by Simo and Hughes [66] and subsequently exploited by Miehe [68], whose effectiveness for modelling quasiincompressible materials (ν → 0.5) has been profusely assessed in the last two decades. In this concern, we recall a particular model where the primary unknowns are: (i) the displacement field u, (ii) the Lagrange multiplier for pressurẽ p, and (iii) the independent kinematic variableJ .
In line with the two previous approaches, we start by getting the residuals for these three primary unknowns by discretising from Eq. (51): The consistent linearization of this system is obtained through the Gateaux directional derivative concept, resulting in By deriving the residuals, we reach the expression for the Jacobian matrices: What is observed in this formulation is that in addition to the normal formulation, we have to introduce shape functions for the interpolations for the mixed variables related to the pressure Np and the dilatation NJ . However, since they do not have to satisfy the continuity between the elements, we can suppose that their values, Np and NJ , have a constant scalar value of 1.
As outlined above, this mixed finite element is herewith reformulated in order to accommodate gradient-enhanced damage models, taking Eq. (51) as the basis for its derivation. Therefore, we obtain that the baseline three-field mixed formulation is coupled with the non-local gradient-enhanced damage model, leading to a system of four residual equations.
As inter-element continuity is not required for both the pressure and the dilatation DOFs, they can be removed from the system of equations by employing a standard static condensation process, reaching the system proposed in Eq. (98), obtaining both the residual and stiffness contributions likẽ R d = R d + K K dp Rp − K dp (KpJ ) −1 RJ (130) where K = (KpJ ) −1 KJJ (KJp) −1 .  [107]. Furthermore, this instance does not consider the gradient-enhanced CDM approach, and we aim to validate the mixed u-p-J formulation.

Numerical examples
The next numerical experiment consists of a series of parametric studies carried out on a plate with a hole to study the influence of the regularisation properties of the model. We adopt the example in Sect. 5.2 of the work by Waffenschmidt et al. [42], replicating the dependence on the mechanical behavior for the compressible samples and analyzing the performance for nearly incompressible materials.
The last experiment aims to simulate a challenging example for a large deformation problem vulnerable to both shear and volumetric pathologies. For this, based on the proposed example in Sect. 4.1 in the work by Reese et al. [107], we propose a notched cylindrical shell subjected to an extreme bending load.

Benchmark example-nearly incompressible block under compression
Taken from Reese et al. [107], it is modelled the quarter of a cubic brick (symmetry conditions are considered on the planes X and Y) with 1,000 brick elements under a compressive status, i.e., a displacement of u z = −0.5 mm is applied on the area in bold of Fig. 1   In order to quantify the differences in performance among the three proposed frameworks, it is plotted the deformed configuration (Fig. 2) and the reaction force-displacement curves (Fig. 3) for every algorithm. Emphasizing the main advantages and disadvantages of them, it is itemized for every technique, in particular, the results obtained, also comparing with one sample run with ABAQUS C3D8H hybrid elements: • Q1 formulation: It can be envisaged that the performance here is not satisfying overall. Starting from the left image in Fig. 2, it is observed a high deformation on the outer surfaces of the block due to the over-stiffening of the structure caused by the large deformation process. Confirming this aspect with the force-displacement curve in Fig. 3, it is deduced that the main reason for this behavior is the volumetric locking that exhibits the solid by employing this algorithm, whose force-displacement curve highly overestimates the mechanical answer of the material, compared with the result obtained for the ABAQUS hybrid element. • Q1E24 Formulation: There were not found any results for this scheme due to convergence issues since the correction that the EAS method does on the global N-R scheme is considerably high in order to find an approximate solution for the software under incompressibility.
To tackle this issue, it was tried to run the simulation with different element sizes, but without any improvement in the convergence of the scheme, displaying the unsuitability of this formulation to model incompressible specimens. • Q1P0 Formulation: Without any doubt, the mixed u-p-J FE formulation proves to comply as the most robust performance in order to avoid the volumetric locking pathology. Its deformed configuration (right image, Furthermore, it shows a strong correlation in the mechanical behavior (see Fig. 3), yet considerably better performance in terms of convergence, with the hybrid two-field ABAQUS elements, as this latter shows, through midtesting, premature failure (Fig. 4). These results exhibit that the mixed u-p-J formulation has a worth-mentioning potentiality in quantifying the response of the incompressible sample and, at the same time, avoiding locking pathologies, unlike Q1 discretization and the ABAQUS C3D8H element, thus providing its solid capability in modelling nearly incompressible hyperelastic materials under large deformation states.

Test details
As for the second numerical example, the deformation of a plate with a hole under a traction load is considered. The geometry and boundary conditions for the full model are exhibited in Fig. 5. In line with the previous example, considering symmetry conditions on both X, Y, and Z planes, specifically, only one eighth part of the sample is modelled for the sake of reducing the computational cost. The bottommost surface is considered to be clamped, whereas the topmost one is subjected to a vertical displacement u z up until the post-peak behavior of the sample. The material properties employed for this approach are plotted in Table 1 for a Poisson ratio of ν = 0.25. Among the parameters for the non-local damage: • The β d parameter has been calibrated in order to ensure a solid convergence for the local N-R monolithic scheme. • γ d has the value of the unity to guarantee a non-local damage framework if c d > 0.
These compressible specimens will be run with algorithms encompassing both Q1Q1 and Q1Q1E24 formulations, not considering Q1Q1P0 as they are not suitable to model problems with ν < 0.45. The parameters related to the damage law i.e., both the damage threshold κ d and the saturation η d magnitudes along with the non-local regularisation parameter c d will be modified throughout these series of tests where the first aim will be the validation of the proposed CDM models with the constitutive behavior of the elements from ABAQUS, on one side, and the failure pattern, on the other

Quantitative and qualitative validation
First, in order to validate the CDM approaches, it is required to verify the mechanical performance of the proposed algorithms with the theoretical ABAQUS elements. To establish a quantitative standpoint with them, the damage parameters of both κ d and η d have to be adjusted properly. According to [42], higher values of η d do accelerate the onset of the damage process, being deviated from the purely neo-Hookean response plotted by the ABAQUS elements at lower stages of the corresponding loading process. On the contrary, really small saturation parameters approximate the curve to the theory but do delay the damage onset considerably after the deformation range considered. In addition to this, by considering a damage threshold of κ d = 0 MPa, i.e., the damage onset happens upon loading, not augmenting the stiffness of the mechanical curve. Thus, we have fixed a value of η d = 0.1 MPa −1 and κ d = 0 MPa for these validation tests with a discretization of 9,525 hexahedral elements and a regularisation parameter of c d = 1000 MPa −1 mm 2 . Having conducted all of them successfully, we display the reaction force-displacement curves for the samples with standard CDM (Q1Q1) and CDM + EAS (Q1Q1E24) formulations in comparison with the purely neo-Hookean response run with ABAQUS C3D8 elements and it is envisaged the same performance for each scheme, see Fig. 6. Both formulations display a solid equivalence with the ABAQUS elements' curve on the first stages of testing and, upon damage propagation, they both manage to capture the post-peak softening of the curve. Even though we get the same quantitative results, the most advantageous model for these configurations without shear locking is the Q1Q1 formulation, as the CDM + EAS model (Q1Q1E24) considers incompatible strains, which add up computational cost to the problem. Having tested the verification quantitatively with ABAQUS samples, for a qualitative address, we plot the evolution of crack propagation for three different displacements, see Figs. 7a-c. In order to get a clear damage pattern, we have calibrated the non-local damage parameters for a final failure around u z = 20 mm, employing for this: η d = 1 MPa −1 , κ d = 1 MPa and c d = 500 MPa −1 mm 2 . Following an expected behavior, the first isocontour at Fig. 7a reveals its nucleation near the notched region (a stress concentrator) only to be continued by a mode I propagation, see Fig. 7b. In the end, it is observed that upon reaching the end of the width of the plate, the crack grows in the direction of the height of the specimen, see Fig. 7c. Therefore, with these micro- graphs of a foreseeable failure pattern, we can also verify the proposed CDM approach qualitatively, obtaining similar results for Q1Q1 and Q1Q1E24 formulations. It is important to highlight that a tolerance in the function f d < 0.05 has been established to avoid ill-conditioning in the equations; the reader is referred to [42] for further information on this aspect.

Mesh objectivity
Another required calibration for CDM approaches concerns the verification for mesh independence, as the condition for mesh-objectivity is associated with the damage evolution being independent of the discretisation for high enough regularisation parameters c d . Therefore, we run a test for three different hexahedral meshes employing the non-local damage properties exhibited in Table 2: one coarse (810 elements), one medium-sized (9,595 elements) and one considerably refined (75,360 elements) mesh. The three simulations are compared by plotting the damage function isocontour f d for a displacement of u z = 8 mm in Fig. 8, and they turn out to be practically identical for both formulations Q1Q1 and Q1Q1E24, thus providing that the dissipation energy is independent of the size of the element.

Parametric study on the regularisation parameter c d
Having proven that the results reasonably adjust to the theoretical nonlinear elastic behavior and that the proposed frameworks are mesh-objective, with the medium-sized discretisation (9,525 hexahedral elements), we have carried out a parametric analysis by leaving the degree of regularisa-  Table 2 Non-local damage properties employed for the tests for verification of mesh objectivity tion parameter c d as a free variable. As in [42], in order to quantify the dependence of the mechanical behavior to this parameter, a parametric analysis is made with a range of c d = {10, 100, 500, 1000} MPa −1 mm 2 for the three different proposed formulations with the mechanical properties in Table 1. The employed parameters for the degradation law are η d = 1 MPa −1 and κ d = 1 MPa. Looking to prove the validity of both standard CDM and CDM + EAS approach to model this specimen, we exhibit the reaction-force displacement curves for values of c d = 1000 MPa −1 mm 2 , c d = 500 MPa −1 mm 2 and c d = 100 MPa −1 mm 2 plotted in Figs. 9a-c, respectively. In fact, we can establish that both the standard CDM approach referred to as the current configuration, and the combined CDM + EAS approach, referred to as the reference configuration, do manage to capture the full mechanical behavior of the specimen until the failure of the sample, even though this last one requires more computational cost due to the calculation of the incompatible strains.
Staying with the most computationally efficient formulation i.e., Q1Q1, we realize a further comparison for the results with different c d in Figs. 10 and 11. For the first series of images (Fig. 10), we have plotted several contour plots of the damage function f d for different c d and, nucleating from the notched region, the major difference is observed in the range of the affected region, being the gradient wider for lower values of c d , and the minimum value of f d , decreasing as the c d parameter is augmented. Furthermore, by plotting the reaction force-displacement (Fig. 11a) and the minimum f d value (Fig. 11b), we basically obtain that the failure of the specimen gets delayed by just increasing this regularisation parameter, increasing the maximum force and softening the collapse of the graph, being this result in line with what is exhibited in [42]. In addition to this, it is worth noting that the result for c d = 10 MPa −1 mm 2 falls short of reaching the peak behavior of the curve as it displays convergence issues, being this due that the mesh is not discretised enough for this low value of c d , meaning that for this mesh, the limit value of c d for total convergence falls between this value and c d = 100 MPa −1 mm 2 .

Test details
The next numerical example that we have considered is an incompressible version of the previous example, i.e., we increase the value of the bulk modulus K up to a value ∼ 10 7 MPa, associated with a Poisson ratio of ν = 0.49999. With only that subtle change, the same experiments that were carried out for the Sect. 5.2 are repeated, keeping the same geometry as Fig. 5 and the same material properties plotted in Table 1, except the aforementioned K . For these series of tests, the schemes employed consist in both Q1Q1 and Q1Q1P0 formulations, as the EAS technique has been checked not to be appropriate to model nearly incompressible specimens (see Sect. 5.1). In line with the previous section for compressible samples (Sect. 5.2), we start with the verification of the CDM models by comparing them with the hybrid elements from ABAQUS (Sect. 5.3.2). Upon verification, we repeat the parametric analysis of the previous sections where we study the dependence of c d on the mechanical behavior of the now nearly incompressible specimens (Sect. 5.3.3), where we end up elucidating which is the better formulation to conduct these numerical experiments.

Quantitative and qualitative validation
For the comparison between the standard (Q1Q1) and the three-field mixed CDM formulations (Q1Q1P0) with the referential two-field ABAQUS C3D8H element, required for computations with ν > 0.475, we plot the force- It is observed that although the graphs manage to capture the post-peak behavior, it is the Q1Q1P0 formulation the better performer for this problem, as at the early stages of the curve, it solidly matches the ABAQUS referential graph, unlike the Q1Q1 which slightly overestimates the trajectory due to volumetric locking pathologies caused by ν ∼ 0.5. The failure pattern for both formulations resembles the one plotted in Fig. 7, which are not shown here for the sake of brevity but do provide the qualitative check for the testing of nearly incompressible specimens.

Parametric analysis on the regularisation parameter c d
Considering the parameters for the degradation law to be η d =   Fig. 13. According to what was aforementioned in Sect. 5.3.2, standard CDM formulation provides an over-stiffened curve due to the incompressibility of the model that leads to the phenomenon of volumetric locking that does not address this scheme. For that reason, the most robust scheme to model this problem is deduced to be the Jacobian-pressure mixed framework. Even though that Q1Q1 theory provides a more prolonged softening of the curve during crack propagation, it is deduced from these results that the Q1Q1P0 element covers the volumetric locking pathology by considering both pressure and dilatation terms, see Eqs. (130)-(131), and that is proven by the reduction in the stiffness of the constitutive response of the sample shown in all the examples in Fig. 13.
Subsequently, it is exhibited the (a) reaction forcedisplacement curves and (b) the evolution of the minimum value of the damage function f d associated with every displacement in Fig. 14, where is envisaged an overall analogous pattern than the one plotted for the compressible cases conducted with the standard CDM formulation, see Fig. 11, i.e., both peak force and displacement are augmented monotonically with the increase of the regularisation parameter. In addition to this, the slope for the decrease of the minimum value of f d for changing displacements is also increased. Although not so robust on convergence as the Q1Q1 formulation for compressible specimens, with the quantitative and qualitative results provided in this section, Q1Q1P0 has provided to be a more solid formulation to model these incompressible experiments.

Test details
With the objective of addressing the further potential of the present frameworks (so far, the Q1Q1 and the Q1Q1E24 formulations have been proven to model large deformation problems for compressible specimens, while the Q1Q1P0 scheme have performed very robustly with volumetric locking pathologies in an incompressible problem), we aim to model the challenging problem of a cylindrical shell under a bending load. For this, a notched cylindrical specimen with the geometry and BCs exhibited in Fig. 15 is analyzed under extreme bending conditions to capture the capabilities of both EAS and three-field mixed FE formulations in modelling this large deformation problem prone to show locking pathologies. In line with the previous example, we have modelled both series of compressible and nearly incompressible mate-   Table 3 for specimens with a discretisation of 12,328 hexahedral elements. In order to avoid boundary effects, for the extremities of the cylindrical shell (where the fixed and the displacement conditions are applied), the damage saturation parameter η d has been increased, so the damage pattern is not affected by these phenomena. The final deformed configuration with the mesh for the bending experiments is displayed in Fig. 16, exhibiting the amount of deformation in the sample experiments.
In line with the plate tests, we have started this subsection by simulating this experiment to establish a comparison between the formulations suitable for compressible specimens modelling, i.e., Q1Q1 and Q1Q1E24, with the referential ABAQUS elements in Sect. 5.4.2, additionally presenting the failure isocontour for this problem.
Subsequently, a parametric analysis on addressing the damage saturation parameter η d is carried out in Sect. 5.4.3, where the main pursued aim is the comparison of both algorithms in displaying their potentiality to model this large deformation complex problem.

Quantitative and qualitative validation
To compare with a referential standpoint, we have conducted several tests with a damage saturation parameter η d = 0.1 MPa −1 with Q1Q1 and Q1Q1E24 formulations, along with the theoretical ABAQUS elements with (C3D8I) and without (C3D8) incompatible deformation modes.
The displayed force-displacement curves, exhibited in Fig. 17, do reveal that the standard CDM approach (Q1Q1) follows the trajectory described by the ABAQUS C3D8 element in the early stages. However, both elements are affected by parasitic shear locking, showing an over-stiffening of the curve caused by this pathology. Concerning the CDM + EAS formulation (Q1Q1E24), it is observed a solid correlation Fig. 17 Load-displacement curves for the cylindrical shell test conducted for Q1Q1 and Q1Q1E24 formulations compared with the ABAQUS theoretical curves with (C3D8I) and without considering incompatible modes (C3D8). The employed saturation parameter has been η d = 0.1 MPa −1 between the mechanical performance of this element with the curve of the ABAQUS C3D8I element, which is the one required to model this kind of problems prone to locking i.e., under bending loads. It is worth highlighting that the match of these two curves is not perfect due to the EAS contribution being regularised via f d , affecting the terms in Eqs. (108)-(111), as these are multiplied by the function damage, meaning that the correction gets reduced as soon as the material gets damaged, being the EAS contribution eliminated upon total damage. This regularisation is done in order to avoid the convergence issues that arise by using this formulation, as the correction in the stiffness curve grows considerably throughout the experiment. However, this aspect does not undermine the potentiality of the CDM + EAS approach to model compressible materials under shear locking, as the correlation is very robust with the referential curve.
The damage evolution of the specimen during the test is displayed in Fig. 18. The crack onset occurs in the center of  (Fig. 18a), being propagated afterwards in Mode I in the Z-direction, first to the extremities of the specimen (Fig. 18b) and then, in the direction towards the notch (Fig. 18c). These isocontours have been obtained from a test conducted by employing the for the Q1Q1E24 formulation with a damage saturation parameter of η d = 0.5 MPa −1 , demonstrating to be consistent with analysis until failure of a specimen under extreme bending conditions.

Parametric analysis on the saturation parameter Á d
This section is focused on analyzing the role of the saturation parameter in the mechanical performance of the specimen. Varying this η d parameter in a range of η d = {0.5, 1, 2} MPa −1 for both the Q1Q1 and the Q1Q1E24 formulations. The differences among them for η d = 0.5 MPa −1 (Fig. 19a), η d = 1 MPa −1 (Fig. 19b) and η d = 2 MPa −1 (Fig. 19c) are exhibited by means of the force-displacement curves, where it can be deduced that the CDM + EAS formulation (Q1Q1E24) is verified to be a very effective tool to conduct the problem of extreme bending for compressible notched cylindrical shells, as is the better performer to tackle the shear locking phenomenon that causes an over-stiffening in the curves.
Intending to plot the η d dependence on one single image, we put together the reaction force-displacement curves for the Q1Q1E24 formulation, along with the evolution of the minimum value of f d in Fig. 20. It is deduced that the role of this parameter is similar to the one of c d , progressively delaying the failure of the sample along with augmenting the peak force when this saturation magnitude is increased.

Test details
The last series of experiments consists of the modelling of the previous examples in Sect. 5.4, now for nearly incompressible specimens, carried out in order to test the validity of the mixed u-p-J CDM formulation (Q1Q1P0) to model this tricky problem. Employing the same configuration as in Fig. 15 and the same properties as in Table 3, but this time, changing the bulk modulus K up to a value of ∼ 10 7 MPa, associated with an established Poisson ratio of ν = 0.49999 and in line with what has been realized in the previous sections, we start by comparing the results of Q1Q1 and Q1Q1P0 schemes with the same tests conducted with ABAQUS elements in Sect. 5.5.2. Once this demonstration has been fulfilled, a last parametric analysis addressing the dependence on the saturation parameter η d is conducted in Sect. 5.5.3.

Quantitative and qualitative validation
Displaying the same failure pattern as the one exhibited for compressible specimens (see Fig. 18, which is not included here for the sake of brevity), we focus now on a quantitative viewpoint with the two proposed formulations along with the referential ABAQUS elements, see Fig. 21. First, we can definitely conclude that Q1Q1 formulation is not suitable to model this problem as both the volumetric and shear locking cause a very abrupt rise of the load of the sample, justifying the invalidation of this formulation.
Concerning the mixed u-p-J CDM (Q1Q1P0) formulation, it can be observed that it covers the volumetric locking solidly, as the curve is below the one related to the behavior of the hybrid element ABAQUS C3D8H. However, it falls short of covering the shear locking, as this graph is considerably less compliant than the hybrid ABAQUS element that considers incompatible deformation modes (C3D8IH), i.e., the one required to run this simulation. Observing how robustly did perform the Q1Q1E24 formulation in the compressible cases analyzed in Sect. 5.4, the extension of the present formulation of Q1Q1P0 to Q1Q1P0E24 is proposed to model these specimens and will be addressed in future work.

Parametric analysis on the saturation parameter
Considering the same range for the damage saturation parameter η d = {0.5, 1, 2} MPa −1 , at last, we plot in a single curve all the load-displacement curves for different η d conducted by Q1Q1P0 (see Fig. 22a), along with the minimum value of the degradation law f d for every displacement (see Fig. 22b). What is envisaged in this last representation is an analogy of what was exhibited before for the compressible specimens: with the reduction in η d , the softening in the quantitative response associated with the mechanical performance is delayed. Therefore, even though there is no solid correlation with the theory, the Q1Q1P0 formulation has proven to be the one among the three conducted schemes to proportionate conclusive results for the large deformation bending tests applied on cylindrical shells.

Conclusions
Within this article, a duo of gradient-enhanced continuum damage formulations has been developed to model failure in compressible and incompressible isotropic hyperelastic specimens, with the main focus of addressing the volumetric and shear locking pathologies. Among these experiments, the samples have been tested to a wide range of different types of loading, including traction, compression, and bending, comparing the two novel schemes, Q1Q1E24 and Q1Q1P0, with the Q1Q1 referential CDM formulation.
Making a one-by-one analysis of the performance of the aforementioned schemes, first, the spatial standard CDM framework (Q1Q1), based on the one proposed by [42], has been validated in this work to be a remarkable instrument in modelling damage in compressible structures subject to extensive pulling. However, it has exhibited both volumetric and shear locking in incompressible specimens under compression and bending status, respectively. Therefore, in order to overcome these locking pathologies, the 24-incompatible modes EAS technique Q1Q1E24 and the mixed displacement-Jacobian-pressure FE formulation Q1Q1P0 have been proposed for their use. Q1Q1E24 scheme is proven to solve the shear locking phenomenon in compressible samples, while the employment of Q1Q1P0 formulation, by considering both the pressure and the dilatation as separate DOFs, has been demonstrated to be a more than a remarkable instrument in modelling complex incompressible cases, correcting the volumetric locking pathology.
We also have demonstrated that Q1Q1E24 element performs poorly in predicting damage in incompressible materials, while Q1Q1P0 overestimates the curve in structures subject to bending loads, i.e., displaying shear locking. Therefore, to tackle this, further approaches combining EAS with the mixed Jacobian-pressure formulations will be proposed in the future. In addition to this, an extension to model anisotropic hyperelastic almost incompressible materials is proposed here for future work, being envisaged as a compelling but challenging extension of the proposed frameworks, considering the distortion that such kind of materials exhibit during pure pressure loading [108].
Funding Open access funding provided by Scuola IMT Alti Studi Lucca within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.