Stabilized mixed formulation for phase-field computation of deviatoric fracture in elastic and poroelastic materials

In the numerical approximation of phase-field models of fracture in porous media with the finite element method, the problem of numerical locking may occur. The causes can be traced both to the hydraulic and to the mechanical properties of the material. In this work we present a mixed finite element formulation for phase-field modeling of brittle fracture in elastic solids based on a volumetric-deviatoric energy split and its extension to water saturated porous media. For the latter, two alternative mixed formulations are proposed. To be able to use finite elements with linear interpolation for all the field variables, which violates the Ladyzenskaja–Babuska–Brezzi condition, a stabilization technique based on polynomial pressure projections, proposed and tested by previous authors in fluid mechanics and poromechanics, is introduced. We develop an extension of this stabilization to phase-field mixed models of brittle fracture in porous media. Several numerical examples are illustrated, to show the occurrence of different locking phenomena and to compare the solutions obtained with different unstable, stable and stabilized low order finite elements.


Introduction
Porous media are materials characterized by a heterogeneous internal structure consisting of a solid phase, which confers stiffness to the material, and empty spaces, called pores, which may be filled by one or more fluids [26]. Due to this complicated internal microstructure, for engineering purposes it is more convenient to model the behavior of these materials at the so-called macroscopic scale, i.e. substituting the real structure of the material by ideally superimposed continua which occupy the entire domain at the same time [14]. At this scale the classical balance equations characteristic of continuum mechanics can be applied without taking into account the discontinuity at the interfaces between the real constituents. In order to derive the macro-B Claudio Gavagnin c.gavagnin@tu-braunschweig.de 1 Dipartimento di Ingegneria Civile, Edile e Ambientale, Università degli Studi di Padova, Via F. Marzolo 9,35131 Padua, Italy 2 Institute of Applied Mechanics, Technische Universität Braunschweig, Pockelsstr. 3, 38106 Braunschweig, Germany scopic equations describing the behavior of these substitute continua, several theories have been developed in the field of the mechanics of porous media. These theories can be classified into two main categories [26]: one approach starts from the mechanical description of the behavior of the real constituents at the microscopic scale, with the derivation of balance equations in which interaction forces and discontinuities at the interfaces between the constituents are taken into account. The macroscopic equations are then derived using an averaging process based on the integration of the microscopic quantities over a control volume, called representative volume element. These averaging theories, based on the upscaling of microscopic quantities, are known as hybrid mixture theories [20][21][22]37]. A second possible approach is to start the derivation of the balance equations for the single phases directly at the macroscopic level, and is based on the fundamental concept of volume fraction, defined as the ratio of the volume of the constituents to the volume of the control space [26]. The volume fraction allows the smearing operation of the intrinsic mechanical properties of the single phases over the entire domain . Phenomenological approaches, such as Biot's theory [3,4] and mixture theories [34][35][36] restricted by the volume fraction concept belong to this category [11,18,29,30]. A comprehensive description of the balance equations that govern the (thermo)-hydromechanical behavior of porous media can be found in [26].
Concerning the modeling of fracture, a fundamental contribution to the study of the problem of evolution of preexisting cracks in elastic materials is the energetic approach proposed by Griffith [19]. Griffith's criterion states that a preexisting fracture in an elastic body can propagate only if the release of potential energy upon an infinitesimal increment of the fracture surface equals the surface energy associated to the increment itself. A variational formulation of Griffith's energetic criterion was first proposed by Francfort and Marigo [16]. The phase-field approach to fracture [9] is based on the regularization of this variational formulation. A review on several phase-field formulations of brittle fracture present in the literature can be found in [1].
In order to couple the phase-field approach to fracture with porous media mechanics, several approaches can be followed. A variational formulation of the coupled problem, limited to the saturated case, has been proposed in [27]. Other possible strategies to realize the coupling can be found in [8,23,28]. In this work we follow the approach proposed in [13], where the coupling between the two problems is realized including a dependency on the phase-field parameter into the constitutive law for the effective stress, which is defined as the portion of the stress directly correlated with the elastic deformations of the solid matrix of the porous medium [32].
In the numerical approximation of phase-field models of fracture in porous media with the finite element (FE) method, the problem of volumetric locking may occur. The causes can be traced both to the hydraulic and the mechanical properties of the material. Regarding the hydraulic properties, it is well known in the literature [7,26,39] that, in the early stage of the consolidation process of porous media with low permeability, oscillating solutions for the pore pressure can be found when using FEs that violate the so-called Ladyzenskaja-Babuska-Brezzi (LBB) condition [12]. Regarding the mechanical properties, it has been shown in [33] that, when using a deviatoric-volumetric energy split in the phase-field simulation of brittle fractures subjected to compressive loading, the phase-field solution may feature a much larger thickness of the localization band than that expected based on the chosen characteristic length, or the simulation may even not converge at all. To the best of the authors' knowledge, the combination of these two phenomena has not been studied yet.
This paper focuses on the study and the development of FE formulations which mitigate the aforementioned locking phenomena. For this purpose, we introduce a mixed u− p−d (displacements/pressure/phase-field) formulation for phasefield modeling of fracture in elastic solids and its extension to water saturated porous media. The latter results in two alternative mixed formulations of the problem. In the FE discretization, the higher degree of interpolation for the displacement field due to the LBB condition, in combination with a fine discretization required in order to resolve the length scale inherent to the phase-field approach, leads to a very large number of degrees of freedom (DOFs). Therefore, to allow the use of linear interpolation for all the fields, a stabilization technique based on polynomial pressure projections is applied. This stabilization has been proposed in the literature for the solution of the incompressible Stokes equations [6,15], and successfully applied to the solution of the Darcy problem [5], the coupled problem of saturated fluid flow in a deforming porous medium [38] and for the stabilization of the fluid pressure fields in phase-field modeling of fracture in partially saturated porous media [24]. This technique has two interesting properties: the stabilization term is local (thus requiring only a few additional lines in the FE code) and its influence is controlled by a material parameter [38]. In this work we propose an extension of this stabilization to phase-field mixed models of brittle fracture based on a volumetric-deviatoric energy split. This extension stems from the analogy between the equations governing the problem of consolidation in water saturated porous media under undrained conditions and the mixed u − p formulation for incompressible elasticity.
The paper is organized as follows: Sect. 2 focuses on modeling fully saturated porous media and introduces the mixed formulation and the stabilization technique for this case. In Sect. 3 the same stabilization technique is extended to phasefield modeling of brittle fracture in elastic solids and, in Sect. 4, the two problems are combined. Here two alternative mixed formulations, based on the mean effective stress and on the mean total stress, respectively, are proposed. Several numerical examples are analyzed to compare the solutions obtained with unstable, stable and low order stabilized FEs. Conclusions are outlined in Sect. 5.

Governing equations
The governing equations for the coupled problem of water flow in a deforming porous medium are briefly derived.
The following assumptions are made: isothermal and fully saturated conditions (u − p w , displacement/water pressure formulation), incompressibility of both constituents, and geometric linearity. With these assumptions, we can write the strong form of the equilibrium equation as and the strong form of the mass balance equation of solid and water as where σ t is the Cauchy total stress tensor, g the vector of gravity acceleration, u the solid matrix displacement, v ws the relative water velocity with respect to the solid matrix, and is the density of the mixture, with ρ s and ρ w indicating the intrinsic density of the solid and of the water phase, respectively, and n indicating the porosity. Equations (1) and (2) are defined on the domain of the body. A superimposed dot denotes derivation with respect to time, while ∇· is the divergence operator. Based on Terzaghi's effective stress principle [32], the Cauchy total stress tensor σ t can be decomposed as where p w is the water pore pressure, σ e is the so-called effective stress, namely the portion of the total stress directly related to the deformation of the solid matrix, and I is the second order identity tensor. To solve the system of differential Eqs. (1) and (2) some constitutive laws and the strain-displacement relation need to be introduced. Concerning the mechanical behavior of the solid matrix, we assume a linear elastic constitutive relation between the effective stress tensor σ e and the infinitesimal strain tensor ε, namely σ e = C : ε (5) where C is the fourth-order elastic tensor. Note that ε = ∇ s u, where ∇ s is the symmetric part of the gradient operator.
For the water flow we introduce, as constitutive equation, Darcy's law, which, for a saturated porous medium, takes the form where k w is the isotropic hydraulic conductivity (expressed in m/s). Introducing (5), (6) and (4) into (1) and (2) we obtain the system of governing equations which can be solved for the primary variables u and p w . In order for the system of Eqs. (7) and (8) to be a wellposed problem, we need to specify some boundary conditions (BCs) on the boundary = ∂ , together with some initial conditions (ICs) at the time t = 0 on ∪ . The ICs and the BCs are expressed as where u 0 and p w 0 are the initial values of the displacement and water pressure, u andp w are the imposed values of the displacement and water pressure on the Dirichlet boundary D u ∩ D p , t andq are the imposed values of the traction and the water flux on the Neumann boundary N u ∩ N p , and n is the outward unit normal to .

Weak form and FE discretization
Defining the following spaces for the trial functions u and p w and the test functions w u and w p (10) the solution of the problem (7), (8) and (9) is the pair {u, p w } ∈ T u × T p that solves, for any admissible pair w u , w p ∈ W u × W p , the following weak form of the problem Starting from this weak form, it is now possible to discretize the problem in time and in space. Concerning the discretization in time, we denote with n + 1 and n the quantities at the current and the previous time steps, respectively, with t = t n+1 − t n the size of the current time step and we apply the Backward Euler scheme. We obtain the following discrete-time counterpart of (11): Concerning the discretization in space, we subdivide the domain in a mesh of FEs, and we consider an approximation of the trial and test functions based on polynomial shape functions with local support, namelỹ (13) where() are the approximated trial and weighting functions, () the vectors containing the values of those functions on the mesh nodes, N p = N 1 p N 2 p ... N n n p is the vector containing the shape functions for the water pressure at the n n nodes of the mesh and N u a matrix defined as where (15) with N i u as the shape function for the displacement relative to the node i. Furthermore we introduce the matrices B and b, defined such that Bû = {ε} and bû = ∇ ·ũ, with {ε} being the Voigt (vector) representation of the infinitesimal strain tensor ε. The same shape functions are used to approximate the corresponding trial and weighting functions, recovering the standard Bubnov-Galerkin Method.
We can therefore derive the discrete counterpart of the weak form (12) as The matrix C is the elasticity matrix, which for the plain strain assumption adopted in the numerical simulations presented in this paper, reads where λ and μ are the Lamè coefficients. The quantities (17)- (21) do not depend on the solution at the current time step, so the problem is linear and the solution can be directly calculated solving the system The next step is the choice of the functions N u and N p . It is well known in the literature that the pair (N u , N p ) has to fulfill the so called LBB condition to ensure the stability of the solution in the locally undrained limit (e.g. [39]). The study of the stability of the FE formulation will be treated in the next section.

Stable, unstable and stabilized FE formulations
We consider now the algebraic structure of the problem in the local undrained limit, i.e. when k w t → 0 [38]. In this case the system of Eq. (23) becomes We can notice that the lower diagonal block in the system (24) is now zero, and that the matrix assumes the typical structure of incompressible elasticity problems in solid mechanics. From a physical point of view, this is exactly what is happening: if we consider an arbitrary control volume in our domain, in the undrained limit the water can not flow in or out and, due to its very low compressibility, constrains the solid matrix to develop pure deviatoric deformations. In order to obtain a stable solution, the choice of the spaces W u and W p has to satisfy the LBB condition. The combination of quadratic displacements and linear water pressure interpolations (Taylor-Hood element Q9P w 4) satisfies this condition, while the low order element with linear displacement and linear pressure interpolations (Q4P w 4) does not [26,38,39]. Obviously, the element Q9P w 4, with its 22 DOFs, is computationally more expensive than the Q4P w 4, which has only 12 DOFs. The computational advantage of the linear-displacements/linear-pressure interpolation becomes even more evident in three-dimensional problems, where every hexahedral element has 32 DOFs, while using a quadratic-displacements/linear-pressure interpolation the number of DOFs per hexahedral element rises to 89.
A smart approach to stabilize the linear-displacements/ linear-pressure interpolation has been introduced by Bochev and Dohrmann in [5,6,15] in the context of fluid mechanics, and successively applied by White and Borja in [38] in porous media mechanics. The idea proposed in [15] is to add a stabilization term in the weak formulation of the problem, resulting in a non-zero lower diagonal block in the corresponding discrete equation. Let us first define the element piecewise constant projection operator as where Π is the projection operator applied to the generic variable and V e is the volume of the element domain e . When applied to the water pressure p w , this is nothing but the definition of the average element water pore pressure. In [38] White and Borja proposed the following modified version of (16): where the stabilization vector R stab n+1 has been added. This is defined in [38] as with S pw being the so-called stabilization matrix, defined as where c W B stab is a penalization factor, which regulates the influence of the stabilization matrix S pw .
The idea behind this stabilization is substantially to 'penalize' the oscillations of the pressure around their element average value. This is actually more than a simple penalization: Bochev and Dohrmann demonstrated that this additional term quantifies, and therefore corrects, the deficiency of the linear-displacements/linear-pressure interpolation with respect to the LBB condition [15,38]. There are two main advantages in the use of this stabilization: • The stabilization matrix S pw is a local additional term: it requires only the computation of the average value of the element shape functions, making this stabilization easy to implement and adding a minimal computational cost. • The stabilization coefficient c W B stab acts as a penalization factor in (28) and can be related to a physical parameter characterizing the problem. In [15] this factor is taken as the inverse of the dynamic water viscosity, while in [38] it is chosen as where τ is an additional coefficient introduced to compensate a possible excess of stabilization (τ = 0.04 has been used in [38] for a two-dimensional numerical application).
Being the permeability k w > 0 (even if very small), the problem is always characterized by a transient evolution (even if very slow) and, therefore, the stabilized quantity is the time derivative of the water pressure. Now it is interesting to notice that, in analogy to the idea in [15] for the stabilization of the Stokes equations, White and Borja [38] related the penalization factor to the inverse of the material parameter 2μ. However, while in the Stokes equations the parameter multiplying the velocity is a scalar value and so the choice of the penalization is quite clear (the inverse of that value), in the present case the strain tensor ε is multiplied by the fourthorder elasticity tensor C, containing the two Lamè constants λ and μ. Thus a criterion for the choice of a scalar parameter representative of the tensor C is needed. In analogy with [15] we suggest to use as a stabilization coefficient with n d being the dimension of the problem, where the value C 1111 = C 2222 = ... = C n d n d n d n d has been chosen as scalar parameter representative of the tensor C. We notice that, with our choice, the stabilization coefficient is always smaller than or equal to the coefficient c W B stab , proposed in [38], leading to a weaker stabilization effect. On the other hand, we do not introduce any additional coefficient τ as in (29).
With the introduction of the stabilization matrix S pw , the stabilized counterpart of (23) can be derived, obtaining

Numerical example
In this section we simulate the consolidation process of a vertical column of water saturated soil [38], known as Terzaghi's problem, whose analytical solution is known [31]. The domain of the problem is a layer of poroelastic material of base L and height H , subjected to a constant surface load w applied at the top and vertically constrained at the bottom (see Fig. 1). The water flux is free at the top and set to zero at the bottom. Concerning ICs, the water pressure is set to zero on the entire domain. On the lateral boundaries zero horizontal displacement and water flux are imposed. Plain strain is assumed. The mesh used for the discretization of the domain consists in a column of 20 quadrilateral elements. Our goal is to compare the numerical solutions of the water pressure obtained using different stable and unstable FEs at  [38] the early stage of the consolidation. The problem is hence solved using the Taylor-Hood element Q9P w 4, the low order element Q4P w 4, and the stabilized elements Q4P w 4s-WB (the one proposed by White and Borja with the stabilization coefficient (29)) and Q4P w 4s (with our modified stabilization coefficient (30)). The material parameters are: E = 1 kPa, k w = 10 −5 m/s, ρ s = 0, ρ w = 1000 kg/m 3 , g = 0 m/s 2 [38]. The time step is Δt = 1 s. Two different Poisson's ratios are considered, ν = 0 (λ = 0 kPa, μ = 0.5 kPa) and ν = 0.4 (λ = 1.43 kPa, μ = 0.36 kPa). The latter is a very high value for soils, but is chosen here to amplify the difference between the choices (29) and (30) of the penalization coefficient. Immediately after the application of the load, for a factor k w Δt small enough, the porous medium is subjected to local undrained conditions, behaving as an incompressible material: the water bears the entire load and the analytical solution consists in null vertical displacements and constant water pressure equal to the applied load along the column, with a sharp jump in the solution at the top, where Dirichlet BCs are applied. Under these conditions unstable FEs are expected to show oscillations for the pressure solution. Figure 2a shows that, for the case ν = 0, the solution obtained for the water pressure with the element Q4P w 4 features large oscillations that propagate also far away from the boundary where the load is applied (pressure jump), while the solution obtained with the stable Taylor-Hood element Q9P w 4 shows a small oscillation, that disappears at a small distance from the boundary. Moreover, Fig. 2a shows also that both stabilized elements Q4P w 4s-WB and Q4P w 4s reproduce perfectly the solution of the stable element Q9P w 4. The two solutions coincide as, for ν = 0 (i.e. λ = 0), (29) and (30) deliver the same result.
In Fig. 2b the results obtained using ν = 0.4 are shown. For the element Q4P w 4 the results show again large oscillations along the entire domain, while the Taylor-Hood element, as expected, leads to a stable pressure distribution. In this case, however, the difference in the choice of the stabilization coefficient becomes evident. While the proposed stabilized element Q4P w 4s delivers again exactly the same solution obtained with the stable Q9P w 4, Q4P w 4s-WB leads to no oscillation. This result, which seems to indicate an advantage of Q4P w 4s-WB, is instead an indicator of an excess of stabilization. In fact, we consider as indicator of the correctness of the stabilization coefficient the agreement with the numerical solution obtained with a stable formulation, rather than the agreement with the analytical solution. This excess of stabilization for the element Q4P w 4s-WB was also noticed by the authors themselves, and could be the reason of the additional coefficient τ < 1, introduced in (29) in other to reduce the excess of diffusion obtained in some simulations [38]. Except for the choice of the stabilization coefficient, the stabilization technique proposed in [38] is very effective and, at the same time, easy to implement.

Energy functional
As already metioned in the introduction, the phase-field approach to brittle fracture is based on the work of Bourdin et al. [9], and consists in the regularization of the variational formulation of Griffith's theory of brittle fracture, first proposed in 1998 by Francfort and Marigo [16]. The problem of equilibrium and quasi-static evolution of the phase-field variable d, assumed to vary smoothly between the values d = 0 (intact material) to the value d = 1 (fully broken material), is governed by the minimization of the energy functional where g (d) = (1−d) 2 +η is the so-called degradation function, with 0 < η 1 being a residual stiffness, introduced for numerical stability reasons, Ψ + (ε) and Ψ − (ε) are the so-called 'positive' and 'negative' parts of the undamaged elastic energy Ψ (ε) = 1 2 λtr 2 (ε) + μtr(ε 2 ), whose definition depends on the particular model chosen for the energy split, G c is the fracture toughness, C v is a normalization constant, depending on the model used for the local part of the dissipation function w(d), and l is a characteristic length parameter controlling the amount of regularization. The last integral is introduced to enforce the irreversibility of the evolution of the phase-field variable d [17].
The functional (32) is a general expression, which needs to be particularized by choosing a specific model for the function w(d) and for the the split of the energy Ψ = Ψ + + Ψ − . Concerning the local part of the dissipation function, we choose the so-called AT1 model [10], characterized by the following expression for w(d) and C v : Regarding the choice of the split of the energy density function, we introduce the volumetric-deviatoric energy split proposed by Lancioni and Royer-Carfagni [25]. In this model it is where ε dev = ε − 1 3 tr(ε)I is the deviatoric component of the strain tensor ε, K = λ + 2 3 μ is the bulk modulus of the material, C dev and C vol are the deviatoric and the volumetric parts of the elasticity tensor C, respectively. Inserting (33) and (34) in (32), we obtain We notice at this point that, to model the real behavior of brittle materials, a more adequate choice for a split based on the volumetric-deviatoric decomposition of the elastic energy would be the one proposed by Amor et al. [2] where implying that also the volumetric energy contributes to the fracture evolution, if the volumetric deformation is positive. Anyway, the numerical instabilities due to volumetric locking, which are in focus in this paper, are expected to manifest themselves when the volumetric deformation is negative, in which case the splits (34) and (36) coincide.

Minimization problem
The phase-field formulation of brittle fracture consists in the minimization of E l (u, d) with (u, d) ∈ T u × H 1 ( ). Necessary conditions for this minimization are are the directional derivatives of the functional E(u, d) with respect to u and d in directions v and α, respectively. Applying Green's Lemma to (37) and (38), it is possible to derive the following Euler's equations where σ (ε, d) is the stress tensor, defined as Together with the BCs Equations (39) and (40), i.e. the equilibrium equations and the phase-field evolution equation, respectively, represent the strong form of the variational problem.

Mixed formulation and stabilization
Because of the volumetric-deviatoric energy split introduced in the functional (35), when a fracture occurs the deviatoric stiffness of the material becomes very small in comparison with the volumetric one, which is not affected by the phase-field variable d. Therefore, a discretization that cannot treat the deviatoric and the volumetric part of the deformation independently can lead to instabilities when the fracture localizes. This problem is similar to the well known phenomenon of volumetric locking in the modeling of incompressible materials, and can be solved using a mixed formulation in which, in addition to the displacement u, also the hydrostatic component of the stress p, defined as is treated as a primary variable. The latter is related to the volumetric deformation of the body by the constitutive equation Following the procedure in [33], we now introduce a mixed formulation u − p − d of the phase-field model presented in Sect. 3.1. Taking into account the elastic constitutive relation (45), the volumetric part of the elastic energy can be expressed as a function of p, namely With this definition the functional (35) can be rewritten as The evolution of the pressure p is related to the volumetric deformation by the constitutive equation (45), which can be viewed as a constraint equation, namely This constraint equation can be included in the variational formulation using the method of Lagrange multipliers. Indeed, if we define the following Lagrangian where β is a Lagrange multiplier, introduced to enforce the constraint (48), the variational problem for the resulting mixed formulation consists in the stationarity of L(u, p, d, β) with (u, p, d, β) ∈ T u × L 2 ( ) × H 1 ( ) × L 2 ( ). Necessary conditions for stationarity are From the condition we obtain that This means that the hydrostatic stress p acts as a Lagrange multiplier in the mixed formulation. The Lagrangian (49) becomes therefore and its directional derivatives read for every (v, q, α) ∈ W u × L 2 ( )× H 1 ( ). After the application of Green's lemma to Eqs. (54) and (56), we can derive the strong form of the mixed variational formulation, which consists in the differential equations valid on the domain , together with the BCs

FE discretization
We introduce now an approximation of the trial and test functions based on polynomial shape functions with local support, namelỹ is the vector containing the shape functions for the phase-field at the n n nodes of the mesh. The following discrete version of the weak form of the problem is obtained, where R u , R p and R d are the residual vectors and C dev and C vol are the deviatoric and the volumetric com-ponent of the elasticity matrix C which, in the case of plain strain, read The system of Eqs. (64), (65) and (66) is solved using a staggered approach in which, for every loading step, the discrete coupled u− p problem, consisting in the system of Eqs. (64) and (65), and the discrete phase-field evolution equation (66) are solved independently, using the last calculated value of the other fields as a constant. These procedure is repeated until a convergence condition is satisfied, namely [33] where s + 1 and s are the current and the previous staggered iteration, respectively. Due to the presence of the Macaulay brackets, the residual R d is nonlinear with respect to d, so an iterative procedure has to be used to solve (66) at each staggered iteration within the current loading step. Using the Newton-Raphson method, the solution is obtained solving, at each nonlinear iteration, the system of equations Concerning the coupled u − p problem, Eqs. (64) and (65) are solved with a monolithic approach, which leads to the following linear system of equations: and Q and f u are as defined in Sect. 2.2.
We focus now on the stability of the formulation which, as we saw in the previous section, is ensured by a choice of the shape functions N u and N p satisfying the LBB condition. The combination of quadratic displacement and linear pressure interpolations (Q9P4 for quadrilateral FEs) satisfies this condition, while the combination of linear displacement and linear pressure interpolations (Q4P4 for quadrilateral FEs) does not. If we look at the structure of the system of Eqs. (72), when K → ∞, we can notice a clear analogy with the structure of the system (24), characterizing the numerical solution of the problem of consolidation of a porous media under undrained conditions. This is either the case of materials with very high volumetric stiffness K , or the case of cracked materials modeled with a phase-field approach based on a volumetric-deviatoric split of the energy, characterized by a deviatoric stiffness several order of magnitudes smaller than the volumetric one in correspondence of values d 1 of the phase-field.
In order to stabilize the solution obtained with the elements Q4P4, we use the pressure projection technique introduced in Sect. 2.3. Therefore we add to the system (72) a stabilization matrix S p f defined as where, as shown in Sect. 2, the choice of the stabilization coefficient c p f stab plays an important role in the performance of the stabilized solution. Following (30), the stabilization coefficient, which now depends on a scalar value representative of the damaged deviatoric part of the elasticity tensor C, is defined as With the introduction of the stabilization matrix S p f , the following stabilized counterpart of the mixed formulation (72) is obtained

Numerical examples
We present now two numerical examples in order to test the performance of different FE interpolations. In particular we consider interpolations with the quadrilateral elements Q4/Q4 and Q9/Q4 for the formulation u − d, and with the quadrilateral elements Q4P4/Q4, Q4P4s/Q4 and Q9P4/Q4 for the mixed formulation u − p − d. These symbols are specified in Table 1, including an indication of the field in which, if present, the stabilization is applied. In the notation used to identify the elements, the interpolation of the phasefield d is separated by the symbol "/", to point out that we use a staggered procedure in which the phase-field equation is solved separately from the other equations. We choose a linear interpolation for the phase-field for all the FEs considered, for a better comparison of the results. With regard to the choice of the penalty coefficient introduced to enforce the irreversibility of the phase-field, we use the expression derived in [17], which for the AT1 model and a tolerance TOL ir = 0.01 on the error in the evalution of the fracture energy, leads to the following expression for the penalty parameter: This expression is used in all the phase-field numerical simulations in this paper.

Uniaxial tension test
In this example we study the numerical solution of the uniaxial tension test of a bar, with the dimensions shown in Fig. 3. The mesh is composed by a grid of square elements, with side h el = 0.005 m, and plane strain conditions are assumed. A Dirichlet BC u x = −t, with t being a pseudo-time parameter, is imposed on the left boundary, while the right one is assumed fixed. Concerning the phase-field variable d, homogeneous Dirichlet BCs are imposed on the left and the right boundaries, in order to maintain the symmetry of the phase-field solution, thus promoting the localization of the phase-field in the middle of the domain. The parameters of the model are: E = 1 Pa, ν = 0.4, G c = 1 N/m, l = 0.05 m, η = 10 −8 . The time step is t = 0.05 s. We perform 80 time steps, with a tolerance TOL p f = 10 −8 on the residual for the phase-field Newton-Raphson iterative scheme, and a tolerance TOL stag = 10 −5 for the staggered loop (see Sect. 3.4). The localization of the phase-field variable occurs at t loc = 2.7 s. Figure 4a plots the profiles of the phase-field at t end = 4 s, computed with the standard elements Q9/Q4 and Q4/Q4, and the analytical solution derived in [33]. Looking at the profile obtained with the element Q9/Q4, we can notice that the analytical solution is almost perfectly recovered, with the exception of the small plateau obtained in correspondence of the cusp of the analytical optimal profile, which is due to the numerical approximation and, having the dimension of one FE, tends to disappear for h el → 0. On the other hand, the element Q4/Q4 shows a localization width significantly larger than the analytical one, with a round profile at the peak instead of the cusp that characterizes the analytical solution. This is how the locking phenomenon occurs when the phasefield localizes, caused by the fact that the deviatoric stiffness, degraded by the phase-field, becomes very small in comparison with the volumetric one.
In Fig. 4b we compare the analytical profile with the one obtained with the mixed elements Q4P4/Q4, Q4P4s/Q4 and Q9P4/Q4. In this case we can notice that all the solutions reproduce very well the analytical one. However the profile obtained with the element Q4P4/Q4 is slightly wider than the one obtained with the stabilized element Q4P4s/Q4, whose solution coincides with that of the stable Q9P4/Q4.
In order to compare the general behavior of all the standard and mixed FEs adopted so far, in Fig. 5 we plot the forcedisplacement curve relative to all the numerical solutions. We notice that, after the localization of the phase-field, the reaction force vanishes for all the elements, except for the Q4/Q4. This is due to an excess of residual stiffness, again characteristic of a locking behavior. A confirmation of this can be obtained comparing the maximum values of d at the peak of the damage: while for all the other elements d reaches  According to these results, except for a small difference in the width of the localization profile, the mixed element Q4P4/Q4 seems to behave as well as the elements Q9/Q4, Q4P4s/Q4 and Q9P4/P4. However, the positive effect of the stabilization becomes evident if we look at the profile of the hydrostatic stress p, obtained with the mixed elements, after the localization. In Fig. 6 we can clearly observe that, in correspondence of the crack, the profile obtained with the stabilized element Q4P4s/Q4 is in good agreement with the one obtained with the stable Q9P4/Q4, while the unstable element Q4P4/Q4 leads to large oscillations.

Two-dimensional compression test
In this example we study the development of a shear fracture in a rectangular domain subjected to lateral compression. This problem, with different dimensions, has already been solved in [33], in order to show the stability of the mixed element Q9P4/Q4, compared to the standard one Q4/Q4. We recompute this example with the elements used in the previous section. The mesh is composed by square elements with side h el = 0.005 m, and the dimensions of the domain are shown in Fig. 7 along with the BCs. In this case, on the left side we apply symmetric BCs, so the vertical displacements u y and the phase-field d are not constrained. The problem is solved in control of displacements, with u x = −t as non-homogeneous Dirichlet BC on the right boundary. The parameters of the model are: E = 1 k Pa, ν = 0.3, G c = 1 N/m, l = 0.02 m, η = 10 −8 . The time step is t = 2.5e −4 s. We perform 320 time steps, with a tolerances TOL p f = 10 −8 and TOL stag = 10 −5 (see previ-  Figure 8 shows the contours of the phase-field d obtained at the end of the simulation with the different FEs. Looking at the solution obtained with standard elements, we can again notice that the element Q9/Q4 leads to the expected solution, i.e. a fracture that localizes along an inclined shear band starting from the center of the domain (considering the symmetry). The standard element Q4/Q4, instead, localizes along a thick vertical crack, whose width is no longer controlled by the internal length l. Focusing on the mixed elements, we notice that the stable Q9P4/Q4 and the stabilized Q4P4s/Q4 are able to catch the correct crack pattern and width, while the unstable Q4P4/Q4 shows, as in the previous example, a width slightly broader than the one obtained with the stabilized Q4P4s/Q4, but with a non realistic localization pattern, which starts vertically from the center and then tilts down along an inclined band.

Governing equations
In this section we aim at coupling the model for fluid flow in saturated elastic porous media, formulated in Sect. 2, with the phase-field model for fracture with volumetric-deviatoric energy split, given in Sect. 3. If we combine the phase-field equation (40) with the system of Eqs. (7) and (8), we obtain where the linear elastic constitutive equation (5) for σ e has been replaced by the following constitutive relation which is the same relation used in Sect. 3 for σ . The ICs and BCs of the problem are u = u 0 at t = 0 Following the same procedure of Sect. 2, after the discretization in time and space, we obtain the following expression for the residuals: Again, as explained in Sect. 3.1, the problem is split into two sub-problems, solved for the pair (u, p w ) and for the phase-field d independently. The residual R d is identical to the one obtained in Sect. 3.4 and the same nonlinear solution scheme is adopted. The sub-problem (84) (85) is solved with a monolithic approach, which leads to the linear system In this case, when deviatoric fractures develop under locally undrained conditions, both the instabilities analyzed in Sects. 2 and 3 can appear. The use of a stabilized element Q4P w 4s for the problem (87), together with a linear interpolation (Q4) of the phase-field, is expected to be stable under undrained conditions, and to become unstable when the fracture develops. The Taylor-Hood element Q9P w 4/Q4 is stable with respect to the undrained conditions and, as shown in the numerical example in Sect. 3.5, can lead to acceptable results for the phase-field. However, none of these elements ensures a stable behavior when the phase-field localizes. In the next section we derive two possible mixed formulations, which are stable under both undrained and fractured conditions.

Mixed formulations and stabilization
We focus now only on the poromechanical sub-problem governed by Eqs. (79) and (80), based on the fact that the phase-field equation is solved independently in the staggered procedure adopted. Following the same procedure of Sect. 3.3, we can derive a mixed formulation of equation (79) leaving unchanged equation (80). We obtain the system of differential equations where p e is the hydrostatic component of the effective stress tensor σ e , defined through the effective stress principle (4). Equation (89) is the volumetric part of the elastic constitutive equation (82). Following the usual FE procedure, we obtain the following system of algebraic equations: to be solved along with the discrete counterpart of the phasefield equation (86). Based on the main variables of the problem, we call this the "u − p e − p w − d formulation". Now, looking at the structure of the matrix on the lefthand side of the system of Eq. (91), we can notice that both equations for the fields p e and p w are coupled with the equations for the displacement field u. In the perspective of the stability of the mixed formulation it will be necessary either to use stable FEs, like the quadrilateral element Q9P e 4P w 4 (quadratic-displacement/linear-effective pressure/linear-water pressure), or to develop a low order FE defining a stabilization for both fields, dealing though with the problem of the correct balance between the two stabilizations.
A possible way to reduce the number of couplings is to perform a change of variables, based on the effective stress principle restricted to the volumetric component of the stress, namely where p t is the hydrostatic part of the total stress tensor σ t . After applying the change of variables (92), the system of Eqs. (88), (89) and (90) becomes We can notice that equation (95) is still directly depending from u. But if we take the total time derivative of equation (94) and we combine it with (95) we obtain the system of differential equations This leads to the discrete system of equations: to be solved together with the phase-field equation (86). Based on the main variables of the problem, we call this the "u − p t − p w − d formulation". Now we can notice that p w is not directly coupled with u anymore.
Both the previous formulations are expected to be stable if we use for u a order of interpolation higher than the one used for the pressures. This is the case of the element Q9P e 4P w 4/Q4 for the u − p e − p w − d formulation, and of the element Q9P t 4P w 4/Q4 for the u − p t − p w − d formulations (the order of interpolation for the phase-field is always assumed linear).
It is clear that the number of DOFs of a stable FE, already large for the u − p w − d formulation, becomes even larger for the two developed mixed formulations. For example, if we focus on the interpolation for the poromechanical subproblem, the standard Taylor-Hood element Q9P w 4 implies 22 DOFs, while the mixed Q9P e 4P w 4 (or Q9P t 4P w 4) element has 26 DOFs.
Using the results obtained in the previous two sections for the stabilization of p w and p independently, we can now derive the stabilized versions of the two mixed formulations developed in this section, in order to be able to use FEs with linear interpolations for all the fields and reduce the number of DOFs. As a general guideline, we assume that the stabilization of a particular pressure field is needed when that field is directly coupled with the displacement u. It follows that two equations need to be stabilized in (91), while in (99) only one.
The stabilized version of the discrete system of Eq. (91) is defined as ⎡ where S pw and S p f are the two stabilization matrices introduced in Sects. 2.3 and 3.4, respectively. With these stabilization terms it is now possible to use the element Q4P e 4sP w 4s, with linear interpolations for all the fields (the letter 's' indicates as usual the stabilization of the corresponding field).
Similarly, the stabilized version of the discrete system of Eq. (99) is defined as Also in this case, due to the introduction of the stabilization matrix S p f , it is now possible to use the low order element Q4P t 4sP w 4/Q4.
To conclude we notice that, using the stabilized element Q4P e 4sP w 4s (or Q4P t 4sP w 4) instead of the stable Q9P e 4P w 4 (or Q9P t 4P w 4), the number of DOFs per element reduces   Table 2.

Terzaghi's consolidation problem
In this first example we solve again the Terzaghi's problem analyzed in Sect. 2.4. Our goal is a first comparison, without the influence of the phase-field d, of the results obtained with the two formulations u − p e − p w and u − p t − p w with those obtained with the u − p w formulation. The FEs used in the simulation are the stabilized Q4P e 4sP w 4s and the stable Q9P e 4P w 4 for the u − p e − p w formulation, the stabilized Q4P t 4sP w 4 and the stable Q9P t 4P w 4 for the u − p t − p w formulation, and the Taylor-Hood Q9P w 4 for the u− p w formulation. The parameters are the same of Sect. 2.4, with ν = 0.4. In Fig. 9 the pressure profile after the first time step for the different FEs is shown. Both stable mixed elements Q9P e 4P w 4 and Q9P t 4P w 4 lead to the same solution obtained with the standard Taylor-Hood Q9P w 4, used as reference solution. Coming to the stabilized low order elements Q4P e 4sP w 4s and Q4P t 4sP w 4, we notice that the second one leads to the same results obtained with the reference element Q9P w 4, while the first one shows slightly larger oscillations in the vicinity of the top surface. Moreover, the Q4P t 4sP w 4 turns out to be stable for the water pressure field p w , even if the stabilization is applied only to the equation relative to (a) (b) Fig. 9 Pressure distribution along the column at the first time step for different FEs. a u − p e − p w formulation, b u − p t − p w formulation the total hydrostatic stress field p t , confirming our assumption that the stabilization is needed only for the field directly coupled with the displacement field u. On the other hand, the larger oscillations obtained with the Q4P e 4sP w 4s are probably due to the difficulty to define a correct value for the two stabilization coefficients needed in S pw and S p f , which, due to the coupling of the equations, need to be mutually balanced.

Two-dimensional water saturated compression test
In this second example we recompute the compression test of Sect. 3.5.2 considering, in this case, the material as a water  Table 2. The parameters of the problem are the same of Sect. 3.5.2, with the exception of the number of steps, which is equal to 240. The permeability of the porous medium is k w = 10 −8 m/s (a value quite low, corresponding to a silt). The localization of the phase-field occurs at t loc = 0.05525 s (the presence of the water accelerates the process of formation of the crack).
In Fig. 10 we show the results relative to the u− p w −d formulation. We can see how the element Q4P w 4/Q4 (Fig. 11a) shows both oscillations in the water pressure and a wrong localization profile in the phase-field, as already noticed separately in the previous sections. The element Q4P w 4s/Q4 (Fig.  11b) solves the instabilities in the water pressure field, but the phase-field profile still indicates the presence of locking after the phase-field localization. The Taylor-Hood element Q9P w 4/Q4 (Fig. 11c) behaves well also in the coupled simulation, leading to a good solution for both the water pressure and the phase-field variable.
In Fig. 11 we show the results relative to the mixed u − p e − p w − d and u − p t − p w − d formulations. When a quadratic interpolation is used for the displacement field (elements Q9P e 4P w 4/Q4 and Q9P t 4P w 4/Q4), both formulations are stable and lead to identical results (in Fig. 10c only the results relative to the element Q9P e 4P w 4/Q4 are reported). Looking at the phase-field profiles obtained with the stabilized elements with linear interpolation for the displacement field, we notice that the element Q4P e 4sP w 4s/Q4 (Fig. 10a) shows a crack which starts vertically, and then takes the expected direction, while the right failure pattern is perfectly captured by the element Q4P t 4sP w 4/Q4 (Fig. 10b), which shows results very similar to those obtained with the stable mixed element Q9P t 4P w 4/Q4. As for the water pressure field, the element Q4P t 4sP w 4/Q4 leads again to the same solution obtained with the element Q9P t 4P w 4/Q4, while the element Q4P e 4sP w 4s/Q4 shows slightly wider oscillations in the vicinity of the top surface, as noticed already in the previous example (see Fig. 9a).

Conclusions
In this paper we focused on the problem of the numerical locking, due to a condition of high volumetric stiffness of the solid matrix, that can occur using the FE method in the numerical approximation of a phase-field model of deviatoric fracture in elastic solids and in water saturated poroelastic materials. We showed how the causes of this state of incompressibility can be traced both to the hydraulic and to the mechanical properties of the material.
In the first part we reviewed the problem of the numerical locking in the simulation of the consolidation process of water saturated porous media under undrained conditions. The problem appears when a linear interpolation for both the solid and the fluid fields is used, due to the fact that the LBB condition is not satisfied. A stabilization technique presented in [38], developed using a polynomial pressure-projection technique originally introduced for the incompressible Stokes equation in [15], has been adopted and, additionally, a correction of the penalization coefficient introduced in [38] has been proposed, based on an analogy with the coefficient introduced in [15].
In the second part, the problem of volumetric locking due to the presence of fully developed cracks obtained with a phase-field model based on a decomposition of the elastic energy into volumetric and deviatoric parts has been investigated. It has been noticed how, in this case, the locking is due to the fact that the deviatoric stiffness becomes several order of magnitudes smaller than the volumetric one, causing an enlargement of the localization band of the phase-field variable, and a dissipative residual force in the crack. To solve this problem, a stabilized mixed u − p − d FE formulation has been developed, using the same technique applied in the first part of the paper. Finally, the model has been extended to water saturated porous media, and two alternative mixed and stabilized mixed formulations of the problem, based on the Terzaghi's effective stress principle, have been proposed. These two mixed formulations have been tested in two numerical applications, a one-and a two-dimensional consolidation problem, showing how the stabilized mixed formulation based on the hydrostatic part of the total stress tensor p t leads to better results than the one based on the hydrostatic part of the effective stress tensor p e , in terms of crack pattern, phasefield profile and stability of water pressure field. 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/.