This work develops a numerical technique in which the Particle Finite Element method is applied to simulate solid intrusion problems in Geotechnical Engineering. The work describes the numerical developments that made the method functional and showcases its potential with various application problems.
New numerical developments include: integration schemes for largestrain, elastoplastic constitutive models, novel integration methods for rigid body contact constraints and stabilized mixed formulations for singlephase and twophase problems.
An explicit stress integration scheme is developed for elastoplastic, largestrain constitutive models using a multiplicative split of the deformation gradient. This scheme uses adaptive substepping and a yield violation drift correction technique.
Rigid bodies with prespecified shape and motion may be modeled. Contact constraints are introduced using a penalty method in which tangential behavior is treated with an elastoplastic analogy. This elastoplastic contact model is integrated by means of an implicit integration technique; an alternative scheme using the Implex algorithm is also described and assessed.
Loworder (linear) elements are employed to speed computation. These elements may suffer volumetric locking in quasiincompressible conditions. For soils, such conditions appear under undrained loading or at critical state. To alleviate this problem, mixed formulations are developed and stabilization techniques are applied. Two different threefield mixed formulations for the coupled hydromechanical problem are presented, adding either the effective pressure or the Jacobian as nodal variables to the solid skeleton displacement and water pressure. Stabilization terms are used in the mass conservation equation of the biphasic medium and in the rest of scalar equations. Several mixed formulations are also implemented for the simplified singlephase problem, which approximates saturated soil behavior under undrained conditions.
The first application examples are three total stress problems: indentation of a rigid strip footing, TBar motion and a rough cone penetration test. All are frequently used benchmark problems; they allow to compare the performance of the developed numerical scheme with other techniques. It appears that the numerical strategy followed by PFEM, obtains similar results than those attained with alternative numerical methods with significant savings in computational efforts.
A total stress approach is also used in a parametric study of tube sampling in clay. The parameters explored include sampler geometry (roundtipped or beveled cutting shoes; outer diameter to wall thickness ratios); constitutive parameters, roughness factor and boundary conditions are studied. Outputs are analyzed in terms of classical sampling disturbance measures: the centerline strain path and the specific recovery ratio. The results show good agreement with experimental evidence and question the frequently accepted reference role of Strain Path method solutions.
The set of simulations in which a rough interface behavior is considered has been used to assess the theory proposed by [1] to predict the occurrence of a plug inside of an openended pile. The numerical results obtained here corroborate this theory: a plug inside of the tube is formed once the mobilized forces are equal to those that would mobilize a smooth closedended pile. Indeed, the failure mechanism that prevails during the penetration of a plugged tube is coincident with that of a closedended pile.
The last analysis of this work is the hydromechanical simulation of the cone penetration test in a Modified Cam Clay soil. A parametric analysis covers the effect of the permeability of the soil from drained to undrained conditions and the interface friction angle. The effect of these parameters on the cone resistance, sleeve friction and pore pressure at three potential measurement points is fully characterized. These numerical results are used to assess several techniques to estimate the permeability of soils from CPTu testing. Special attention is paid to onthefly techniques, in which permeability could be directly estimated from the CPTu data stream without the need for any stoppage.
This monograph is a revised version of the doctoral dissertation of the first author, entitled Insertion Problems in Geomechanics with the Particle Finite Element Method. The thesis was defended on November 5th 2018 in Barcelona. The research was conducted under the supervision of Professors M. Arroyo, J. M. Carbonell and A. Gens.
The financial support of the Ministry of Education of Spain through research grants BIA201127217 and BIA201459467R is gratefully appreciated. The support provided by the European Commission funded GEORAMP RISE project (H2020645665 GEORAMP) is also acknowledged.
Many activities in Geotechnical Engineering involve the insertion of a rigid body into the soil. Examples include most foundation solutions, that range from displacement piles that may be driven into the soil by striking them with a hammer, spudcan foundations or even suction caissons in which the water pressure inside of the bucket is decreased, causing the foundation to sink into the sea floor.
Penetration problems are also involved in the processes to estimate or obtain constitutive parameters to design these structures. For instance, to perform laboratory tests on a sample, first it has to be retrieved from the ground by pushing or hammering a tube into the ground. Alternatively, constitutive parameters may be obtained from insitu tests, in which instrumented probes are pushed into the soil and the value of constitutive parameters are estimated from test readings.
In this kind of problem large displacements and deformations of the soil mass are always present. The coupled hydromechanical response of the soil adds further complexity, even in cases where insertion speed is tightly controlled. Analysis of problems of rigid body insertion into soil masses had traditionally relied on highly idealized approaches such as geometrically simple cavity expansion mechanisms [2]. Although much insight is gained from such analyses, a number of basic features of the problem are left aside and, as a consequence, a host of not fully understood empirical corrections and methods have been relied upon for practical applications.
Numerical simulation seems an obvious alternative to advance understanding in this area. However, the numerical simulation of rigid body insertion into soils is a complex task, since the system is full of nonlinearities, contactrelated, materialrelated and also geometrical. This latter kind of nonlinearity was a fundamental obstacle to the Lagrangian or updated Lagrangian formulations of the finite element method (FEM) that are successful in other areas of geotechnical engineering. Strong mesh distortion resulted in large inaccuracies and/or stopped calculation at relatively small displacements [3].
In the last decades several numerical frameworks have been developed to address those problems. Some approaches are not based on continuum mechanics and use instead discrete element methods [4,5,6]; however, continuumbased approaches are dominant, particularly for finegrained soils. Within continuumbased methods the approach most frequently applied to geotechnical insertion problems has been that of Arbitrary LagrangianEulerian formulations (ALE). ALE finite element formulations combine the Lagrangian and Eulerian kinematic descriptions, by separately considering material and computational mesh motions [7].
The application of ALE to insertion problems in soil mechanics may be traced to [8]. Afterwards, three main FEM ALE methods have been increasingly developed and repeatedly applied in this area: the socalled remeshing and interpolation technique by small strain (RITSS) developed by Randolph and coworkers [9,10,11], the socalled efficient ALE approach (EALE) developed by Nazem and coworkers [12,13,14] and the successive builtin implementation of ALE in Abaqus/Explicit, currently known as the Coupled EulerianLagrangian (CEL) method [15] which have been applied to insertion problems by several teams [16,17,18]. A comparative review of these ALE methods has been recently presented by [19].
A second continuumbased numerical framework is that of the Material Point Method (MPM). A set of particles (material points) move within a fixed finite element computational grid. Material points carry all the information (density, velocity, stress, strain, external loads…) which, at each step, is transferred to the grid to solve the mechanical problem. The computed solution allows updating of position and properties of the material points. The application of MPM to geotechnical problems is relatively recent [20,21]. Despite that, several implementations of MPM have been already used to model rigid body insertion into soils [22,23].
The Particle Finite Element Method (PFEM) is a third continuumbased approach that seems suitable to address geotechnical insertion problems. PFEM is actually an updated Lagrangian approach, but one that avoids mesh distortion problems by frequent remeshing. The nodes discretizing the analysis domain are treated as material particles the motion of which is tracked during the numerical solution. Remeshing in PFEM is based in Delaunay tessellations and uses loworder elements. PFEM was first developed to solve fluidstructure interaction problems [24,25,26,27,28] and then extended to other areas, like solidsolid interaction and thermoplastic problems [29,30], erosion problems [26], and Binghamlike rheology models to simulate flowslides [31].
Carbonell et al. [32,33] first applied PFEM to geomechanical problems, extending the method to deal with toolrock interaction problems in small and large scale excavations. In their work, however, material removal at the interface, rather than tool insertion was the focus. The excavated material was treated as a singlephase solid using a damage law as a constitutive model.
Zhang et al. [34,35,36] presented a new PFEM implementation for granular flow applications, using a variational theorem to discretize the governing equations. A singlephase rigid plastic constitutive description of the soil was employed. An example of pipeline insertion into a Tresca soil was presented in [34], however most applications have focused on soil flow problems [35,36].
This work focus on developing a numerical framework capable of simulating insertion problems found in Geotechnical Engineering.
The starting point of this work is an implementation of the Particle Finite Element method (PFEM) in the numerical platform Kratos Multiphysics (Kratos Multiphysics; Carbonell et al., 2010, 2013, 2015; Dadvand et al, 2010; Rodriguez et al, 2016, 2017,a,b); this code is able to simulate the contact between multiple deformable bodies employing classical Solid Mechanics formulations and constitutive equations describing metals. As such, from a numerical point of view, the basic objectives of this work are:
Regarding geotechnical knowledge:
This work is organized in 10 chapters, from which the first and last one are dedicated to the introduction and conclusions, respectively. Chapter 2 presents a general overview of the computational framework; several numerical novelties are presented from Chapter 3 to 5. The geotechnical applications of the developed numerical scheme are presented in Chapters 6 to 9.
In particular, this work is organized as follows:
Chapter 2 discusses the basis of the numerical model employed on the present work. First, the basic features Particle Finite Element method (PFEM) are described; additionally, PFEM is compared with other wellestablished codes used in Geotechnical Engineering. After briefly describing some basic results of nonlinear Solid Mechanics, the balance equations for the onephase problem and coupledhydromechanical simulations are presented. The finite element method is applied to discretize the spatial domain employing an Updated Lagrangian formulation whereas completely implicit timemarching schemes are used.
Chapter 3 is devoted to the numerical treatment of constitutive models at large strains. After revising the large strains elastoplastic theory based on a multiplicative split of the deformation in an elastic and plastic part, a novel explicit scheme is proposed. This chapter also presents the constitutive models employed in this work.
The numerical procedures to impose the contact between a deformable body and a rigid object are described in Chapter 4. Contact constraints are imposed to the solution by the Penalty method. To describe the tangential part of the contact model the socalled elastoplastic analogy is used; twodifferent algorithms to integrate this model are employed: an implicit algorithm that has the same formal structure than the onedimensional return mapping and one based on the Implex algorithm.
The low order finite elements typically employed in PFEM suffer from severe volumetric locking when the material shows a quasiincompressible response. In Chapter 5 several mixed stabilized formulations to mitigate locking for the onephase mechanical problem and the coupled hydromechanical problem at large strains are presented.
Chapter 6 presents the first application, the totalstress analysis soilstructure insertion problems in geomechanics. In particular, this chapter presents several analysis of increasing complexity; namely, the insertion of a rigid strip footing, the TBar and a rough cone penetration test. These three problems have been frequently as benchmark problems to assess the robustness and accuracy of largestrain geotechnical codes; thus, allows to compare the performance of the developed numerical scheme with other techniques.
Chapter 7 describes a series of simulations of the tube sampling process in clay materials using a total stress approach. In particular, a parametric analysis of the problem, in which several sampler geometries, constitutive parameters and interface behaviors are analyzed. Additionally, the occurrence of a plug is assessed in terms of the theory proposed by [1].
Chapter 8 shows the possibilities of the developed numerical scheme to simulate insertion problems unsing a coupled hydromechanical approach. The scheme is first assessed in the loading and dissipation of a footing resting in a poroelastic medium. Afterwards, results of the most challenging cone penetration test in a wide range of hydraulic conditions (ranging from practically undrained to drained) and using (tip and shaft) rough interfaces are presented. Details of the effect of these two parameters to the CPTu reactions (cone resistance, water pressure at the , and sleeve resistance) are given.
Several methods have been proposed to obtain the permeability or the coefficient of consolidation from CPTu readings. In Chapter 9, simulation outputs obtained for different input constitutive parameters and permeabilities are examined to obtain direct estimates of permeability using different methods proposed in the literature; additionally, methods to estimate the hydraulic conductivity during the piezocone penetration are also used. These estimates are then compared with the known input permeability value to assess their reliability.
Finally, in Chapter 10, the main outputs of this work are summarized, a number of conclusions are drawn and ideas for future enhancements and developments are enumerated.
This monograph also includes several appendices:
Appendix A presents a Matlab file to perform symbolic operations. In particular, this code is used to obtain an estimate for the stabilization parameter for the onedimensional element with linear (equalorder) shape functions.
Appendix B analyzes the effect of the input values of the [42] hyperelastic model and later modified by [43] on the material response.
The linearization of the mixed formulations introduced in Chapter 5 are presented in Appendix C.
The analysis of closedended piles and the bearing capacity factor of tubes, a byproduct of Chapter 7, are presented in Appendix D.
Some of the simulations of Chapter 8 are recalculated in Appendix E using slightly different input values of the hyperelastic model to assess its effect in a boundary value problem. The sets of employed constitutive parameters are first described, from a elementary point of view, in Appendix B.
Finally, in Appendix F, results of threedimensional analysis are reported. In addition, in this appendix, the hyperelastic plastic and hypoelastic plastic constitutive models are compared in a number of representative numerical simulations.
This monograph, which contains some unpublished material, is based on the following publications wherein some ideas and figures have appeared previously.
Publications in journals:
Pubications in conference proceedings:
This chapter is devoted to present the basics of the numerical model. First, the Particle Finite Element Method, the computational method used in this work, is presented. Subsequently, some basic definitions and results of continuum mechanics that will be latter used for the development of the governing equations are reviewed. Finally, the finite element formulation for the singlephase and the coupled fluidsaturated porous media at large strain is presented.
The Particle Finite Element Method (PFEM) is a numerical method well suited for mechanical problems involving large displacements, large deformations, intermittent separation and/or fusion of bodies. Soft porous materials, such as soils, suffer these kinds of mechanical transformations during many activities of engineering interest. Relevant examples for the case of soils include probing, sampling, pile installation, excavation and ploughing.
PFEM originated to address problems of fluid mechanics [24,25,26,27,28], including those of fluid interaction with rigid bodies. Later it was extended to deformable singlephase solids [29,44] and to the contact between multiple deformable bodies in thermoplastic problems [30,40,41]. Subsequently, several PFEM extensions have addressed geomaterials: Carbonell et al. [32,33] used PFEM to simulate ground excavation problems. The interaction between soil and structures has been also simulated using hypoplastic formulations to describe the constitutive behavior of soft soils [45]. Flowlike landslides have been also studied using PFEM, but considering a singlephase material description: [36] employed a rigid plastic constitutive response for the soil whereas [31] used a nonNewtonian modified Bingham law. Larese et al. [46] presented a strategy to simulate the free surface flow over and throughout a rockfill; PFEM is adopted for the evaluation of the structural response, whereas an Eulerian fixedmesh approach is employed for the fluid.
In PFEM the continuum is modelled using an Updated Lagrangian formulation; that is, a Lagrangian description of the motion is used and all variables are referred to the last known configuration. A mesh discretization of the domain is generated in order to solve the governing equations in the standard FEM fashion. Nodes in that mesh are treated as material particles whose motion is tracked during the solution.
The quality of the numerical solution depends on the discretization chosen. The original idea of the PFEM was to improve the mesh quality by performing a retriangulation of the domain only when needed. Usually that is performed according to some criteria associated with element distortion. Mesh distortion is corrected and improved naturally with the Particle Finite Element Method (PFEM), because retriangulation is based on Delaunay tessellations that maximize the minimum angle of all the triangles in the tessellation. Therefore, thin, stretched elements are avoided while still capturing large changes in the continuum domain without global remeshing and mesh to mesh interpolation. The process can be easily extended to 3D using tetrahedral elements.
Additionally, hadaptive techniques are employed to obtain a better discretization of the domain. New particles are introduced in areas where large gradients in the flow variables are detected or where a high plastic dissipation is generated. These zones must be refined because the number of particles may become too low to obtain an accurate solution. On the contrary, due to high shear deformations, particles may locally concentrate in the same region of the domain. To overcome this difficulty particles that are closer than a characteristic distance are removed.
Conforming meshes are used to preserve mass in the remeshing process: the boundary of the domain remains unchanged so the volume of the discretization does not vary. The mass of the domain also depends on the density and thus on the transfer process of Gauss point variables. This, however, is relatively unimportant for soil applications, since highest mesh distortion is usually associated with shearing and the attainment of incompressible critical state conditions.
Although it is not strictly necessary (see [34,35,36]), low order finite elements are typically used in PFEM: linear triangles in twodimensional models and linear tetrahedron in three dimensions. Linear interpolated elements have several advantages due to their simplicity: particles usually define exclusively the mesh nodes and no extra interpolations are needed after remeshing. Also, the computational cost is lower compared to highorder elements. A relative drawback is that stabilized mixed formulations are generally required instead.
The interpolation of state variables plays a crucial role in the accuracy of the results. When new particles are inserted in the domain, flow variables are linearly interpolated from those of the previous mesh element. To avoid excessive smoothing of internal variables, information is transferred directly from the previous Gauss points to the new ones. In this work, two different interpolation procedures for Gauss point variables are used. The first method is a nearest interpolation procedure; hence, new integration points inherit the information of the closes Gauss point of the previous mesh. In the second method, internal variables are interpolated over the whole domain and the values at the new integration points are obtained by employing a least square approach [47,48]. Both strategies ensure that information is not altered in elements that are preserved during the meshing process. These strategies are described in Section 2.1.1.
A typical solution algorithm of PFEM is conceptually illustrated in Figure 1 for a general case of fragmentation and deformation of a solid mass under the action of an external object; note, however, that fracture is not considered in this work. For clarity purposes we will define the collection or cloud of nodes (C) belonging to the analysis domain, the volume (V) defining the analysis domain and the mesh (M) discretizing the domain.
The simulation involves the following steps [24]:
Further details of the numerical implementation can be found in [32,33,38] and [30,40,41].
PFEM has some similitudes with other well established methods used in geomechanics, particularly variants of the Arbitrary LagrangianEulerian approach: the socalled Efficient Aribtrary LagrangianEulerian (EALE) developed by [12,13] and the remeshing and interpolation technique by small strain (RITSS) [9,10]. A fundamental difference with EALE lies on mesh treatment. On one hand, in EALE, the number of nodes, elements and the topology of the Finite Element mesh are preserved during the analysis: that is, between solution steps the bounday and interior nodes are recolocated computing a complementary elastic problem [12]. On the contrary, the original idea of PFEM is to minimize nodal position changes during computation while constantly updating mesh topology using a Delaunnay's Tesselation. Later, to reduce the dependence of the solution on the initial discretization, additional adaptive techniqueds have been introduced: (i) insertion and removal of nodes based on plastic dissipation and (ii) improving the mesh quality through Laplacian smoothing [33,30]. Despite these new features, modification of original nodal position in PFEM remains relatively infrequent.
The remeshing aspects of PFEM makes it closer to RITTS, particularly in its first implementations [9,50]. These periodical remeshing was also taking place using retriangulation, adaptive techniques and mesh smoothing. The algorithmic details are, however, different, particularly in respect to retriangulation. The degree of adaptivity of some newer versions of RITSS, that rely on Abaqus for mesh generation is perhaps more limited, as the role of users experience was emphasized by [19]. The algorithms employed for transfer of informaiton between successive meshes in PFEM are also different from those employed in RITSS.
Another significant difference is that whereas both RITSS and EALE use highorder elements to solve the governing equations (quadratic elements in RITSS and quadratic and higher in EALE); in PFEM linear elements are typically used. This, however, is not the only possible choice: Zhang et al. [34,35,36] used mixed elements with high order interpolants.
Other differences between the code presented here and previous approaches may also be a matter of choice, for instance early implementation of RITSS relied on small strain finite element formulations, but there is no major obstacle for it.
A comparative review of these two ALE methods along with the Coupled EulerianLagrangian has been recently presented by [19]. In a set of benchmark problems, the results of RITSS and EALE agreed very well with each other, suggesting that differences adopted in the remeshing part of the code provide similar accuracy.
One of the main difficulties encountered in adaptive methods is the transfer of information between different discretizations [51,52]. In PFEM information is transferred between the previous mesh and the newly constructed one at every remeshing stage. Nodal variables (displacements, mean pressure, water pressure, ...) are mapped into the new mesh using the previous mesh shape functions.
Note that for the relatively simple constitutive models employed here (that are presented in Chapter 3), the only internal variables that need mapping are the elastic left Cauchy Green tensor that, through the hyperelastic model uniquely determines the Cauchy Stress tensor and plastic history variables.
Two different transfer operators have been used. The algorithms selected share a common trait: when an element of the old mesh is preserved, the new value of the variable coincides with the value on the previous mesh. The first method consists simply on copying to each element of the new mesh the information of the element of the previous mesh whose centroid is nearest to the centroid of the new element. In the second strategy, a least square interpolation procedure is used [47]: let and be some internal variable in the new and previous mesh; a piecewice constat interpolation of these variables over the whole domain is constructed, and , where the interpolation functions and are equal to the unity in elements and of the new and old mesh and zero elsewhere.
The value of the internal variable in the new mesh is computed solving the following problem:

(2.1) 
After carrying out the minimization, the following explicit expression is found for the new value of the internal variable for one integration point elements:

(2.2) 
Therefore, the value of the internal variable on an element of the new mesh is the mean value of the variable of the elements of the previous mesh that overlapped that position, averaged by the area of overlapping. This algorithm is implemented using the supermesh concept developed by [48].
Only the first transfer algorithm (nearest point interpolation) ensures that the new state is admissible (using the second algorithm, the stress state may be outside of the yield surface); none of the two guarantee that the deformable body is in equilibrium. The first problem is directly tackled applying a yield surface drift correction algorithm if the stress state lies outside the yield surface. Possible errors due to outofbalance forces after remeshing are ignored and a timestep is advanced before mechanical equilibrium is again imposed.
These two algorithms are numerically assessed and compared in Section 6.2, that presents the results of the simulation of a rigid strip footing on clay employing a total stress approach. It is found that effect of the interpolation technique on the results is limited; thus, the nearest element interpolation is selected in the majority of the simulations due to its relatively smaller computational cost.
The objective of this section is to review some basic definitions and results of nonlinear continuum mechanics needed for the development of the formulation of the balance equations for the singlephase mechanical problem, the coupled hydromechanical problem and then the formulation of the elastoplastic constitutive equations. Extended details may be found in [53,54,55], among others.
Let define the reference configuration of a continuum body and a particle of this body is labeled by . The position of this particle at time is:

(2.3) 
where is the mapping from the initial configuration to the current one. Let us define now the displacement of a material point as the difference between the current position and its original position. The displacement vector, , can be written as:

(2.4) 
The material velocity and material acceleration, denoted by and respectively, are defined as:

(2.5) 

(2.6) 
Meanwhile, a spatial or Eulerian description is also possible:

(2.7) 

(2.8) 
where the first term is known as the local derivative and the second one is the convective time derivative.
Figure 2: Sketch of the reference and current configuration of a deformable body [56]. 
The deformation gradient, , is the Jacobian matrix of the motion:

(2.9) 
To further illustrate the deformation gradient, let us consider two different particles, and in the reference configuration, whose infinitesimal relative position in the reference configuration is , see Figure 2. Then, the relative position in the deformed configuration, , may be obtained as:

(2.10) 
Furthermore, a differential volume in the reference configuration, , is related to the same differential volume in the deformed configuration, , as:

(2.11) 
where the Jacobian determinant is frequently referred to as the Jacobian: .
Considering a differential of area in the reference configuration, , where is the unit normal, the following expression for the transformation is given by Nanson's formula:

(2.12) 
Several deformation measures are used in nonlinear continuum mechanics. Two of the most common ones are the right CauchyGreen tensor, , and the left CauchyGreen tensor, :

(2.13) 

(2.14) 
where the first one is defined in the reference configuration whereas the second one is defined in the deformed configuration. In the initial state or in the absence of deformation, , both deformation measures are equal to the identity tensor.
According to the polar decomposition theorem, the deformation gradient may be decomposed multiplicatively into a stretching and a rotational part:

(2.15) 
where is a proper orthogonal tensor, called the rotational tensor, and and are two symmetric positive definite tensors. These tensors fulfill the following properties:

(2.16) 
where and are, respectively, the right and left CauchyGreen tensors. The eigenvalues of and are coincident, whereas its eigenvectors differ:

(2.17) 
where are the eigenvalues of and , also known as principal streches, whereas and are the eigenvectors of and respectively.
Taking into account the decomposition of the previous tensors, both CauchyGreen tensors may be redefined as:

(2.18) 
Let us start by computing the rate of the deformation gradient, Equation 2.9:

(2.19) 
that is, the time derivative of the deformation gradient is equal to the material velocity gradient, .
Applying the chain rule to the previous expression:

(2.20) 
then, the spatial velocity gradient, , may be obtained as:

(2.21) 
In the development of finite element stiffness matrices and in the constitutive equations, the temporal derivative of several strain measures will be employed. For instance, the temporal derivative of the right CauchyGreen tensor may be useful:

(2.22) 
where is the spatial rate of deformation tensor.
The temporal variation of the Jacobian determinant is:

(2.23) 
To define the GreenLagrange strain and the Almansi strain, let us compute the variation of length of the segment of Figure 2:

(2.24) 
From this previous expressions, the GreenLagrange strain, , and the Almansi strain, , are defined as:

(2.25) 
In terms of the eigenvectors and eigenvalues of the left and right CauchyGreen tensors, Equation 2.18, these strains may be obtained as:

(2.26) 
The unit elongation along a direction or is:

(2.27) 
Consequently, the unit elongation may be further related to the GreenLagrange and Almansi strain tensors through the following expression:

(2.28) 
The definition of the stretch has already appeared without a proper definition, for instance in Equation 2.17. The stretch along an arbitrary direction, , is related to the unit elongation through:

(2.29) 
Another popular strain measure is the (spatial) Hencky strain:

(2.30) 
(a)  (b) 
Figure 3: Comparison of the Almansi, GreenLagrange and Hencky strain in terms of the unit elongation in the small strain range, (a), and at finite strain range, (b). 
In order to compare all these strain measures, Figure 3 depicts the dependency of these strain measures in terms of the unit elongation for a onedimensional medium. In the small strain regime (unit elongations below 10 %), all these strains agree well and predict similar strain values (Figure 3(a)). However, discrepancies appear at the large deformation regime: as noted by its definition, Equation 2.26, the maximum value in extension of the Almansi strain is 0.5. On the other hand, the maximum value in compression of the GreenLagrange strain is 0.5. Potentially, the Hencky strain may take any real value.
The former observations are restricted to onedimensional cases. It should be noted that in a general threedimensional continuum, the Almansi and the Hencky strain share eigenvectors, that are different from those of the GreenLagrange tensor. By definition, the GreenLagrange tensor assumes a material description and, thus, acts to an element defined in the reference configuration. On the contrary, the other two strain tensors have an spatial description and acts to an element of the deformed configuration [56].
There exist a large number of stress measures. The most intuitive one is the Cauchy stress tensor, , that is defined in the current deformed configuration . The relation between the traction vector, , and the surface normal vector, , is:

(2.31) 
From the local balance of angular momentum, it can be obtained that the Cauchy stress tensor is symmetric.
The Kirchhoff stress tensor, , is also symmetric and defined in the current configuration. One reason for its use is that, in many equations, the Cauchy stress appears together with the Jacobian and the use of it simplifies formulas:

(2.32) 
To define the Second PiolaKirchhoff stress tensor, , let us express the force acting on a differential of area:

(2.33) 
performing a pullback of the force and introducing the Nanson's formula:

(2.34) 
which leads to the definition of the Second PiolaKirchhoff as:

(2.35) 
that is a symmetric stress tensor defined in the reference configuration.
The time derivative of stress tensors is of significance for the statement of incremental forms of constitutive equations. For the second PiolaKirchhoff, that is referred to the initial configuration, the derivative with respect to time is given by the material time derivative:

(2.36) 
However, the temporal derivative of stress tensors related to the deformed configuration are more complex. For instance, the derivative of the Cauchy stress tensor reads:

(2.37) 
This section is devoted to frame indifference and objective transformation. Although the discussion of which deformation measures and stress measures transform objectively seems quite arid and unmotivated, the importance of the previous discussion will become apparent in the chapter devoted to constitutive relations (Chapter 3), that need stress and strain measures independent of rigid body motions.
A spatial tensor field is said to transform objectively under superposed rigid body motions if it transforms according to the standard rules of tensor analysis [53].
Given a motion , a superposed rigid motion of the current place is the map defined as:

(2.38) 
where is a function of time and is a proper orthogonal transformation depending only on time. The superposed motion is called rigid because, for any given two points , , the following property holds

(2.39) 
in other words, the mapping preserves the distances.
Thus, the deformation gradient becomes:

(2.40) 
Then, the left CauhyGreen tensor, Equation 2.14, reads:

(2.41) 
which transforms objectively.
On the contrary, material objects remain unaltered to the superposed rigid body motions; for instance, the right CauchyGreen tensor, Equation 2.13:

(2.42) 
The spatial velocity gradient, Equation 2.21 is given by:

(2.43) 
which does not transform objectively because of the additional last term. However, its symmetric part, the rate of deformation tensor, transforms objectively since:

(2.44) 
Demonstrating that the Cauchy stress tensor is objective is quite straightforward: on physical grounds, it is natural to postulate that the tractions are objective, whereas unit outward normal is objective, . Then:

(2.45) 
it follows from the previous expression that . Since this equation holds for any normal, it follows that:

(2.46) 
that is, the Cauchy stress tensor transforms objectively.
However, the temporal derivative of the Cauchy stress tensor is not objective:

(2.47) 
that is clearly nonobjective since it is not equal to .
Several objective rates of the Cauchy stress and the Kirchhoff stress tensor that is, modified time derivatives constructed to preserve objectivity have been proposed in the literature. Among others, one that will appear repeatedly in this text is the Lie Derivative, that is defined as:

(2.48) 
using the expression of the derivative of the inverse and arranging some terms:

(2.49) 
In this section, the formulation for the singlephase mechanical problem is presented. The section starts by enunciating the governing equation the linear momentum balance equation in an Updated Lagrangian fashion. After briefly describing the weak form, the balance equation is discretized with low order shape functions and an implicit timestepping algorithm. The nonlinear system of equations is solved with a NewtonRaphson method; thus, the required linearization is derived.
A quasistatic linear momentum formulation in an updatedLagrangian form (i.e. expressing all quantities and their derivatives in the deformed configuration), may be written as:

(2.50) 
where is the Cauchy stress tensor, stands for the appropriate constitutive equation for path dependent materials (large strains elastoplastic constitutive equations based on the multiplicative split [57] are used here, see Chapter 3), is the total deformation gradient and represents the set of internal variables of the model. stands for the initial displacement, is the external body force vector, is the gravity vector and () defines the boundary of the domain where displacements and tractions are prescribed.
Meanwhile, the mass conservation principle requires that the mass of the continuous media domain remains constant. The mass balance equation in its Lagrangian form reads:

(2.51) 
Finally, the local balance of angular momentum yields that the Cauchy stress tensor is symmetric:

(2.52) 
The derivation of the weak form of the problem defined in Equation 2.50 starts by multiplying the local balance equation by a test function, , and integrating it in the domain:

(2.53) 
Integrating by parts the term related to the Cauchy stress gradient, applying the divergence theorem, introducing the Neumann boundary conditions and using index notation yield:

(2.54) 
This equation coincides, in a formal sense, with that obtained in the infinitessimal strain theory. However, using in the large deformation theory, stress and virtual strains have to be evaluated with respect to the current configuration. Additionally, it is integrated over the deformed configuration. As such, the nonlinearities do appear, however, they are hidden.
In order to obtain the finite element discrete equations of the weak form, first the nodal variables are approximated with the Finite Element shape functions

(2.55) 
where is the finite element approximation of displacement whereas are the nodal values. are the shape functions whereas is the number of nodes. The same procedure is used to discretize the test functions, .
As such, the problem is then: Find the function such that:

(2.56) 
Since are arbitrary, from the previous expression independent equations may be obtained, where is the number of spatial dimensions of the problem. As such, the matrix form of the Galerkin weak form is:

(2.57) 
where the internal and external forces read:

(2.58) 

(2.59) 
where has the same formal structure than the small deformation straindisplacement matrix [58] and corresponds to the Voigt notation of tensor .
As already mentioned, a quasistatic loading is assumed, so that inertial terms have been neglected. Then, the Galerkin form of the linear momentum balance equation, Equation 2.57, is time independent. Then, the global problem may be directly stated, for time , as: find such that:

(2.60) 
Equation 2.60 defines a nonlinear system of equations, whose residual is given by:

(2.61) 
where the superscript is obviated for simplicity. Then, to evaluate numerically the solution, an iterative nonlinear solver is required. In this work, the NewtonRaphson method is used. Writing an iterative correction as:

(2.62) 
where the subscript stands for the number of iteration. Using a Taylor's approximation of the residual, Equation 2.61,

(2.63) 
where is the socalled tangent or Jacobian matrix, that is defined as:

(2.64) 
Assuming that is null, the iterative correction, , is given by the following system of linear equations:

(2.65) 
It is worth noting that the term may be also obtained by using the directional or Gateaux derivative along the direction [54]. That is:

(2.66) 
Let us linearize the governing equations of the problem, Equation 2.60. However, instead of working on the discrete set of equations, which appear after the introduction of the finite element discretization, let us work with the continuum counterpart 2.54. In the case of elastic materials, the linearization of the discrete equations and the discretization of the linearization of the continuum equations are equivalent [57,54]. However, when inelastic materials are considered, both forms are no longer equivalent: in the first approach the material constitutive tensor depends upon the stress integration algorithm (the socalled consistent tangent matrix) whereas in the second approach the Jacobian matrix depends on the continuum constitutive tensor. As such, using the consistent constitutive matrix, second order convergence in the solution of the global problem may be achieved; thus, virtually representing a saving of computational time.
In order to linearize the discrete equation of the balance of linear momentum, Equation 2.54, let us first work with the internal forces term. First, a pullback of the internal forces term is performed; this way, all the quantities are referred to the reference configuration:

(2.67) 
Then:

(2.68) 
Then, a pushforward of this expression yields:

(2.69) 
where the directional derivative of the deformation gradient, , is expressed as:

(2.70) 
whereas the term related to the second PiolaKirchhoff may be related to the Lie derivative of the Kirchhoff stress and, also, to the spatial constitutive tensor:

(2.71) 
where is the spatial constitutive tensor.
So finally,

(2.72) 
where .
Surprisingly, the terms related to the external loads are typically not linearized; however, as it will be shown here, the term due to the prescribed tractions may induce a geometrical nonlinear term.

(2.73) 
The pullback of the term related to the gravitational loads reveals that there is no dependency on the displacements; that is

(2.74) 
Figure 4: Convective coordinates on a parametrized surface. [59]. 
On the other side, a pullback of the term related to the imposed tractions, making use of a parametrized description of the boundary, see Figure 4, reads:

(2.75) 
where . By using the chain rule, this term reads:

(2.76) 
According to [59], this directional derivative may be expressed as:

(2.77) 
where .
Finally, introducing the expression 2.77 to 2.76 and performing a pushforward:

(2.78) 
Roughly speaking, this term stands for the variation of the virtual work associated to the variation of area and direction of the contour where tractions are imposed. Although it has not been explicitly stated, it is assumed that the imposed traction is independent of the boundary. To exemplify this fact, let us assume that an external load is imposed in an area at the initial state; then, the imposed force, is equal to . Assuming that the body deforms but the normal to the surface where Neumann condition remains unaltered (), the area of the surface is and the total applied load at the Neumann boundary is . As such, the force in the Neumann boundary varies; since it has been assumed that the traction that is applied at the boundary is independent of the deformation of the body. Of course, this is a modeling decision and other options are possible.
Finally, the matrix form of these expressions is obtained by inserting the definition of the virtual dispalcements and the displacements, Equation 2.55:

(2.79) 
where the terms involed in the computation of the nonlinear stiffness matrix, and , may be described, for two dimensional analysis, as [60,54]:

(2.80) 

(2.81) 
where is the derivative of the local shape function 1 with respect to coordinate 1 and is the deformed radial coordinate. It should be noted that the expressions provided in Equations 2.80 and 2.81 are particularized for a twodimensional axisymmetric case. For a planestrain situation, only the first four columns in Equation 2.80 and the first four columns and rows in Equation 2.81 should be considered. The definition of these matrices for a threedimensional analysis may be found in [60].
The mechanics of porous media is relevant in a broad range of engineering applications. Several models have been proposed under the small strains assumption, that slightly differ in the treatment of the dynamic terms of the formulations, see, for instance, the socalled (solid displacementDarcy's velocitywater pressure) and the simplified of [61]. The governing equations have been obtained by a broad range of theories, from physical approach to mixture theory, but identical equations are found irrespectively of the method to derive them.
The extension of the poromechanical problem at large strains can be traced back to [62], that proposed a three dimensional finite element formulation for consolidation problems; finite strains constitutive models based on a Jaumann rate description were used.
Borja and Alarcón [63] developed a formulation based on the mixture theory combined with the volume fraction concept, where the constituents are considered incompressible. Constitutive models are described using a multiplicative decomposition of the deformation gradient in an elastic and plastic part.
An important extension of the formulation is due to [64], that introduced the fluid compressibility in terms of a volumetric logarithmic strain, a measure that is energy conjugated to the fluid pressure; then a free energy function for the fluid phase is defined. By doing so, the formulation is able to model the porous medium at nearly saturated conditions.
In this section, the computational framework for the analysis of fluidsaturated porous media is presented. First, the governing equations the linear momentum and mass balance equations of the mixture are derived at large strains. After discretizing the equations using the finite element method in space and a totally implicit timeintegration algorithm, the consistent linearization required for the nonlinear solver is presented. Finally, the methods to stabilize the formulation are highlighted; this includes a novel technique to estimate the value of the stabilization parameter for hydromechanical problems.
In this subsection the equation for the mass balance of the mixture is derived.
Let and be the density of the solid particles and the fluid, respectively, and the porosity in the current configuration. Then, the mass of solid and liquid, and , in an arbitrary deformed domain, , may be defined by the following integrals:

(2.82) 
Assuming that the mass is conserved along a material volume, then material derivative of the integral is null:

(2.83) 
where and are the material derivatives with respect to the solid and fluid phase and and are the solid and fluid velocities.
Note that in the previous expression the solid (or skeleton) phase is used as a reference frame for the porous media as a whole. As such, all the quantities are referenced to this configuration, leading to a modified Eulerian description of the fluid motion. Then, all the spatial gradients are taken with respect to the solid phase motion.
Since integration takes place in an arbitrary volume, rearranging terms and introducing the definition of the material derivative with respect to the solid phase, , the last expression may be recasted as:

(2.84) 
where the Darcy velocity has been introduced: .
Expanding the derivatives and dividing by the density, the following expressions may be obtained:

(2.85) 
Finally, the following expression for the mass conservation of the porous media is obtained summing both equations:

(2.86) 
In the rest of the document it is assumed that the solid particles (the grains) are incompressible:

(2.87) 
In order to define a constitutive model for the fluid, Larsson and Larsson [64] introduced a logarithmic strain for each constituent, that, for the case of the fluid constituent reads:

(2.88) 
where is a reference fluid density.
Based on thermomechanical considerations, Larsson and Larsson [64] proposed several constitutive models relating the (Cauchy) water pressure and the logarithmic strain. The simplest model is to assume that the water pressure is proportional to the logarithmic volumetric strain:

(2.89) 
where is the water bulk modulus.
As a consequence:

(2.90) 
Several extensions of the Darcy's Law at large strains have been proposed. For instance, Borja and Alarcón [63] proposed to postulate the Darcy's Laws as:

(2.91) 
where is the Kirchhoff water pressure, is the permeability tensor, is the gravity and . This is done since, in the referred work, the Kirchhoff water pressure is a nodal variable and because, accidentally, it is assumed that , which, in the general case, is not true.
On the other hand, in this work the expression proposed by [64] is used:

(2.92) 
where and is the permeability tensor with velocity units.
In order to enhance the discussion of the constitutive model for the permeability tensor, let us perform a Piola transformation of the Darcy's velocity [62,64]:

(2.93) 
where is the material gradient of the water pressure. As a consequence of the previous expression, we can relate the permeability in the reference configuration, , and in the current configuration as:
After defining the Darcy's Law and showing the relation between the spatial and material definition of the permeability tensor, let us discuss the constitutive models for the permeability tensor.
At this point, two different hypothesis may be performed: On the one hand the material permeability tensor, , may be considered constant; this hypothesis may be interesting for media with anisotropic permeability tensors since the tensor rotates with the deformation of the solid phase [64]. On the other hand, the most common approach is to assume that the spatial permeability tensor, , is constant. In the majority of the simulations, the second hypothesis is assumed.
Not only a constant permeability tensor (in the reference or the deformed configuration) are used in this work; the more sophisticated KozenyCarman Law is also used [65]. In this law, the permeability tensor depends on the volumetric strain of the biphasic porous media through the void ratio. Then, the spatial permeability tensor is defined as:

(2.94) 
where and are the permeability and void ratio at the initial state and is the void ratio, that is defined according to:

(2.95) 
It must be pointed out that the residual form of the balance mass is not affected on the assumption made on the permeability; however, the tangent matrix to the discrete finite element method depends on the assumptions made on the permeability constitutive tensor. However, these terms have not been fully derived in this work; for instance, when the KozenyCarman law is used, the terms that appear due to the dependency of the permeability on the volumetric strains are not considered.
Once the constitutive equations relevant to the fluid phase have been introduced, let us further elaborate the expression for the balance of the mass of the mixture, Equation 2.86. Introducing the definition of the water phase compressibility (Equation 2.90) and assuming that the solid phase is incompressible, the following expression is found:

(2.96) 
Or, alternatively, introducing the definition of the Darcy's Law, Equation 2.92, and expressing the gradient of the water density in terms of the water pressure:

(2.97) 
It is assumed that the term , since very large water bulk modulus will be used, so the spatial gradient of the density is nearly zero. Additionally, it is assumed that the term is almost constant. Finally, the mass balance equation of the mixture reads:

(2.98) 
The local form, using an spatial description of the domain, of the balance of linear momentum equation reads:

(2.99) 
where is the mixture density and stands for the total Cauchy stress.
According to the principle of effective stress, the total Cauchy stress tensor, , is equal to the sum of the pore pressure, , and the effective stress:

(2.100) 
where the effective Cauchy stress depends on the solid skeleton deformation through the deformation gradient, , and a set of history dependent parameters, .
The porosity may be further related to the solid skeleton motion. In an infinitesimal volume of the mixture, , the initial void volume is whereas the volume of the solid fracion is . As deformation takes place, the volume of the solid matrix varies according to . Since the mass of solid particles is the same in the deformed configuration and the solid phase is assumed incompressible, the volume of the solid fraction is whereas the volume of the voids is . As such, the porosity may be described by the following expression:

(2.101) 
Then, the mixture density, , can be further elaborated to:

(2.102) 
Let us state the system of equations governing the mechanics of the byphasic porous media in an updatedLagrangian form for quasistatic cases along with the corresponding boundary conditions:

(2.103) 
where is the material time derivative with respect to the solid phase. The boundary of the domain is divided in two parts, (), where fixed water pressure and prescribed water flow are imposed.
The weak form of the linear momentum balance equation reads:

(2.104) 
Integrating by parts the term related to the effective Cauchy stress divergence and water pressure gradient, applying the divergence theorem, introducing the Neumann boundary conditions and using index notation yields:

(2.105) 
where is the second order identity tensor in index notation.
The weak form of the mass balance equation, Equation 2.98, reads:

(2.106) 
where is the virtual water pressure.
Note that integration takes place over the reference domain; this is not the only possibility, but it has been done this way because, as it will be shown later, some of the finite element matrices are constant and it keeps the same rational as most mixed formulations for the onephase problem (the socalled formulations). On the other hand, Borja and Alarcón [63] integrate the mass balance equation directly over the current configuration (in other words, Equation 2.106 is multiplied by the Jacobian, ), whereas Larsson and Larsson [64] state the mass balance equation in terms of fluid content (i.e, Equation 2.106 scaled by ).
Again, the divergence theorem is applied to the weak form, that yields:

(2.107) 
After obtaining the weak form of the balance equations, Equations 2.105 and 2.107, the semidiscrete equations of the hydromechanical formulation are obtained. First, let us introduce the interpolants:

(2.108) 
where is the finite element approximation of the field whereas are the nodal values. and are the shape functions whereas is the number of nodes. As it can be inferred, the same order shape functions are used for both, solid phase displacements, , and water pressure, . Since the shape functions are defined in the reference configuration and, thus, time independent, , the following property holds: .
The semidiscrete equations of the hydromechanical formulation, Equations 2.105 and 2.107, are given by:

(2.109) 
where the definition of the Darcy's law has been introduced and the matrices and vectors of the previous expression are defined as:

(2.110) 

(2.111) 

(2.112) 

(2.113) 

(2.114) 
Since quasistatic conditions are assumed, the temporal discretization of the linear momentum balance equation is quite straightforward. However, the mass balance has terms related to the temporal derivative of the displacements (the velocity) and the water pressure. A completely implicit integration scheme is used; then:

(2.115) 
Then, introducing the temporal discretization to the semidiscrete equations, Equation 2.109, the problem reads:

(2.116) 
Both balance equations, defined in Equation 2.116, are solved in a monolithical approach. Then, the residual takes the form:

(2.117) 
where:

(2.118) 

(2.119) 
Then, applying the NewtonRaphson nonlinear solver, the linear system of equations equivalent to Equation 2.65 has the form:

(2.120) 
As in the previous case, the linearization is performed for the continuum equations and then the spatial discretization will be introduced to obtain the discrete stiffness matrices. In this case, the balance equations with the temporal discretization already introduced are linearized:

(2.121) 
where the absence of subscript at displacements and water pressure stands for the value at time , and .
In order to obtain the linearization of the balance of linear momentum equation, let us first write all the derivatives in terms of the reference configuration and introduce the expression of the density of the mixture, Equation 2.102:

(2.122) 
where the term due to the surface tractions has not been included since its linearization has been presented previously.
First, let us perform the linearization with respect to displacements. Going term by term, the linearization of the term due to the effective second PiolaKirchhoff, , is quite straightforward:

(2.123) 
this term reflects the material constitutive model response. On the other hand, to obtain a term with the same formal structure than the geometrical stiffness matrix of the singlephase formulation let us now compute the derivative due to the last deformation gradient that appears in the integral of the internal forces:

(2.124) 
The only two terms that are still missing in the linearization are those related to the Jacobian and to the inverse of the Right CauchyGreen tensor:

(2.125) 
similarly to the temporal derivative of the right CauchyGreen tensor, Equation 2.22:

(2.126) 
as a consequence:

(2.127) 
Then, introducing this previous expression and the linearization of the Jacobian, :

(2.128) 
On the other hand, differently from the singlephase problem, the external load term due to the gravity also produces a second order term due to the variation of the mixture density:

(2.129) 
where it has been assumed that the water density is almost constant: that is, the term should be small since it is divided by the water bulk modulus, that is very large in comparison with all the other parameters.
Finally, the linearization with respect to the water pressure is quite straightforward:

(2.130) 
So, finally, the linearization of the linear momentum balance equation, defined in Equation 2.118, reads:

(2.131) 
In order to calculate the linearization of the mass balance equation, let us first perform a pushback of Equation 2.121, so that all the quantities are referred to the initial configuration:

First, let us calculate the linearization with respect to displacements:

where the first and last term of the first line of the previous equation are independent of displacements and, consequently, its derivative is null. Introducing the definition of the Darcy's velocity:

(2.132) 
First, let us work on the linearization of the first term of the previous equation:

(2.133) 
The inverse of a matrix derivative may be computed as: . Then, introducing this result to Equation 2.133 and performing a pushforward:

(2.134) 
Once the linearization of the first term has been presented, let us proceed with the second term of Equation 2.132:

(2.135) 
where, as commented previously, it has been assumed that the spatial permeability tensor, , is constant; other assumptions made on the permeability tensor may induce different large displacement terms. So, the following expression may be obtained:

(2.136) 
With respect to the water pressure, the same result is obtained either from the Updated or Total Lagrangian expression of the balance equation; then:

which has the same formal structure than the small strains counterpart.
So, finally, let us state the final form of the linearization of the mass balance equation, defined in Equation 2.119:

(2.137) 
where the first order terms have a similar formal structure than the small strains counterpart (except for the inverse of the Jacobian) and additional geometrical terms appear.
Finally, the matrix expression of the terms defined in Equation 2.120 may be obtained as:

(2.138) 

(2.139) 

(2.140) 

(2.141) 
where is the Voigt notation of the tensor , and the definition of terms involved in the computation of the nonlinear stiffness matrix, and , may be found in Equations 2.80 and 2.81 for twodimensional problems and, for threedimensional cases, in [60]. Due to the complexity, the geometrical terms of the mass balance equation that appear in Equation 2.140, , have not been obtained in matrix format; its expression using indicial notation is presented in Equation 2.137. In the code, these terms have been implemented using indicial notation.
Undrained conditions in watersaturated soils with nearly incompressible constituents results in quasiincompressible behavior. Therefore, as the problem approaches zero permeability and incompressilibity of the soil and water constituents, the system of the discretized equations describing the formulation has the similar structure to that found when using a mixed formulation of Solid Mechanics problems [66]. Depending on the interpolants used to discretize both nodal variables, the problem may become illposed from a mathematical point of view, which may result in nonuniqueness and meshdependence of the solution, the wellknown volumetric locking.
Volumetric locking introduce numerical stiffening and spurious high spatial variability in the solution, eventually leading to numerical instability The reason behind this behavior is the noncompliance of the BabuskaBrezzy conditions or the patch test due to an improper finitedimensional space in the finite element discretization.
It is well known that in the hydromechanical problem, the pore pressure field tends to exhibit oscillations, which tend to increase when the time step is reduced [67]. As shown in the onedimensional (Terzhagi's equation) analysis of the discrete equations, the use of implicit timeintegration scheme produce a stiffness matrix whose eigenvalues are positive and, consequently, the overall algorithm is unconditionally stable. However, the water pressure field may exhibit oscillations if the time step is smaller than:

(2.142) 
where is the coefficient of consolidation and the element size. By numerical examples, it can be shown that this limit also applies to two and threedimensional analysis [67].
To avoid this problem two strategies are common: either to use more complex, but stable, finite elements with different order interpolation of displacement and water pressure fields, or to apply stabilization procedures to originally unstable finite elements [68].
Since equalorder, loworder finite element shape functions are used to discretize both, displacements and water pressure, the only possible technique that could be used to mitigate volumetric locking is the use of stabilization techniques.
In this work, two simple stabilization techniques are used: the Polynomial Pressure Projection [69] and the Fluid Pressure Laplacian [58]. Both techniques have been used in the literature [70,71,72]. By using these techniques, the discrete governing equations are modified by introducing a new term in the mass balance equation 2.116:

(2.143) 
where is the stabilization matrix, whose definition depends on the employed stabilization technique.
The Polynomial Pressure Projection (PPP) stabilization method has been originally developed to stabilize Stokes equations [73,69] and has also been applied to stabilize the mixed formulation of Solid Mechanics (among others, [30]) and in Soil Mechanics to stabilize formulations similar to the one presented here [70,72].
The PPP has two main ingredients:
The method is obtained by modifying the mixed variational equation (i.e. the pressure continuity equation) by using local polynomial pressure projections of the pressure variable. The application of the projections in conjunction with minimization of the problem field mismatch, eliminates the inconsistency of equalorder approximations and leads to a stable variational formulation. Unlike other stabilization methods, the Polynomial Pressure Projection does not require the calculation of higherorder derivatives. It uses a projection on a discontinuous space and, as a consequence, can be implemented at element level avoiding the need of mesh dependent parameters. The implementation of this stabilization scheme reduces to a simple modification of the weak continuity equation (the incompressibility constraint).
Given a function , the projection operator is defined by

(2.144) 
where is the best approximation of the pressure in the space of polynomials of order .
Then, the stabilization term reads:

(2.145) 
where is the stabilization parameter.
Finally, the discrete stabilization matrix, , may be expressed:

(2.146) 
where are the set of polynomials introduced in Equation 2.144; in the case of linear triangles, these local element polynomials are .
In stabilization techniques, there exist a free parameter, that depends on the size of the element size, material constitutive parameters and the time increment [58]. However, obtaining an estimate for this parameter can be very complex. In this work, to estimate the stabilization factor the same technique that was proposed by [74] to obtain the critical timestep for implicit methods is used.
The rational of the method proposed by [74] is quite simple and is based on the solution of the onedimensional consolidation problem by using a finite element mesh whose nodes are equally spaced. It is assumed that all the nodes of the mesh have the same initial water pressure; at one extreme of the mesh the water pressure is increased whereas the other extreme of the mesh has null Neumann boundary conditions. In order to avoid oscillations on the computed water pressure, the solution has to fulfill two conditions:
By using this rational, the authors obtained that the critical time step should be:

(2.147) 
where is the coefficient of consolidation, is the constrained modulus and is the permeability in velocity units, whereas
It should be noted that this time constraint is equivalent to the one presented by [67] using a similar method.
Assuming small strains, one dimensional conditions, null boundary loads, incompressibility of the water and linear elasticity, the problem statement, Equation 2.103, reduces to:
