## Abstract

### Insertion problems in Geomechanics with the Particle Finite Element method

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 large-strain, elasto-plastic constitutive models, novel integration methods for rigid body contact constraints and stabilized mixed formulations for single-phase and two-phase problems.

An explicit stress integration scheme is developed for elasto-plastic, large-strain 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 pre-specified shape and motion may be modeled. Contact constraints are introduced using a penalty method in which tangential behavior is treated with an elasto-plastic analogy. This elasto-plastic contact model is integrated by means of an implicit integration technique; an alternative scheme using the Implex algorithm is also described and assessed.

Low-order (linear) elements are employed to speed computation. These elements may suffer volumetric locking in quasi-incompressible 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 three-field mixed formulations for the coupled hydro-mechanical 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 single-phase problem, which approximates saturated soil behavior under undrained conditions.

The first application examples are three total stress problems: indentation of a rigid strip footing, T-Bar 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 (round-tipped 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 open-ended 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 closed-ended pile. Indeed, the failure mechanism that prevails during the penetration of a plugged tube is coincident with that of a closed-ended pile.

The last analysis of this work is the hydro-mechanical 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 on-the-fly techniques, in which permeability could be directly estimated from the CPTu data stream without the need for any stoppage.

## Preface

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 BIA-2011-27217 and BIA-2014-59467-R is gratefully appreciated. The support provided by the European Commission funded GEO-RAMP RISE project (H2020-645665- GEO-RAMP) is also acknowledged.

# 1 Introduction

## 1.1 Introduction

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-, spud-can 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 in-situ 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 hydro-mechanical 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 non-linearities, contact-related, material-related and also geometrical. This latter kind of non-linearity 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, continuum-based approaches are dominant, particularly for fine-grained soils. Within continuum-based methods the approach most frequently applied to geotechnical insertion problems has been that of Arbitrary Lagrangian-Eulerian 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 so-called remeshing and interpolation technique by small strain (RITSS) developed by Randolph and co-workers [9,10,11], the so-called efficient ALE approach (EALE) developed by Nazem and co-workers [12,13,14] and the successive built-in implementation of ALE in Abaqus/Explicit, currently known as the Coupled Eulerian-Lagrangian (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 continuum-based 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 continuum-based 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 low-order elements. PFEM was first developed to solve fluid-structure interaction problems [24,25,26,27,28] and then extended to other areas, like solid-solid interaction and thermo-plastic problems [29,30], erosion problems [26], and Bingham-like rheology models to simulate flowslides [31].

Carbonell et al. [32,33] first applied PFEM to geomechanical problems, extending the method to deal with tool-rock 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 single-phase 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 single-phase 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].

## 1.2 Aim and objectives

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:

• Implement a coupled hydro-mechanical formulation at large strain.
• Develop and implement a robust and accurate explicit stress integration algorithm for elasto-plastic models at finite strain.
• Incorporate geotechnical interface models to the contact algorithms.
• Mitigate the volumetric locking of low-order elements for the one-phase and coupled hydromechanical formulations.

Regarding geotechnical knowledge:

• Assess the effect of the cutting shoe-geometry, wall thickness and contact roughness in the deformation path of tube sampling problems.
• Enhance the knowledge on the cone penetration test by analyzing the influence of permeability and contact roughness on the cone readings and subsequent dissipation.

## 1.3 Outline of the monograph

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 well-established codes used in Geotechnical Engineering. After briefly describing some basic results of non-linear Solid Mechanics, the balance equations for the one-phase problem and coupled-hydromechanical simulations are presented. The finite element method is applied to discretize the spatial domain employing an Updated Lagrangian formulation whereas completely implicit time-marching schemes are used.

Chapter 3 is devoted to the numerical treatment of constitutive models at large strains. After revising the large strains elasto-plastic 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 so-called elasto-plastic analogy is used; two-different algorithms to integrate this model are employed: an implicit algorithm that has the same formal structure than the one-dimensional 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 quasi-incompressible response. In Chapter 5 several mixed stabilized formulations to mitigate locking for the one-phase mechanical problem and the coupled hydro-mechanical problem at large strains are presented.

Chapter 6 presents the first application, the total-stress analysis soil-structure insertion problems in geomechanics. In particular, this chapter presents several analysis of increasing complexity; namely, the insertion of a rigid strip footing, the T-Bar and a rough cone penetration test. These three problems have been frequently as benchmark problems to assess the robustness and accuracy of large-strain 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 ${\textstyle u_{2}}$, 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 one-dimensional ${\textstyle \mathbf {u} -p_{w}}$ element with linear (equal-order) 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 closed-ended piles and the bearing capacity factor of tubes, a by-product 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 three-dimensional analysis are reported. In addition, in this appendix, the hyper-elastic plastic and hypo-elastic plastic constitutive models are compared in a number of representative numerical simulations.

## 1.4 Publications

This monograph, which contains some unpublished material, is based on the following publications wherein some ideas and figures have appeared previously.

Publications in journals:

• Monforte, L., J. M. Carbonell, M. Arroyo & A. Gens (2017). Performance of mixed formulations for the particle finite element method in soil mechanics problems. Computational Particle Mechanics. 4(3):269-284.
• Monforte, L., M. Arroyo, J. M. Carbonell & A. Gens (2017). Numerical simulation of undrained insertion problems in geotechnical engineering with the Particle Finite Element Method (PFEM). Computers and Geotechnics. 82:144-156.
• Monforte, L., M. Arroyo, J. M. Carbonell & A. Gens (2018). Coupled effective stress analysis of insertion problems in geotechnics with the Particle Finite Element Method. Computers and Geotechnics. 101:114-129.
• Monforte, L., M. Arroyo, J. M. Carbonell & A. Gens (2018). Hydraulic conductivity from CPTu on-the-fly: a numerical evaluation. Géotechnique Letters. 8(4) doi.org/10.1680/jgele.18.00108.

Pubications in conference proceedings:

• Monforte, L., M. Arroyo, A. Gens & J. M. Carbonell (2014). Explicit finite deformation stress integration of the elasto-plastic constitutive equations. In The 14th International Conference of the International Association for Computer Methods and Advances in Geomechanics (14IACMAG).
• Monforte, L., M. Arroyo, A. Gens & J. M. Carbonell (2015). Integration of elasto-plastic constitutive models in finite deformation: an explicit approach. In XII International Conference on Computational Plasticity. Fundamentals and applications.
• Monforte, L., J. M. Carbonell, M. Arroyo & A. Gens (2015). Numerical simulation of penetration problems in geotechnical engineering with the particle finite element method (PFEM). In IV International Conference on Particle-based Methods – Fundamentals and applications.
• Monforte, L. (2016). Numerical analysis of penetration problems in clay with the Particle Finite Element Method. In 25th European Young Geotechnical Engineers Conference.
• Gens, A., M. Arroyo, J. Butlanska, L. Monforte, J. M. Carbonell, M. Ciantia & C. O'Sullivan (2016). Simulation of the cone penetration test: discrete and continuum approaches. In 5th International Conference on Geotechnical and Geophysical Site Characterization.
• Gens, A., M. Arroyo, J. M. Carbonell, M. Ciantia & L. Monforte (2017). Simulation of penetration problems in geomechanics. In 15th International Conference on Computational Plasticity - Fundamentals and Applications.
• Monforte, L., M. Arroyo, A. Gens & J. M. Carbonell (2017). G-PFEM: a Particle Finite Element Method platform for geotechnical applications. In Alert Workshop.
• Monforte, L., M. Arroyo, A. Gens & J. M. Carbonell (2018). Numerical study of penetration problems in fine grained soils with the Particle Finite Element Method. In 4th International Symposium on Computational Geomechanics.
• Monforte, L., M. Arroyo, A. Gens & C. Parolini (2018). Permeability estimates from CPTu: a numerical study. In CPT’18 – International Symposium on Cone Penetration Testing.
• Monforte, L., M. Arroyo, A. Gens & J. M. Carbonell (2018). Three-dimensional analysis of penetration problems using G-PFEM. In 9th NUMGE Conference on Numerical Methods in Geotechnical Engineering.

# 2 Numerical model

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 single-phase and the coupled fluid-saturated porous media at large strain is presented.

## 2.1 Particle Finite Element method

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 single-phase solids [29,44] and to the contact between multiple deformable bodies in thermo-plastic 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]. Flow-like landslides have been also studied using PFEM, but considering a single-phase material description: [36] employed a rigid plastic constitutive response for the soil whereas [31] used a non-Newtonian 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 fixed-mesh 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 re-triangulation 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 re-triangulation 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, h-adaptive 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 two-dimensional 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 high-order 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.

 Figure 1: Sequence of steps to update in time a “cloud” of nodes representing a soil mass that is progressively fragmented the action of an external rigid footing using the PFEM. In the boundaries the particles are fixed.

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

1. Begin the computation of each time step with a cloud of points defining the analysis domain. For instance ${\textstyle C_{n}}$ denotes the cloud at time ${\textstyle t=t_{n}}$, see Figure 1.
2. Identify the boundaries defining the analysis domain ${\textstyle V_{n}}$. This is an essential step as some boundaries may be severely distorted during the solution and, conceptually, may include separation and re-entering of nodes. The domain boundary may be identified with the ${\textstyle \alpha {-}shape}$ method [49] or that of the previous step, using conforming meshes.
3. Discretize the continuum domains with a finite element (FE) mesh ${\textstyle M_{n}}$.
4. Solve the Lagrangian equations of motion in the domain. Compute the state variables at the next (updated) configuration for ${\textstyle t+\Delta t}$: displacements, pressure, water-pressure, strains, stresses and internal variables, etc.
5. Move the mesh nodes to a new position ${\textstyle C_{n+1}}$ where ${\textstyle n+1}$ denotes the time ${\textstyle t_{n}+\Delta t}$, in terms of the time increment size. This step is typically a consequence of the solution process of step 4.
6. Go back to step 1 and repeat the solution process for the next time step to obtain ${\textstyle C_{n+2}}$.

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 Lagrangian-Eulerian approach: the so-called Efficient Aribtrary Lagrangian-Eulerian (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, ${\textstyle h}$-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 high-order 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 Eulerian-Lagrangian 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.

### 2.1.1 Mapping between evolving meshes

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 ${\textstyle T_{i}}$ and ${\textstyle {\overline {T}}_{k}}$ be some internal variable in the new and previous mesh; a piecewice constat interpolation of these variables over the whole domain is constructed, ${\textstyle T=T_{i}w_{i}}$ and ${\textstyle {\overline {T}}={\overline {T}}_{k}{\overline {w_{k}}}}$, where the interpolation functions ${\textstyle w_{i}}$ and ${\textstyle {\overline {w}}_{k}}$ are equal to the unity in elements ${\textstyle i}$ and ${\textstyle k}$ 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:

 ${\displaystyle T_{i}=\arg \left(\min _{T_{i}}\left(\int _{\Omega }\left(T_{i}w_{i}-{\overline {T}}_{k}{\overline {w}}_{k}\right)^{2}d\Omega \right)\right)}$
(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:

 ${\displaystyle T_{i}={\frac {{\overline {T}}_{k}\int _{\Omega }{\overline {w}}_{k}w_{i}\,d\Omega }{\int _{\Omega }w_{i}\,d\Omega }}}$
(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 super-mesh 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 out-of-balance forces after remeshing are ignored and a time-step 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.

## 2.2 Continuum mechanics

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 single-phase mechanical problem, the coupled hydro-mechanical problem and then the formulation of the elasto-plastic constitutive equations. Extended details may be found in [53,54,55], among others.

### 2.2.1 Motion and deformation

#### Description of motion and time derivatives

Let ${\textstyle {\mathcal {B}}\subset \mathbb {R} ^{3}}$ define the reference configuration of a continuum body and a particle of this body is labeled by ${\textstyle \mathbf {X} \in {\mathcal {B}}}$. The position of this particle at time ${\textstyle t}$ is:

 ${\displaystyle \mathbf {x} ={\boldsymbol {\varphi }}(\mathbf {X} ,t)}$
(2.3)

where ${\textstyle {\boldsymbol {\varphi }}}$ 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, ${\textstyle \mathbf {u} (\mathbf {X} ,t)}$, can be written as:

 ${\displaystyle \mathbf {u} (\mathbf {X} ,t)={\boldsymbol {\varphi }}(\mathbf {X} ,t)-\mathbf {X} =\mathbf {x} -\mathbf {X} }$
(2.4)

The material velocity and material acceleration, denoted by ${\textstyle \mathbf {V} (\mathbf {X} ,t)}$ and ${\textstyle \mathbf {A} (\mathbf {X} ,t)}$ respectively, are defined as:

 ${\displaystyle \mathbf {V} (\mathbf {X} ,t)={\frac {\partial {\boldsymbol {\varphi }}(\mathbf {X} ,t)}{\partial t}}}$
(2.5)
 ${\displaystyle \mathbf {A} (\mathbf {X} ,t)={\frac {\partial \mathbf {V} (\mathbf {X} ,t)}{\partial t}}={\frac {\partial ^{2}{\boldsymbol {\varphi }}(\mathbf {X} ,t)}{\partial ^{2}t}}}$
(2.6)

Meanwhile, a spatial or Eulerian description is also possible:

 ${\displaystyle \mathbf {v} (\mathbf {x} ,t)=\mathbf {v} ({\boldsymbol {\varphi }}(\mathbf {X} ,t),t)}$
(2.7)
 ${\displaystyle \mathbf {a} (\mathbf {x} ,t)=\mathbf {a} ({\boldsymbol {\varphi }}(\mathbf {X} ,t),t)={\frac {\mathrm {D} \,\mathbf {v} (\mathbf {x} ,t)}{\mathrm {D} \,t}}={\frac {\partial \,\mathbf {v} (\mathbf {x} ,t)}{\partial \,t}}+\nabla \mathbf {v} (\mathbf {x} ,t)\cdot \mathbf {v} (\mathbf {x} ,t)}$
(2.8)

where the first term is known as the local derivative and the second one is the convective time derivative.

#### 2.2.1.2 Description of deformation

 Figure 2: Sketch of the reference and current configuration of a deformable body [56].

The deformation gradient, ${\textstyle \mathbf {F} }$, is the Jacobian matrix of the motion:

 ${\displaystyle \mathbf {F} (\mathbf {X} ,t)={\frac {\partial {\boldsymbol {\varphi }}(\mathbf {X} ,t)}{\partial \mathbf {X} }}}$
(2.9)

To further illustrate the deformation gradient, let us consider two different particles, ${\textstyle P}$ and ${\textstyle Q}$ in the reference configuration, whose infinitesimal relative position in the reference configuration is ${\textstyle \mathrm {d} \mathbf {X} }$, see Figure 2. Then, the relative position in the deformed configuration, ${\textstyle \mathrm {d} \mathbf {x} }$, may be obtained as:

 ${\displaystyle \mathrm {d} \mathbf {x} ={\frac {\partial \mathbf {x} }{\partial \mathbf {X} }}\cdot \mathrm {d} \mathbf {X} =\mathbf {F} \cdot \mathrm {d} \mathbf {X} }$
(2.10)

Furthermore, a differential volume in the reference configuration, ${\textstyle \mathrm {d} \Omega _{0}}$, is related to the same differential volume in the deformed configuration, ${\textstyle \mathrm {d} \Omega _{t}}$, as:

 ${\displaystyle \mathrm {d} \Omega _{t}=\det(\mathbf {F} (\mathbf {X} ,t))\,\mathrm {d} \Omega _{0}}$
(2.11)

where the Jacobian determinant is frequently referred to as the Jacobian: ${\textstyle J=\det(\mathbf {F} (\mathbf {X} ,t))}$.

Considering a differential of area in the reference configuration, ${\textstyle \mathrm {d} \mathbf {A} =\mathbf {N} \,\mathrm {d} {A}}$, where ${\textstyle \mathbf {N} }$ is the unit normal, the following expression for the transformation is given by Nanson's formula:

 ${\displaystyle {\begin{cases}\mathrm {d} \mathbf {a} =\mathbf {n} \,\mathrm {d} {a}=\det(\mathbf {F} )\,\mathbf {F} ^{-T}\cdot \mathrm {d} \mathbf {A} \\\mathrm {d} a=\det(\mathbf {F} )\,\|\mathbf {F} ^{-T}\cdot \mathbf {N} \|\,\mathrm {d} A\end{cases}}}$
(2.12)

Several deformation measures are used in non-linear continuum mechanics. Two of the most common ones are the right Cauchy-Green tensor, ${\textstyle \mathbf {C} }$, and the left Cauchy-Green tensor, ${\textstyle \mathbf {b} }$:

 ${\displaystyle \mathbf {C} (\mathbf {X} ,t)=\mathbf {F} ^{T}\cdot \mathbf {F} }$
(2.13)

 ${\displaystyle \mathbf {b} (\mathbf {x} ,t)=\mathbf {F} \cdot \mathbf {F} ^{T}}$
(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, ${\textstyle \mathbf {F} =1}$, 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:

 ${\displaystyle \mathbf {F} =\mathbf {R} \cdot \mathbf {U} =\mathbf {V} \cdot \mathbf {R} }$
(2.15)

where ${\textstyle \mathbf {R} }$ is a proper orthogonal tensor, called the rotational tensor, and ${\textstyle \mathbf {U} }$ and ${\textstyle \mathbf {V} }$ are two symmetric positive definite tensors. These tensors fulfill the following properties:

 ${\displaystyle {\begin{cases}\mathbf {R} \cdot \mathbf {R} ^{T}=1\\\mathbf {U} =\mathbf {C} ^{\frac {1}{2}}\\\mathbf {V} =\mathbf {b} ^{\frac {1}{2}}\end{cases}}}$
(2.16)

where ${\textstyle \mathbf {C} }$ and ${\textstyle \mathbf {b} }$ are, respectively, the right and left Cauchy-Green tensors. The eigenvalues of ${\textstyle \mathbf {U} }$ and ${\textstyle \mathbf {V} }$ are coincident, whereas its eigenvectors differ:

 ${\displaystyle \mathbf {R} =\sum _{i}\mathbf {n} _{i}\otimes \mathbf {N} _{i}}$ ${\displaystyle \mathbf {U} =\sum _{i}\lambda _{i}\,\,\mathbf {N} _{i}\otimes \mathbf {N} _{i}}$ ${\displaystyle \mathbf {V} =\sum _{i}\lambda _{i}\,\,\mathbf {n} _{i}\otimes \mathbf {n} _{i}}$
(2.17)

where ${\textstyle \lambda _{i}}$ are the eigenvalues of ${\textstyle \mathbf {U} }$ and ${\textstyle \mathbf {V} }$, also known as principal streches, whereas ${\textstyle \mathbf {n} _{i}}$ and ${\textstyle \mathbf {N} _{i}}$ are the eigenvectors of ${\textstyle \mathbf {U} }$ and ${\textstyle \mathbf {V} }$ respectively.

Taking into account the decomposition of the previous tensors, both Cauchy-Green tensors may be redefined as:

 ${\displaystyle \mathbf {C} =\sum _{i}\left(\lambda _{i}\right)^{2}\,\,\mathbf {N} _{i}\otimes \mathbf {N} _{i}}$ ${\displaystyle \mathbf {b} =\sum _{i}\left(\lambda _{i}\right)^{2}\,\,\mathbf {n} _{i}\otimes \mathbf {n} _{i}}$
(2.18)

#### Rate of deformation tensor

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

 ${\displaystyle {\dot {\mathbf {F} }}={\frac {\partial \mathbf {F} }{\partial t}}={\frac {\partial }{\partial t}}\left({\frac {\partial {\boldsymbol {\varphi }}(\mathbf {X} ,t)}{\partial \mathbf {X} }}\right)={\frac {\partial }{\partial \mathbf {X} }}\left({\frac {\partial {\boldsymbol {\varphi }}(\mathbf {X} ,t)}{\partial t}}\right)={\frac {\partial \mathbf {V} (\mathbf {X} ,t)}{\partial \mathbf {X} }}=\mathbf {L} }$
(2.19)

that is, the time derivative of the deformation gradient is equal to the material velocity gradient, ${\textstyle \mathbf {L} }$.

Applying the chain rule to the previous expression:

 ${\displaystyle {\frac {\partial \mathbf {V} (\mathbf {X} ,t)}{\partial \mathbf {X} }}={\frac {\partial \mathbf {v} }{\partial \mathbf {x} }}\cdot {\frac {\partial \mathbf {x} }{\partial \mathbf {X} }}}$
(2.20)

then, the spatial velocity gradient, ${\textstyle \mathbf {l} }$, may be obtained as:

 ${\displaystyle \mathbf {l} =\nabla \mathbf {v} (\mathbf {x} ,t)={\dot {\mathbf {F} }}\cdot \mathbf {F} ^{-1}}$
(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 Cauchy-Green tensor may be useful:

 ${\displaystyle {\dot {\mathbf {C} }}={\dot {\mathbf {F} }}^{T}\cdot \mathbf {F} +\mathbf {F} ^{T}\cdot {\dot {\mathbf {F} }}=\mathbf {F} ^{T}\cdot \mathbf {l} ^{T}\cdot \mathbf {F} +\mathbf {F} ^{T}\cdot \mathbf {l} \cdot \mathbf {F} =2\,\mathbf {F} ^{t}\cdot \mathbf {d} \cdot \mathbf {F} }$
(2.22)

where ${\textstyle \mathbf {d} ={\frac {1}{2}}[\mathbf {l} +\mathbf {l} ^{T}]}$ is the spatial rate of deformation tensor.

The temporal variation of the Jacobian determinant is:

 ${\displaystyle {\dot {J}}={\frac {\partial \det(\mathbf {F} )}{\partial t}}=J\,\mathrm {tr} \left(\mathbf {l} \right)}$
(2.23)

#### Description of Strain

To define the Green-Lagrange strain and the Almansi strain, let us compute the variation of length of the segment of Figure 2:

 ${\displaystyle (\mathrm {d} s)^{2}-(\mathrm {d} S)^{2}=\mathrm {d} \mathbf {x} \cdot \mathrm {d} \mathbf {x} -\mathrm {d} \mathbf {X} \cdot \mathrm {d} \mathbf {X} =\mathrm {d} \mathbf {X} \cdot \left(\mathbf {F} ^{T}\cdot \mathbf {F} -1\right)\cdot \mathrm {d} \mathbf {X} =}$ ${\displaystyle =\mathrm {d} \mathbf {x} \cdot \left(1-\mathbf {F} ^{-T}\cdot \mathbf {F} ^{-1}\right)\cdot \mathrm {d} \mathbf {x} }$
(2.24)

From this previous expressions, the Green-Lagrange strain, ${\textstyle \mathbf {E} (\mathbf {X} ,t)}$, and the Almansi strain, ${\textstyle \mathbf {e} (\mathbf {x} ,t)}$, are defined as:

 ${\displaystyle \mathbf {E} (\mathbf {X} ,t)={\frac {1}{2}}\left(\mathbf {F} ^{T}\cdot \mathbf {F} -1\right)}$ ${\displaystyle \mathbf {e} (\mathbf {x} ,t)={\frac {1}{2}}\left(1-\mathbf {F} ^{-T}\cdot \mathbf {F} ^{-1}\right)}$
(2.25)

In terms of the eigenvectors and eigenvalues of the left and right Cauchy-Green tensors, Equation 2.18, these strains may be obtained as:

 ${\displaystyle \mathbf {E} (\mathbf {X} ,t)=\sum _{i}{\frac {1}{2}}\left(\lambda _{i}^{2}-1\right)\,\,\mathbf {N} _{i}\otimes \mathbf {N} _{i}}$ ${\displaystyle \mathbf {e} (\mathbf {x} ,t)=\sum _{i}{\frac {1}{2}}\left(1-{\frac {1}{\lambda _{i}^{2}}}\right)\,\,\mathbf {n} _{i}\otimes \mathbf {n} _{i}}$
(2.26)

The unit elongation along a direction ${\textstyle \mathbf {T} }$ or ${\textstyle \mathbf {t} }$ is:

 ${\displaystyle \epsilon _{\mathbf {T} }=\epsilon _{\mathbf {t} }={\frac {\mathrm {d} s-\mathrm {d} S}{\mathrm {d} S}}}$
(2.27)

Consequently, the unit elongation may be further related to the Green-Lagrange and Almansi strain tensors through the following expression:

 ${\displaystyle \epsilon _{\mathbf {T} }=\epsilon _{\mathbf {t} }={\sqrt {1+2\,\mathbf {T} \cdot \mathbf {E} \cdot \mathbf {T} }}-1={\frac {1}{\sqrt {1-2\,\mathbf {t} \cdot \mathbf {e} \cdot \mathbf {t} }}}-1}$
(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, ${\textstyle \lambda _{\mathbf {t} }}$, is related to the unit elongation through:

 ${\displaystyle \lambda _{\mathbf {t} }=\epsilon _{\mathbf {t} }+1}$
(2.29)

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

 ${\displaystyle {\boldsymbol {\epsilon }}={\frac {\log(\mathbf {b} )}{2}}=\sum _{i}\log(\lambda _{i})\,\,\mathbf {n} _{i}\otimes \mathbf {n} _{i}}$
(2.30)
 (a) (b) Figure 3: Comparison of the Almansi, Green-Lagrange 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 one-dimensional 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${\textstyle _{2}}$, the maximum value in extension of the Almansi strain is 0.5. On the other hand, the maximum value in compression of the Green-Lagrange strain is -0.5. Potentially, the Hencky strain may take any real value.

The former observations are restricted to one-dimensional cases. It should be noted that in a general three-dimensional continuum, the Almansi and the Hencky strain share eigenvectors, that are different from those of the Green-Lagrange tensor. By definition, the Green-Lagrange tensor assumes a material description and, thus, acts to an element ${\textstyle \mathrm {d} \mathbf {X} }$ defined in the reference configuration. On the contrary, the other two strain tensors have an spatial description and acts to an element ${\textstyle \mathrm {d} \mathbf {x} }$ of the deformed configuration [56].

### 2.2.2 Stress measures

#### Stress measures

There exist a large number of stress measures. The most intuitive one is the Cauchy stress tensor, ${\textstyle {\boldsymbol {\sigma }}}$, that is defined in the current deformed configuration ${\textstyle {\boldsymbol {\varphi }}({\mathcal {B}})}$. The relation between the traction vector, ${\textstyle \mathbf {t} }$, and the surface normal vector, ${\textstyle \mathbf {n} }$, is:

 ${\displaystyle \mathbf {t} ={\boldsymbol {\sigma }}^{T}\cdot \mathbf {n} }$
(2.31)

From the local balance of angular momentum, it can be obtained that the Cauchy stress tensor is symmetric.

The Kirchhoff stress tensor, ${\textstyle \mathbf {\boldsymbol {\tau }} }$, 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:

 ${\displaystyle \mathbf {\boldsymbol {\tau }} =J\,\mathbf {\boldsymbol {\sigma }} }$
(2.32)

To define the Second Piola-Kirchhoff stress tensor, ${\textstyle \mathbf {S} }$, let us express the force acting on a differential of area:

 ${\displaystyle \mathrm {d} \mathbf {f} =\mathbf {t} \,\mathrm {d} a={\boldsymbol {\sigma }}^{T}\cdot \mathbf {n} \,\mathrm {d} a}$
(2.33)

performing a pull-back of the force and introducing the Nanson's formula:

 ${\displaystyle \mathrm {d} \mathbf {f} _{0}=\mathbf {F} ^{-1}\cdot \mathrm {d} \mathbf {f} =\mathbf {F} ^{-1}\cdot {\boldsymbol {\sigma }}^{T}\cdot \mathbf {n} \,\mathrm {d} a=\mathbf {F} ^{-1}\cdot {\boldsymbol {\sigma }}^{T}\cdot \mathbf {F} ^{-T}\cdot \mathbf {N} \,\mathrm {d} A=\mathbf {S} ^{T}\cdot \mathbf {N} \,\mathrm {d} A}$
(2.34)

which leads to the definition of the Second Piola-Kirchhoff as:

 ${\displaystyle \mathbf {S} =J\,\mathbf {F} ^{-1}\cdot {\boldsymbol {\sigma }}\cdot \mathbf {F} ^{-T}=\mathbf {F} ^{-1}\cdot {\boldsymbol {\tau }}\cdot \mathbf {F} ^{-T}}$
(2.35)

that is a symmetric stress tensor defined in the reference configuration.

#### Temporal derivatives of stress measures

The time derivative of stress tensors is of significance for the statement of incremental forms of constitutive equations. For the second Piola-Kirchhoff, that is referred to the initial configuration, the derivative with respect to time is given by the material time derivative:

 ${\displaystyle {\dot {\mathbf {S} }}={\frac {\partial \mathbf {S} (\mathbf {X} ,t)}{\partial t}}}$
(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:

 ${\displaystyle {\dot {\boldsymbol {\sigma }}}={\frac {\partial {\boldsymbol {\sigma }}(\mathbf {x} ,t)}{\partial t}}+{\frac {\partial {\boldsymbol {\sigma }}(\mathbf {x} ,t)}{\partial \mathbf {x} }}\cdot \mathbf {v} (\mathbf {x} ,t)}$
(2.37)

### 2.2.3 Objective transformation

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 ${\textstyle {\boldsymbol {\varphi }}(\mathbf {X} ,t)}$, a superposed rigid motion of the current place is the map ${\textstyle {\boldsymbol {\varphi }}^{+}(\mathbf {x} ,t)}$ defined as:

 ${\displaystyle {\boldsymbol {\varphi }}^{+}=\mathbf {c} (t)+\mathbf {Q} (t)\cdot \mathbf {x} }$
(2.38)

where ${\textstyle \mathbf {c} (t)}$ is a function of time and ${\textstyle \mathbf {Q} (t)}$ is a proper orthogonal transformation depending only on time. The superposed motion is called rigid because, for any given two points ${\textstyle \mathbf {x} _{1}}$, ${\textstyle \mathbf {x} _{2}\in \varphi ({\mathcal {B}},t)}$, the following property holds

 ${\displaystyle \mathbf {x} _{1}^{+}-\mathbf {x} _{2}^{+}=\mathbf {Q} \left(\mathbf {x} _{1}-\mathbf {x} _{2}\right)\implies \|\mathbf {x} _{1}^{+}-\mathbf {x} _{2}^{+}\|=\|\mathbf {x} _{1}-\mathbf {x} _{2}\|}$
(2.39)

in other words, the mapping preserves the distances.

 ${\displaystyle \mathbf {F} ^{+}=\mathbf {Q} (t)\cdot \mathbf {F} }$
(2.40)

Then, the left Cauhy-Green tensor, Equation 2.14, reads:

 ${\displaystyle \mathbf {b} ^{+}=\mathbf {F} ^{+}\cdot \mathbf {F} ^{+T}=\mathbf {Q} \cdot \mathbf {b} \cdot \mathbf {Q} ^{T}}$
(2.41)

which transforms objectively.

On the contrary, material objects remain unaltered to the superposed rigid body motions; for instance, the right Cauchy-Green tensor, Equation 2.13:

 ${\displaystyle \mathbf {C} ^{+}=\mathbf {F} ^{+T}\cdot \mathbf {F} ^{+}=\mathbf {F} ^{T}\cdot \mathbf {Q} ^{T}\cdot \mathbf {Q} \cdot \mathbf {F} =\mathbf {F} ^{T}\cdot \mathbf {F} =\mathbf {C} }$
(2.42)

The spatial velocity gradient, Equation 2.21 is given by:

 ${\displaystyle \mathbf {l} ^{+}=\mathbf {Q} (t)\cdot \mathbf {l} \cdot \mathbf {Q} ^{T}(t)+{\dot {\mathbf {Q} }}(t)\cdot \mathbf {Q} ^{T}(t)}$
(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:

 ${\displaystyle \mathbf {d} ^{+}=\mathbf {Q} (t)\cdot \mathbf {d} \cdot \mathbf {Q} ^{T}(t)}$
(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, ${\textstyle \mathbf {t} ^{+}=\mathbf {Q} \cdot \mathbf {t} }$ whereas unit outward normal is objective, ${\textstyle \mathbf {n} ^{+}=\mathbf {Q} \cdot \mathbf {n} }$. Then:

 ${\displaystyle {\begin{cases}\mathbf {t} ={\boldsymbol {\sigma }}\cdot \mathbf {n} \\\mathbf {t} ^{+}={\boldsymbol {\sigma }}^{+}\cdot \mathbf {n} ^{+}=\mathbf {Q} \cdot \mathbf {t} \end{cases}}}$
(2.45)

it follows from the previous expression that ${\textstyle {\boldsymbol {\sigma }}^{+}\cdot \mathbf {Q} \cdot \mathbf {n} =\mathbf {Q} \cdot {\boldsymbol {\sigma }}\cdot \mathbf {n} }$. Since this equation holds for any normal, it follows that:

 ${\displaystyle {\boldsymbol {\sigma }}^{+}=\mathbf {Q} (t)\cdot {\boldsymbol {\sigma }}\cdot \mathbf {Q} ^{T}(t)}$
(2.46)

that is, the Cauchy stress tensor transforms objectively.

However, the temporal derivative of the Cauchy stress tensor is not objective:

 ${\displaystyle {\dot {\boldsymbol {\sigma }}}^{+}=\mathbf {Q} (t)\cdot {\dot {\boldsymbol {\sigma }}}\cdot \mathbf {Q} ^{T}(t)+\left({\dot {\mathbf {Q} }}(t)\cdot \mathbf {Q} ^{T}(t)\right)\cdot {\boldsymbol {\sigma }}^{+}-{\boldsymbol {\sigma }}^{+}\cdot \left({\dot {\mathbf {Q} }}(t)\cdot \mathbf {Q} ^{T}(t)\right)}$
(2.47)

that is clearly non-objective since it is not equal to ${\textstyle \mathbf {Q} (t)\cdot {\dot {\boldsymbol {\sigma }}}\cdot \mathbf {Q} ^{T}(t)}$.

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:

 ${\displaystyle {\mathcal {L}}_{v}{\boldsymbol {\tau }}=\mathbf {F} \cdot {\frac {\partial }{\partial t}}\left(\mathbf {F} ^{-1}\cdot {\boldsymbol {\tau }}\cdot \mathbf {F} ^{-T}\right)\cdot \mathbf {F} ^{T}=\mathbf {F} \cdot {\frac {\partial \mathbf {S} }{\partial t}}\cdot \mathbf {F} ^{T}}$
(2.48)

using the expression of the derivative of the inverse and arranging some terms:

 ${\displaystyle {\mathcal {L}}_{v}{\boldsymbol {\tau }}={\dot {\boldsymbol {\tau }}}-\mathbf {l} \cdot {\boldsymbol {\tau }}-{\boldsymbol {\tau }}\cdot \mathbf {l} }$
(2.49)

## 2.3 Single-phase mechanical media

In this section, the formulation for the single-phase 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 time-stepping algorithm. The non-linear system of equations is solved with a Newton-Raphson method; thus, the required linearization is derived.

### 2.3.1 Strong form of the balance equations

A quasi-static linear momentum formulation in an updated-Lagrangian form (i.e. expressing all quantities and their derivatives in the deformed configuration), may be written as:

 ${\displaystyle {\begin{cases}\nabla \cdot {\boldsymbol {\sigma }}+\mathbf {b} ={\boldsymbol {0}}&{\hbox{ in }}\Omega _{t}\times (0,T)\\{\boldsymbol {u}}({\boldsymbol {X}},t=0)={\boldsymbol {u}}_{0}&{\hbox{ in }}\Omega _{0}\\{\boldsymbol {u}}({\boldsymbol {X}},t)={\overline {\boldsymbol {u}}}&{\hbox{ in }}\Gamma _{u}\times (0,T)\\{\boldsymbol {n}}\cdot {\boldsymbol {\sigma }}={\overline {\boldsymbol {t}}}&{\hbox{ in }}\Gamma _{\overline {t}}\times (0,T)\\\end{cases}}}$
(2.50)

where ${\textstyle {\boldsymbol {\sigma }}={\hat {\boldsymbol {\sigma }}}({\boldsymbol {F}},V)}$ is the Cauchy stress tensor, ${\textstyle {\hat {\boldsymbol {\sigma }}}}$ stands for the appropriate constitutive equation for path dependent materials (large strains elasto-plastic constitutive equations based on the multiplicative split [57] are used here, see Chapter 3), ${\textstyle {\boldsymbol {F}}}$ is the total deformation gradient and ${\textstyle V}$ represents the set of internal variables of the model. ${\textstyle \mathbf {u} _{0}}$ stands for the initial displacement, ${\textstyle \mathbf {b} =\rho \,\mathbf {g} }$ is the external body force vector, ${\textstyle \mathbf {g} }$ is the gravity vector and ${\textstyle \partial \Omega _{t}=\Gamma _{u}\cup \Gamma _{\overline {t}}}$ (${\textstyle \Gamma _{u}\cap \Gamma _{\overline {t}}=\emptyset }$) defines the boundary of the domain where displacements ${\textstyle {\overline {\boldsymbol {u}}}}$ and tractions ${\textstyle {\overline {\boldsymbol {t}}}}$ 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:

 ${\displaystyle \rho (\mathbf {X} ,t=0)=\rho (\mathbf {X} ,t)\,\det \left(\mathbf {F} (\mathbf {X} ,t)\right)}$
(2.51)

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

 ${\displaystyle {\boldsymbol {\sigma }}={\boldsymbol {\sigma }}^{T}}$
(2.52)

### 2.3.2 Weak form

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, ${\textstyle \mathbf {w} }$, and integrating it in the domain:

 ${\displaystyle \int _{\Omega _{t}}\mathbf {w} \cdot \left(\nabla \cdot {\boldsymbol {\sigma }}+\mathbf {b} \right)\,\,\mathrm {d} \Omega _{t}=0}$
(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:

 ${\displaystyle \int _{\Omega _{t}}{\frac {\partial w_{j}}{\partial x_{i}}}\sigma _{ij}\,\,\mathrm {d} \Omega _{t}=\int _{\Omega _{t}}w_{j}b_{j}\,\,\mathrm {d} \Omega _{t}+\int _{\Gamma _{\overline {t}}}w_{j}\cdot {\overline {t}}_{j}\,\,\mathrm {d} \gamma _{t}}$
(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 non-linearities do appear, however, they are hidden.

### 2.3.3 Finite element discretization

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

 ${\displaystyle {\begin{cases}\mathbf {u} \approx \mathbf {u} ^{h}=\mathbf {N} _{u}\cdot {\tilde {\mathbf {u} }}\\\mathbf {w} \approx \mathbf {w} ^{h}=\mathbf {N} _{u}\cdot {\tilde {\mathbf {w} }}\end{cases}}}$
(2.55)

where ${\textstyle \mathbf {u} ^{h}}$ is the finite element approximation of displacement whereas ${\textstyle {\tilde {\mathbf {u} }}}$ are the nodal values. ${\textstyle \mathbf {N} _{u}=[N_{1}1,N_{2}1,...,N_{n}1]}$ are the shape functions whereas ${\textstyle n}$ is the number of nodes. The same procedure is used to discretize the test functions, ${\textstyle \mathbf {w} }$.

As such, the problem is then: Find the function ${\textstyle \mathbf {u} ^{h}(t)}$ such that:

 ${\displaystyle \int _{\Omega _{t}}{\frac {\partial w_{j}^{h}}{\partial x_{i}}}\sigma _{ij}\,\,\mathrm {d} \Omega _{t}-\int _{\Omega _{t}}w_{j}^{h}b_{j}\,\,\mathrm {d} \Omega _{t}-\int _{\Gamma _{\overline {t}}}w_{j}^{h}\cdot {\overline {t}}_{j}\,\,\mathrm {d} \gamma _{t}=0}$
(2.56)

Since ${\textstyle \mathbf {w} ^{h}}$ are arbitrary, from the previous expression ${\textstyle n_{dim}\,n}$ independent equations may be obtained, where ${\textstyle n_{dim}}$ is the number of spatial dimensions of the problem. As such, the matrix form of the Galerkin weak form is:

 ${\displaystyle \mathbf {P} \left({\boldsymbol {\sigma }}\right)=\mathbf {f} ^{ext}}$
(2.57)

where the internal and external forces read:

 ${\displaystyle \mathbf {P} ({\boldsymbol {\sigma }})=\int _{\Omega _{t}}\mathbf {B} ^{T}\cdot {\underline {\boldsymbol {\sigma }}}\;\mathrm {d} \Omega _{t}}$
(2.58)
 ${\displaystyle \mathbf {f} ^{ext}=\int _{\Omega _{t}}\mathbf {N} _{u}^{T}\cdot \mathbf {b} \;\mathrm {d} \Omega _{t}+\int _{\Gamma _{\overline {t}}}\mathbf {N} _{u}^{T}\cdot {\overline {\mathbf {t} }}\;\mathrm {d} \gamma _{t}}$
(2.59)

where ${\textstyle \mathbf {B} }$ has the same formal structure than the small deformation strain-displacement matrix [58] and ${\textstyle {\underline {\boldsymbol {\sigma }}}}$ corresponds to the Voigt notation of tensor ${\textstyle {\boldsymbol {\sigma }}}$.

### 2.3.4 Temporal discretization

As already mentioned, a quasi-static 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 ${\textstyle t_{n+1}}$, as: find ${\textstyle ^{n+1}\mathbf {u} ^{h}}$ such that:

 ${\displaystyle \mathbf {P} \left(^{n+1}{\boldsymbol {\sigma }}\right)=^{n+1}\mathbf {f} ^{ext}}$
(2.60)

### 2.3.5 Linearization

Equation 2.60 defines a non-linear system of equations, whose residual is given by:

 ${\displaystyle \mathbf {R} (^{n+1}\mathbf {u} )=\mathbf {P} \left(^{n+1}{\boldsymbol {\sigma }}\right)-^{n+1}\mathbf {f} ^{ext}}$
(2.61)

where the superscript ${\textstyle {\tilde {\mathbf {u} }}}$ is obviated for simplicity. Then, to evaluate numerically the solution, an iterative non-linear solver is required. In this work, the Newton-Raphson method is used. Writing an iterative correction as:

 ${\displaystyle _{i+1}^{n+1}\mathbf {u} =_{i+1}^{n+1}\mathbf {u} +_{i+1}\delta \mathbf {u} }$
(2.62)

where the subscript ${\textstyle {i+1}}$ stands for the number of iteration. Using a Taylor's approximation of the residual, Equation 2.61,

 ${\displaystyle \mathbf {R} (_{i+1}^{n+1}\mathbf {u} )\approx \mathbf {R} (_{i}^{n+1}\mathbf {u} )+\mathbf {K} (_{i}^{n+1}\mathbf {u} )\cdot _{i+1}\delta \mathbf {u} }$
(2.63)

where ${\textstyle \mathbf {K} (_{i}^{n+1}\mathbf {u} )}$ is the so-called tangent or Jacobian matrix, that is defined as:

 ${\displaystyle \mathbf {K} (\mathbf {u} )={\frac {\partial \mathbf {R} (\mathbf {u} )}{\partial \mathbf {u} }}}$
(2.64)

Assuming that ${\textstyle \mathbf {R} (_{i+1}^{n+1}\mathbf {u} )=\mathbf {0} }$ is null, the iterative correction, ${\textstyle _{i+1}\delta \mathbf {u} }$, is given by the following system of linear equations:

 ${\displaystyle \mathbf {K} (_{i}^{n+1}\mathbf {u} )\cdot _{i+1}\delta \mathbf {u} =-\mathbf {R} (_{i}^{n+1}\mathbf {u} )}$
(2.65)

It is worth noting that the term ${\textstyle \mathbf {K} (\mathbf {u} )\cdot \delta \mathbf {u} }$ may be also obtained by using the directional or Gateaux derivative along the direction ${\textstyle \delta \mathbf {u} }$ [54]. That is:

 ${\displaystyle \delta _{\mathbf {u} }\mathbf {R} (\mathbf {u} )=\left.{\frac {\mathrm {d} \mathbf {R} (\mathbf {u} +\epsilon \delta \mathbf {u} )}{\mathrm {d} \epsilon }}\right|_{\epsilon =0}={\frac {\partial \mathbf {R} (\mathbf {u} )}{\partial \mathbf {u} }}\cdot \delta \mathbf {u} =\mathbf {K} (\mathbf {u} )\cdot \delta \mathbf {u} }$
(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 so-called 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 pull-back of the internal forces term is performed; this way, all the quantities are referred to the reference configuration:

 ${\displaystyle \int _{\Omega _{t}}{\frac {\partial w_{j}}{\partial x_{i}}}\sigma _{ij}\,\,\mathrm {d} \Omega _{t}=\int _{\Omega _{0}}{\frac {\partial w_{j}}{\partial X_{K}}}F_{Ki}^{-1}\sigma _{ij}\,J\,\mathrm {d} \Omega _{0}=\int _{\Omega _{0}}{\frac {\partial w_{j}}{\partial X_{K}}}S_{KL}F_{jL}^{}\,\,\mathrm {d} \Omega _{0}}$
(2.67)

Then:

 ${\displaystyle \delta _{\mathbf {u} }\left(\int _{\Omega _{0}}{\frac {\partial w_{j}}{\partial X_{K}}}S_{KL}F_{jL}^{}\,\mathrm {d} \Omega _{0}\right)=\int _{\Omega _{0}}{\frac {\partial w_{j}}{\partial X_{K}}}\left(\delta _{\mathbf {u} }({S}_{KL})F_{jL}^{}+S_{KL}\delta _{\mathbf {u} }({F}_{jL})\right)\,\mathrm {d} \Omega _{0}}$
(2.68)

Then, a push-forward of this expression yields:

 ${\displaystyle \int _{\Omega _{t}}{\frac {\partial w_{j}}{\partial x_{i}}}F_{iK}\left(\delta _{\mathbf {u} }\left({S}_{KL}\right)F_{jL}^{}+S_{KL}\delta _{\mathbf {u} }\left({F}_{jL}\right)\right)\,{\frac {1}{J}}\,\mathrm {d} \Omega _{t}}$
(2.69)

where the directional derivative of the deformation gradient, ${\textstyle \mathbf {F} }$, is expressed as:

 ${\displaystyle \delta _{\mathbf {u} }\mathbf {F} =\left.{\frac {\mathrm {d} }{\mathrm {d} \epsilon }}\left({\frac {\partial \left(\mathbf {x} +\epsilon \delta \mathbf {u} \right)}{\partial \mathbf {X} }}\right)\right|_{\epsilon =0}={\frac {\partial \delta \mathbf {u} }{\partial \mathbf {X} }}={\frac {\partial \delta \mathbf {u} }{\partial \mathbf {x} }}\cdot {\frac {\partial \mathbf {x} }{\partial \mathbf {X} }}=\nabla \delta \mathbf {u} \cdot \mathbf {F} }$
(2.70)

whereas the term related to the second Piola-Kirchhoff may be related to the Lie derivative of the Kirchhoff stress and, also, to the spatial constitutive tensor:

 ${\displaystyle \mathbf {F} \cdot \delta _{\mathbf {u} }\mathbf {S} \cdot \mathbf {F} ^{T}={\mathcal {L}}_{\delta \mathbf {u} }{\boldsymbol {\tau }}={\boldsymbol {\mathbb {D} }}^{\ast }:\delta \mathbf {u} }$
(2.71)

where ${\textstyle {\boldsymbol {\mathbb {D} }}^{\ast }}$ is the spatial constitutive tensor.

So finally,

 ${\displaystyle \delta _{\mathbf {u} }\left(\int _{\Omega _{t}}{\frac {\partial w_{j}}{\partial x_{i}}}\sigma _{ij}\,\,\mathrm {d} \Omega _{t}\right)=\int _{\Omega _{t}}{\frac {\partial w_{j}}{\partial x_{i}}}\left(\mathbb {D} _{ijkl}+\sigma _{il}\delta _{jk}\right){\frac {\partial \delta u_{k}}{\partial x_{l}}}\,\mathrm {d} \Omega _{t}}$
(2.72)

where ${\textstyle {\boldsymbol {\mathbb {D} }}={\boldsymbol {\mathbb {D} }}^{\ast }/J}$.

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 non-linear term.

 ${\displaystyle \delta _{\mathbf {u} }\left(\int _{\Omega _{t}}w_{i}\,\rho g_{i}\;\mathrm {d} \Omega _{t}+\int _{\Gamma _{\overline {t}}}w_{i}\cdot {\overline {t}}_{i}\;\mathrm {d} \gamma _{t}\right)=\delta _{\mathbf {u} }\left(\int _{\Omega _{0}}w_{i}\,\rho _{0}g_{i}\;\mathrm {d} \Omega _{0}\right)+\delta _{\mathbf {u} }\left(\int _{\Gamma _{\overline {t}}}w_{i}\cdot {\overline {t}}_{i}\;\mathrm {d} \gamma _{t}\right)}$
(2.73)

The pull-back of the term related to the gravitational loads reveals that there is no dependency on the displacements; that is

 ${\displaystyle \delta _{\mathbf {u} }\left(\int _{\Omega _{0}}w_{i}\,\rho _{0}g_{i}\;\mathrm {d} \Omega _{0}\right)=0}$
(2.74)
 Figure 4: Convective coordinates on a parametrized surface. [59].

On the other side, a pull-back of the term related to the imposed tractions, making use of a parametrized description of the boundary, see Figure 4, reads:

 ${\displaystyle \delta _{\mathbf {u} }\left(\int _{\Gamma _{\overline {t}}}w_{i}\,{\overline {t}}_{i}\,\,\mathrm {d} \gamma _{t}\right)=\delta _{\mathbf {u} }\left(\int _{\square }w_{i}{\overline {t}}_{i}\|\mathbf {g} _{1}\times \mathbf {g} _{2}\|\mathrm {d} \theta _{1}\mathrm {d} \theta _{2}\right)=\int _{\square }w_{i}\,\,{\overline {t}}_{i}\,\,\delta _{\mathbf {u} }\left(\|\mathbf {g} _{1}\times \mathbf {g} _{2}\|\right)\,\,\mathrm {d} \theta _{1}\mathrm {d} \theta _{2}}$
(2.75)

where ${\textstyle \mathbf {g} _{\alpha }={\dfrac {\partial {\boldsymbol {\varphi }}}{\partial \theta _{\alpha }}}}$. By using the chain rule, this term reads:

 ${\displaystyle \int _{\square }w_{i}\,\,{\overline {t}}_{i}\,\,{\frac {\mathbf {g} _{1}\times \mathbf {g} _{2}}{\|\mathbf {g} _{1}\times \mathbf {g} _{2}\|}}\cdot \delta _{\mathbf {u} }{\big (}\mathbf {g} _{1}\times \mathbf {g} _{2}{\big )}\,\,\mathrm {d} \theta _{1}\mathrm {d} \theta _{2}}$
(2.76)

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

 ${\displaystyle \delta _{\mathbf {u} }\left(\mathbf {g} _{1}\times \mathbf {g} _{2}\right)=\delta _{\mathbf {u} }(\mathbf {g} _{1})\times \mathbf {g} _{2}+\mathbf {g} _{1}\times \delta _{\mathbf {u} }(\mathbf {g} _{2})=\delta \mathbf {u} _{,1}\times \mathbf {g} _{2}+\mathbf {g} _{1}\times \delta \mathbf {u} _{,2}}$
(2.77)

where ${\textstyle \delta \mathbf {u} _{,\alpha }={\frac {\partial \delta \mathbf {u} }{\partial \theta _{\alpha }}}}$.

Finally, introducing the expression 2.77 to 2.76 and performing a push-forward:

 ${\displaystyle \delta _{\mathbf {u} }\left(\int _{\Gamma _{\overline {t}}}w_{i}\,{\overline {t}}_{i}\,\,\mathrm {d} \gamma _{t}\right)=\int _{\Gamma _{\overline {t}}}w_{i}{\overline {t}}_{i}{\Big (}\delta \mathbf {u} _{,1}\cdot \left(\mathbf {g} _{2}\times \mathbf {n} \right)-\delta \mathbf {u} _{,2}\cdot \left(\mathbf {g} _{1}\times \mathbf {n} \right){\Big )}\mathrm {d} \gamma _{t}}$
(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 ${\textstyle {\overline {t}}\,\mathbf {N} }$ is imposed in an area ${\textstyle A\mathbf {N} }$ at the initial state; then, the imposed force, is equal to ${\textstyle A\,{\overline {t}}\,\mathbf {N} }$. Assuming that the body deforms but the normal to the surface where Neumann condition remains unaltered (${\textstyle \mathbf {n} =\mathbf {N} }$), the area of the surface is ${\textstyle a\,\mathbf {N} }$ and the total applied load at the Neumann boundary is ${\textstyle a\,{\overline {t}}\,\mathbf {N} }$. 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:

 ${\displaystyle \mathbf {K} =\int _{\Omega _{t}}\mathbf {B} ^{T}\cdot {\underline {D}}\cdot \mathbf {B} \,\,\mathrm {d} \Omega _{t}+\int _{\Omega _{t}}\mathbf {B} _{NL}^{T}\cdot {\underline {\underline {\boldsymbol {\sigma }}}}\cdot \mathbf {B} _{NL}\,\,\mathrm {d} \Omega _{t}-}$ ${\displaystyle -\int _{\Gamma _{\overline {t}}}\left(\mathbf {N} _{u}\cdot \mathbf {t} \right)\cdot {\Big (}\mathbf {N} _{u,1}\cdot \left(\mathbf {g} _{2}\times \mathbf {n} \right)-\mathbf {N} _{u,2}\cdot \left(\mathbf {g} _{1}\times \mathbf {n} \right){\Big )}\mathrm {d} \gamma _{t}}$
(2.79)

where the terms involed in the computation of the nonlinear stiffness matrix, ${\textstyle \mathbf {B} _{NL}}$ and ${\textstyle {\underline {\underline {\boldsymbol {\sigma }}}}}$, may be described, for two dimensional analysis, as [60,54]:

 ${\displaystyle \mathbf {B} _{NL}=\left({\begin{array}{ccccccc}N_{1,1}&0&N_{2,1}&0&\cdots &N_{n,1}&0\\N_{1,2}&0&N_{2,2}&0&\cdots &N_{n,2}&0\\0&N_{1,1}&0&N_{2,1}&\cdots &0&N_{n,1}\\0&N_{1,2}&0&N_{2,2}&\cdots &0&N_{n,2}\\N_{1}/r&0&N_{2}/r&0&\cdots &N_{n}/r&0\\\end{array}}\right)}$
(2.80)

 ${\displaystyle {\underline {\underline {\boldsymbol {\sigma }}}}=\left({\begin{array}{cccc|c}\sigma _{11}&\sigma _{12}&0&0&0\\\sigma _{21}&\sigma _{22}&0&0&0\\0&0&\sigma _{11}&\sigma _{12}&0\\0&0&\sigma _{21}&\sigma _{22}&0\\0&0&0&0&\sigma _{33}\\\end{array}}\right)}$
(2.81)

where ${\textstyle N_{1,1}}$ is the derivative of the local shape function 1 with respect to coordinate 1 and ${\textstyle r}$ is the deformed radial coordinate. It should be noted that the expressions provided in Equations 2.80 and 2.81 are particularized for a two-dimensional axisymmetric case. For a plane-strain 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 three-dimensional analysis may be found in [60].

## 2.4 Two-phase hydro-mechanical problem

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 so-called ${\textstyle \mathbf {u} -\mathbf {w} -p_{w}}$ (solid displacement-Darcy's velocity-water pressure) and the simplified ${\textstyle \mathbf {u} -p_{w}}$ 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 poro-mechanical 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 fluid-saturated 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 time-integration algorithm, the consistent linearization required for the non-linear 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 hydro-mechanical problems.

### 2.4.1 Mass balance equation

In this subsection the equation for the mass balance of the mixture is derived.

Let ${\textstyle \rho _{s}}$ and ${\textstyle \rho _{w}}$ be the density of the solid particles and the fluid, respectively, and ${\textstyle \varphi }$ the porosity in the current configuration. Then, the mass of solid and liquid, ${\textstyle m_{s}}$ and ${\textstyle m_{w}}$, in an arbitrary deformed domain, ${\textstyle \Omega _{i}}$, may be defined by the following integrals:

 ${\displaystyle m_{s}=\int _{\Omega _{i}}\rho _{s}\,(1-\varphi )\,\mathrm {d} \Omega }$ ${\displaystyle m_{w}=\int _{\Omega _{i}}\rho _{w}\,\varphi \,\mathrm {d} \Omega }$
(2.82)

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

 ${\displaystyle {\frac {\mathrm {d} m_{s}}{\mathrm {d} t}}=0\implies \int _{\Omega _{i}}\left({\frac {\partial \left(\rho _{s}\left(1-\varphi \right)\right)}{\partial t}}+{\frac {\partial \left(\rho _{s}\left(1-\varphi \right)v_{i}\right)}{\partial x_{i}}}\right)d\Omega =0}$ ${\displaystyle {\frac {{\hbox{D}}m_{w}}{{\hbox{D}}t}}=0\implies \int _{\Omega _{i}}\left({\frac {\partial \left(\rho _{w}\,\varphi \right)}{\partial t}}+{\frac {\partial \left(\rho _{w}\,\varphi \,{\overline {v_{i}}}\right)}{\partial x_{i}}}\right)d\Omega =0}$
(2.83)

where ${\textstyle {\frac {\mathrm {d} }{\mathrm {d} t}}}$ and ${\textstyle {\frac {\hbox{D}}{{\hbox{D}}t}}}$ are the material derivatives with respect to the solid and fluid phase and ${\textstyle \mathbf {v} }$ and ${\textstyle {\overline {\mathbf {v} }}}$ 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, ${\textstyle {\frac {\mathrm {d} }{\mathrm {d} t}}={\frac {\partial }{\partial t}}+\mathbf {v} \cdot \nabla }$, the last expression may be recasted as:

 ${\displaystyle {\frac {\mathrm {d} {\big (}\rho _{s}\left(1-\varphi \right){\big )}}{\mathrm {d} t}}+\rho _{s}\left(1-\varphi \right)\nabla \cdot \mathbf {v} =0}$ ${\displaystyle {\frac {\mathrm {d} {\big (}\rho _{w}\varphi {\big )}}{\mathrm {d} t}}+\rho _{w}\varphi \,\nabla \cdot \mathbf {v} +\nabla \cdot {\big (}\rho _{w}\,\mathbf {v} ^{d}{\big )}=0}$
(2.84)

where the Darcy velocity has been introduced: ${\textstyle \mathbf {v} ^{d}=\varphi \left({\overline {\mathbf {v} }}-\mathbf {v} \right)}$.

Expanding the derivatives and dividing by the density, the following expressions may be obtained:

 ${\displaystyle -{\frac {\mathrm {d} \varphi }{\mathrm {d} t}}+{\frac {\left(1-\varphi \right)}{\rho _{s}}}{\frac {\mathrm {d} \rho _{s}}{\mathrm {d} t}}+\left(1-\varphi \right)\nabla \cdot \mathbf {v} =0}$ ${\displaystyle {\frac {\mathrm {d} \varphi }{\mathrm {d} t}}+{\frac {\varphi }{\rho _{w}}}{\frac {\mathrm {d} \rho _{w}}{\mathrm {d} t}}+\varphi \nabla \cdot \mathbf {v} +{\frac {1}{\rho _{w}}}\nabla \cdot (\rho _{w}\mathbf {v} ^{d})=0}$
(2.85)

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

 ${\displaystyle {\frac {\left(1-\varphi \right)}{\rho _{s}}}{\frac {\mathrm {d} \rho _{s}}{\mathrm {d} t}}+{\frac {\varphi }{\rho _{w}}}{\frac {\mathrm {d} \rho _{w}}{\mathrm {d} t}}+\nabla \cdot \mathbf {v} +{\frac {1}{\rho _{w}}}\nabla \cdot (\rho _{w}\mathbf {v} ^{d})=0}$
(2.86)

#### Compressibility of the constituents

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

 ${\displaystyle {\frac {\mathrm {d} \rho _{s}}{\mathrm {d} t}}=0}$
(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:

 ${\displaystyle \epsilon _{v}^{f}=-\ln \left({\frac {\rho _{w}}{\rho _{0}^{w}}}\right)=\ln \left({\frac {\rho _{0}^{w}}{\rho _{w}}}\right)}$
(2.88)

where ${\textstyle \rho _{0}^{w}}$ is a reference fluid density.

Based on thermo-mechanical 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:

 ${\displaystyle p_{w}=K_{w}\epsilon _{v}=-K_{w}\ln \left({\frac {\rho _{w}}{\rho _{0}^{w}}}\right)\implies \rho _{w}=\rho _{0}^{w}\exp \left(-{\frac {p_{w}}{K_{w}}}\right)}$
(2.89)

where ${\textstyle K_{w}}$ is the water bulk modulus.

As a consequence:

 ${\displaystyle {\frac {\mathrm {d} \rho _{w}}{\mathrm {d} t}}=-{\frac {\rho _{w}}{K_{w}}}{\frac {\mathrm {d} p_{w}}{\mathrm {d} t}}}$ ${\displaystyle \nabla \rho _{w}=-{\frac {\rho _{w}}{K_{w}}}\nabla p_{w}}$
(2.90)

#### Darcy's law

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:

 ${\displaystyle \mathbf {v} ^{d}={\boldsymbol {K}}\cdot \left({\frac {\nabla \pi _{w}}{Jg\rho _{w}}}+{\frac {\mathbf {g} }{g}}\right)}$
(2.91)

where ${\textstyle \pi _{w}=J\,p_{w}}$ is the Kirchhoff water pressure, ${\textstyle {\boldsymbol {K}}}$ is the permeability tensor, ${\textstyle \mathbf {g} }$ is the gravity and ${\textstyle g=\|\mathbf {g} \|}$. This is done since, in the referred work, the Kirchhoff water pressure is a nodal variable and because, accidentally, it is assumed that ${\textstyle \nabla J={\boldsymbol {0}}}$, which, in the general case, is not true.

On the other hand, in this work the expression proposed by [64] is used:

 ${\displaystyle \mathbf {v} ^{d}={\boldsymbol {k}}\cdot \left(\nabla p_{w}+\rho _{w}\mathbf {g} \right)}$
(2.92)

where ${\textstyle {\boldsymbol {k}}={\frac {{\boldsymbol {k}}^{\prime }}{g\rho _{w}}}}$ and ${\textstyle {\boldsymbol {k}}^{\prime }}$ 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]:

 ${\displaystyle \mathbf {V} ^{d}=J\,\mathbf {F} ^{-1}\cdot \mathbf {v} ^{d}=J\,\mathbf {F} ^{-1}\cdot {\boldsymbol {k}}\cdot \mathbf {F} ^{-T}\cdot \left(\nabla _{\mathbf {X} }p_{w}+\rho _{w}\,\mathbf {F} ^{T}\cdot \mathbf {g} \right)}$
(2.93)

where ${\textstyle \nabla _{\mathbf {X} }p_{w}}$ is the material gradient of the water pressure. As a consequence of the previous expression, we can relate the permeability in the reference configuration, ${\textstyle {\boldsymbol {k}}_{0}}$, and in the current configuration as: ${\textstyle {\boldsymbol {k}}={\frac {1}{J}}\,\mathbf {F} \cdot {\boldsymbol {k}}_{0}\cdot \mathbf {F} ^{T}}$

#### Constitutive models for the permeability tensor

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, ${\textstyle {\boldsymbol {k}}_{0}}$, 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, ${\textstyle {\boldsymbol {k}}}$, 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 Kozeny-Carman 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:

 ${\displaystyle {\begin{cases}{\boldsymbol {k}}=k(e)\,1\\[6pt]k(e)=C\,{\frac {e^{3}}{1+e}}&{\hbox{ where }}C=k_{r}{\frac {1+e_{0}}{e_{0}^{3}}}\end{cases}}}$
(2.94)

where ${\textstyle k_{r}}$ and ${\textstyle e_{0}}$ are the permeability and void ratio at the initial state and ${\textstyle e}$ is the void ratio, that is defined according to:

 ${\displaystyle e={\frac {\varphi }{1-\varphi }}}$
(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 Kozeny-Carman law is used, the terms that appear due to the dependency of the permeability on the volumetric strains are not considered.

#### Final form of the mass balance equation

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:

 ${\displaystyle {\frac {-\varphi }{K_{w}}}{\frac {\mathrm {d} p_{w}}{\mathrm {d} t}}+\nabla \cdot \mathbf {v} +\nabla \cdot \mathbf {v} ^{d}+{\frac {\mathbf {v} ^{d}}{\rho _{w}}}\cdot \nabla \rho _{w}=0}$
(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:

 ${\displaystyle {\frac {-\varphi }{K_{w}}}{\frac {\mathrm {d} p_{w}}{\mathrm {d} t}}+\nabla \cdot \mathbf {v} +\nabla \cdot {\big (}{\boldsymbol {k}}\cdot \left(\nabla p_{w}+\rho _{w}\mathbf {g} \right){\big )}-{\frac {\mathbf {v} ^{d}}{K_{w}}}\cdot \nabla p_{w}=0}$
(2.97)

It is assumed that the term ${\textstyle {\dfrac {\mathbf {v} ^{d}}{K_{w}}}\cdot \nabla p_{w}\approx 0}$, 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 ${\textstyle {\dfrac {\varphi }{K_{w}}}}$ is almost constant. Finally, the mass balance equation of the mixture reads:

 ${\displaystyle {\frac {-1}{K_{w}}}{\frac {\mathrm {d} p_{w}}{\mathrm {d} t}}+\nabla \cdot \mathbf {v} +\nabla \cdot {\big (}{\boldsymbol {k}}\cdot \left(\nabla p_{w}+\rho _{w}\mathbf {g} \right){\big )}=0}$
(2.98)

### 2.4.2 Balance of linear momentum

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

 ${\displaystyle \nabla \cdot {\boldsymbol {\sigma }}=\rho \,\mathbf {g} }$
(2.99)

where ${\textstyle \rho =\varphi \rho _{w}+(1-\varphi )\rho _{s}}$ is the mixture density and ${\textstyle {\boldsymbol {\sigma }}}$ stands for the total Cauchy stress.

According to the principle of effective stress, the total Cauchy stress tensor, ${\textstyle {\boldsymbol {\sigma }}}$, is equal to the sum of the pore pressure, ${\textstyle p_{w}}$, and the effective stress:

 ${\displaystyle {\boldsymbol {\sigma }}={\boldsymbol {\sigma }}^{\prime }(\mathbf {F} ,V)+p_{w}1}$
(2.100)

where the effective Cauchy stress depends on the solid skeleton deformation through the deformation gradient, ${\textstyle \mathbf {F} }$, and a set of history dependent parameters, ${\textstyle V}$.

The porosity may be further related to the solid skeleton motion. In an infinitesimal volume of the mixture, ${\textstyle \mathrm {d} V}$, the initial void volume is ${\textstyle \varphi _{0}\,\mathrm {d} V}$ whereas the volume of the solid fracion is ${\textstyle (1-\varphi _{0})\,\mathrm {d} V}$. As deformation takes place, the volume of the solid matrix varies according to ${\textstyle \mathrm {d} v=J\,\mathrm {d} V}$. 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 ${\textstyle (1-\varphi _{0})\,\mathrm {d} V}$ whereas the volume of the voids is ${\textstyle \mathrm {d} v-(1-\varphi _{0})\,\mathrm {d} V}$. As such, the porosity may be described by the following expression:

 ${\displaystyle \varphi ={\frac {\mathrm {d} v-(1-\varphi _{0})\,\mathrm {d} V}{J\,\mathrm {d} V}}=1-{\frac {1-\varphi _{0}}{J}}}$
(2.101)

Then, the mixture density, ${\textstyle \rho }$, can be further elaborated to:

 ${\displaystyle \rho =\varphi \rho _{w}+(1-\varphi )\rho _{s}={\frac {(1-\varphi _{0})\rho _{s}+\varphi _{0}\rho _{w}}{J}}+{\frac {J-1}{J}}\rho _{w}={\frac {\rho _{0}}{J}}+{\frac {J-1}{J}}\rho _{w}}$
(2.102)

### 2.4.3 Strong form of the balance equations

Let us state the system of equations governing the mechanics of the byphasic porous media in an updated-Lagrangian form for quasi-static cases along with the corresponding boundary conditions:

 ${\displaystyle {\begin{cases}\nabla \cdot {\boldsymbol {\sigma }}+\mathbf {b} ={\boldsymbol {0}}&{\hbox{ in }}\Omega _{t}\times (0,T)\\{\dfrac {-1}{K_{w}}}{\dot {p}}_{w}+\nabla \cdot \mathbf {v} +\nabla \cdot \mathbf {v} ^{d}=0&{\hbox{ in }}\Omega _{t}\times (0,T)\\{\boldsymbol {u}}({\boldsymbol {X}},t=0)={\boldsymbol {u}}_{0}&{\hbox{ in }}\Omega _{0}\\p_{w}({\boldsymbol {X}},t=0)={p_{w}}_{0}&{\hbox{ in }}\Omega _{0}\\{\boldsymbol {u}}({\boldsymbol {X}},t)={\overline {\boldsymbol {u}}}&{\hbox{ in }}\Gamma _{u}\times (0,T)\\{\boldsymbol {n}}\cdot {\boldsymbol {\sigma }}={\overline {\boldsymbol {t}}}&{\hbox{ in }}\Gamma _{\overline {t}}\times (0,T)\\p_{w}({\boldsymbol {X}},t)={{\overline {p}}_{w}}&{\hbox{ in }}\Gamma _{p_{w}}\times (0,T)\\-{\boldsymbol {n}}\cdot {\boldsymbol {v}}^{d}={\overline {g}}&{\hbox{ in }}\Gamma _{g}\times (0,T)\end{cases}}}$
(2.103)

where ${\textstyle {\dot {p}}_{w}={\frac {dp_{w}}{dt}}}$ is the material time derivative with respect to the solid phase. The boundary of the domain is divided in two parts, ${\textstyle \partial \Omega _{t}=\Gamma _{p_{w}}\cup \Gamma _{g}}$ (${\textstyle \Gamma _{p_{w}}\cap \Gamma _{g}=\emptyset }$), where fixed water pressure ${\textstyle {\overline {p}}_{w}}$ and prescribed water flow ${\textstyle {\overline {g}}}$ are imposed.

### 2.4.4 Weak form

#### Balance of linear momentum

The weak form of the linear momentum balance equation reads:

 ${\displaystyle \int _{\Omega _{t}}\mathbf {w} \cdot \left(\nabla \cdot {\boldsymbol {\sigma }}+\mathbf {b} \right)\,\,\mathrm {d} \Omega _{t}=\int _{\Omega _{t}}\mathbf {w} \cdot \left(\nabla \cdot {\boldsymbol {\sigma }}^{\prime }+\nabla p_{w}+\mathbf {b} \right)\,\,\mathrm {d} \Omega _{t}=0}$
(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:

 ${\displaystyle \int _{\Omega _{t}}{\frac {\partial w_{j}}{\partial x_{i}}}\left(\sigma _{ij}^{\prime }+p_{w}\,\delta _{ij}\right)\mathrm {d} \Omega _{t}=\int _{\Omega _{t}}w_{j}b_{j}\,\,\mathrm {d} \Omega _{t}+\int _{\Gamma _{\overline {t}}}w_{j}\,{\overline {t}}_{j}\,\,\mathrm {d} \gamma _{t}}$
(2.105)

where ${\textstyle \delta _{ij}=(1)_{ij}}$ is the second order identity tensor in index notation.

#### Balance of mass

The weak form of the mass balance equation, Equation 2.98, reads:

 ${\displaystyle \int _{\Omega _{0}}q\left({\frac {-1}{K_{w}}}{\dot {p}}_{w}+\nabla \cdot \mathbf {v} +\nabla \cdot \mathbf {v} ^{d}\right)\mathrm {d} \Omega _{0}=}$ ${\displaystyle =\int _{\Omega _{t}}q\left({\frac {-1}{K_{w}}}{\dot {p}}_{w}+\nabla \cdot \mathbf {v} +\nabla \cdot \mathbf {v} ^{d}\right){\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}=0}$
(2.106)

where ${\textstyle q}$ 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 one-phase problem (the so-called ${\textstyle \mathbf {u} -p}$ 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, ${\textstyle J}$), whereas Larsson and Larsson [64] state the mass balance equation in terms of fluid content (i.e, Equation 2.106 scaled by ${\textstyle J\,\rho _{w}}$).

Again, the divergence theorem is applied to the weak form, that yields:

 ${\displaystyle \int _{\Omega _{t}}q\left({\frac {-1}{K_{w}}}{\dot {p}}_{w}+\nabla \cdot \mathbf {v} \right){\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}-\int _{\Omega _{t}}\nabla q\cdot \mathbf {v} ^{d}\,\,{\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}=\int _{\Gamma _{g}}q\,\,{\overline {g}}\,\,{\frac {1}{J}}\,\,\mathrm {d} \gamma _{t}}$
(2.107)

### 2.4.5 Finite element discretization

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

 ${\displaystyle {\begin{cases}\mathbf {u} \approx \mathbf {u} ^{h}=\mathbf {N} _{u}\cdot {\tilde {\mathbf {u} }}\\\mathbf {w} \approx \mathbf {w} ^{h}=\mathbf {N} _{u}\cdot {\tilde {\mathbf {w} }}\\p_{w}\approx p_{w}^{h}=\mathbf {N} \cdot {\tilde {\mathbf {p} }}_{w}\\q\approx q^{h}=\mathbf {N} \cdot {\tilde {\mathbf {q} }}\\\end{cases}}}$
(2.108)

where ${\textstyle \square ^{h}}$ is the finite element approximation of the field ${\textstyle \square }$ whereas ${\textstyle {\tilde {\square }}}$ are the nodal values. ${\textstyle \mathbf {N} =[N_{1},N_{2},...,N_{n}]}$ and ${\textstyle \mathbf {N} _{u}=[N_{1}1,N_{2}1,...,N_{n}1]}$ are the shape functions whereas ${\textstyle n}$ is the number of nodes. As it can be inferred, the same order shape functions are used for both, solid phase displacements, ${\textstyle \mathbf {u} }$, and water pressure, ${\textstyle p_{w}}$. Since the shape functions are defined in the reference configuration and, thus, time independent, ${\textstyle \mathbf {N} =\mathbf {N} (\mathbf {X} )}$, the following property holds: ${\textstyle {\dot {p}}_{w}^{h}=\mathbf {N} \cdot {\dot {\tilde {\mathbf {p} }}}_{w}}$.

The semi-discrete equations of the hydromechanical formulation, Equations 2.105 and 2.107, are given by:

 ${\displaystyle {\begin{cases}\mathbf {P} ({\boldsymbol {\sigma }}^{\prime })+\mathbf {Q} \cdot {\tilde {\mathbf {p} }}_{w}=\mathbf {f} ^{ext}\\[8pt]\mathbf {Q} ^{\star T}\cdot {\dot {\tilde {\mathbf {u} }}}-{\dfrac {1}{K_{w}}}\mathbf {M} \cdot {\dot {\tilde {\mathbf {p} }}}_{w}-\mathbf {H} \cdot {\tilde {\mathbf {p} }}_{w}=\mathbf {f} ^{p_{w}}\\\end{cases}}}$
(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:

 ${\displaystyle \mathbf {Q} =\int _{\Omega _{t}}\mathbf {B} ^{T}\cdot {\underline {1}}\cdot \mathbf {N} \;\mathrm {d} \Omega _{t}}$
(2.110)
 ${\displaystyle \mathbf {Q} ^{\star T}=\int _{\Omega _{t}}\mathbf {N} ^{T}\cdot {\underline {1}}\cdot \mathbf {B} {\dfrac {1}{J}}\;\mathrm {d} \Omega _{t}}$
(2.111)
 ${\displaystyle \mathbf {M} =\int _{\Omega _{t}}\mathbf {N} ^{T}\cdot \mathbf {N} {\dfrac {1}{J}}\;\mathrm {d} \Omega _{t}}$
(2.112)
 ${\displaystyle \mathbf {H} =\int _{\Omega _{t}}(\nabla \mathbf {N} )^{T}\cdot \mathbf {k} \cdot (\nabla \mathbf {N} ){\dfrac {1}{J}}\;\mathrm {d} \Omega _{t}}$
(2.113)
 ${\displaystyle \mathbf {f} ^{p_{w}}=\int _{\Omega t}(\nabla \mathbf {N} )^{T}\cdot \mathbf {k} \cdot \mathbf {g} \,\,{\dfrac {\rho _{w}}{J}}\;\mathrm {d} \Omega _{t}+\int _{\Gamma _{g}}\mathbf {N} ^{T}{\overline {g}}{\dfrac {1}{J}}\,\mathrm {d} \gamma _{t}}$
(2.114)

### 2.4.6 Temporal discretization

Since quasi-static conditions are assumed, the temporal discretization of the linear momentum balance equation is quite straight-forward. 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:

 ${\displaystyle {\dot {\tilde {\mathbf {u} }}}\approx {\frac {1}{\Delta t}}\left(^{n+1}{\tilde {\mathbf {u} }}-^{n}{\tilde {\mathbf {u} }}\right)={\frac {\Delta {\tilde {\mathbf {u} }}}{\Delta t}}}$ ${\displaystyle {\dot {\tilde {\mathbf {p} }}}_{w}\approx {\frac {1}{\Delta t}}\left(^{n+1}{\tilde {\mathbf {p} }}_{w}-^{n}{\tilde {\mathbf {p} }}_{w}\right)={\frac {\Delta {\tilde {\mathbf {p} }}_{w}}{\Delta t}}}$
(2.115)

Then, introducing the temporal discretization to the semi-discrete equations, Equation 2.109, the problem reads:

 ${\displaystyle {\begin{cases}\mathbf {P} ({^{n+1}{\boldsymbol {\sigma }}}^{\prime })+\mathbf {Q} \cdot ^{n+1}{\tilde {\mathbf {p} }}_{w}=\mathbf {f} ^{ext}\\[8pt]\mathbf {Q} ^{\star T}\cdot \Delta {\tilde {\mathbf {u} }}-{\dfrac {1}{K_{w}}}\mathbf {M} \cdot \Delta {\tilde {\mathbf {p} }}_{w}-\Delta t\,\mathbf {H} \cdot ^{n+1}{\tilde {\mathbf {p} }}_{w}=\Delta t\,\mathbf {f} ^{p_{w}}\\\end{cases}}}$
(2.116)

### 2.4.7 Linearization

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

 ${\displaystyle \mathbf {R} (^{n+1}\mathbf {u} ,^{n+1}\mathbf {p} _{w})={\begin{bmatrix}\mathbf {R} ^{mec}\\\mathbf {R} ^{mass}\end{bmatrix}}}$
(2.117)

where:

 ${\displaystyle \mathbf {R} ^{mec}=\mathbf {P} ({^{n+1}{\boldsymbol {\sigma }}}^{\prime })+\mathbf {Q} \cdot ^{n+1}{\tilde {\mathbf {p} }}_{w}-\mathbf {f} ^{ext}}$
(2.118)

 ${\displaystyle \mathbf {R} ^{mass}=\mathbf {Q} ^{\star T}\cdot \Delta {\tilde {\mathbf {u} }}-{\dfrac {1}{K_{w}}}\mathbf {M} \cdot \Delta {\tilde {\mathbf {p} }}_{w}-\Delta t\,\mathbf {H} \cdot ^{n+1}{\tilde {\mathbf {p} }}_{w}-\Delta t\,\mathbf {f} ^{p_{w}}}$
(2.119)

Then, applying the Newton-Raphson non-linear solver, the linear system of equations equivalent to Equation 2.65 has the form:

 ${\displaystyle {\begin{bmatrix}\mathbf {K} _{uu}&\mathbf {K} _{up_{w}}\\\mathbf {K} _{p_{w}u}&\mathbf {K} _{p_{w}p_{w}}\\\end{bmatrix}}\cdot {\begin{bmatrix}\delta \mathbf {u} \\\delta \mathbf {p} _{w}\end{bmatrix}}=-\mathbf {R} }$
(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:

 ${\displaystyle \int _{\Omega _{t}}{\frac {\partial w_{j}}{\partial x_{i}}}\left(\sigma _{ij}^{\prime }+p_{w}\,\delta _{ij}\right)\mathrm {d} \Omega _{t}=\int _{\Omega _{t}}w_{j}b_{j}\,\,\mathrm {d} \Omega _{t}+\int _{\Gamma _{\overline {t}}}w_{j}\cdot {\overline {t}}_{j}\,\,\mathrm {d} \gamma _{t}}$ ${\displaystyle \int _{\Omega _{t}}q\left({\frac {-1}{K_{w}}}\Delta {p}_{w}+\nabla \cdot \Delta \mathbf {u} \right){\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}-\Delta t\int _{\Omega _{t}}\nabla q\cdot \mathbf {v} ^{d}\,\,{\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}=\Delta t\int _{\Gamma _{g}}q\,\,{\overline {g}}\,\,{\frac {1}{J}}\,\,\mathrm {d} \gamma _{t}}$
(2.121)

where the absence of subscript at displacements and water pressure stands for the value at time ${\textstyle t_{n+1}}$, ${\textstyle \Delta \mathbf {u} =^{n+1}\mathbf {u} -^{n}\mathbf {u} =\mathbf {u} -^{n}\mathbf {u} }$ and ${\textstyle \Delta p_{w}=^{n+1}p_{w}-^{n}p_{w}}$.

#### Linearization of the balance of linear momentum equation

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:

 ${\displaystyle \int _{\Omega _{t}}{\frac {\partial w_{j}}{\partial x_{i}}}\left(\sigma _{ij}^{\prime }+p_{w}\,\delta _{ij}\right)\mathrm {d} \Omega _{t}-\int _{\Omega _{t}}w_{j}\rho \,g_{j}\,\,\mathrm {d} \Omega _{t}=}$ ${\displaystyle =\int _{\Omega _{0}}{\frac {\partial w_{j}}{\partial X_{k}}}\left(S_{kl}^{\prime }+J\,p_{w}C_{kl}^{-1}\,\right)F_{jl}\,\,\mathrm {d} \Omega _{0}-\int _{\Omega _{0}}w_{j}{\big (}\rho _{0}+\left(J-1\right)\rho _{w}{\big )}\,g_{j}\,\,\mathrm {d} \Omega _{0}}$
(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 Piola-Kirchhoff, ${\textstyle \mathbf {S} ^{\prime }}$, is quite straightforward:

 ${\displaystyle \int _{\Omega _{0}}{\frac {\partial w_{j}}{\partial X_{k}}}\delta _{\mathbf {u} }\left(S_{kl}^{\prime }\right)F_{jl}\,\,\mathrm {d} \Omega _{0}=\int _{\Omega _{n+1}}{\frac {\partial w_{j}}{\partial x_{i}}}\mathbb {D} _{ijkl}{\frac {\partial \delta u_{k}}{\partial x_{l}}}\,\mathrm {d} \Omega _{t}}$
(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 single-phase formulation let us now compute the derivative due to the last deformation gradient that appears in the integral of the internal forces:

 ${\displaystyle \int _{\Omega _{0}}{\frac {\partial w_{j}}{\partial X_{k}}}\left(S_{kl}^{\prime }+J\,p_{w}C_{kl}^{-1}\,\right)\delta _{\mathbf {u} }\left(F_{jl}\right)\,\,\mathrm {d} \Omega _{0}=}$ ${\displaystyle =\int _{\Omega _{n+1}}{\frac {\partial w_{j}}{\partial x_{i}}}\left(\sigma _{il}^{\prime }+p_{w}\delta _{il}\right)\delta _{jk}{\frac {\partial \delta u_{k}}{\partial x_{l}}}\,\mathrm {d} \Omega _{t}}$
(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 Cauchy-Green tensor:

 ${\displaystyle \delta _{\mathbf {u} }\left(\mathbf {C} ^{-1}\right)=-\mathbf {C} ^{-1}\cdot \delta _{\mathbf {u} }\left(\mathbf {C} \right)\cdot \mathbf {C} ^{-1}}$
(2.125)

similarly to the temporal derivative of the right Cauchy-Green tensor, Equation 2.22:

 ${\displaystyle \delta _{\mathbf {u} }\left(\mathbf {C} \right)=\,\mathbf {F} ^{T}\cdot {\big (}\nabla \delta \mathbf {u} +(\nabla \delta \mathbf {u} )^{T}{\big )}\cdot \mathbf {F} }$
(2.126)

as a consequence:

 ${\displaystyle \delta _{\mathbf {u} }\left(\mathbf {C} ^{-1}\right)=-\,\mathbf {F} ^{-1}\cdot {\big (}\nabla \delta \mathbf {u} +(\nabla \delta \mathbf {u} )^{T}{\big )}\cdot \mathbf {F} ^{-T}}$
(2.127)

Then, introducing this previous expression and the linearization of the Jacobian, ${\textstyle \delta _{\mathbf {u} }(J)=J\,\nabla \cdot \delta u}$:

 ${\displaystyle \int _{\Omega _{0}}{\frac {\partial w_{j}}{\partial X_{k}}}{\Big (}S_{kl}^{\prime }+\,p_{w}\delta _{\mathbf {u} }\left(JC_{kl}^{-1}\right)\,{\Big )}F_{jl}\,\,\mathrm {d} \Omega _{0}=\int _{\Omega _{n+1}}{\frac {\partial w_{j}}{\partial x_{i}}}p_{w}\left(\delta _{ij}\delta _{kl}-2I_{ijkl}^{4S}\right){\frac {\partial \delta u_{k}}{\partial x_{l}}}\,\mathrm {d} \Omega _{t}}$
(2.128)

On the other hand, differently from the single-phase problem, the external load term due to the gravity also produces a second order term due to the variation of the mixture density:

 ${\displaystyle \delta _{\mathbf {u} }\left(\int _{\Omega _{0}}w_{j}{\big (}\rho _{0}+\left(J-1\right)\rho _{w}{\big )}\,g_{j}\,\,\mathrm {d} \Omega _{0}\right)=\int _{\Omega _{t}}w_{j}\,g_{j}\rho _{w}\delta _{kl}{\frac {\partial \delta u_{k}}{\partial x_{l}}}\,\,\,\mathrm {d} \Omega _{t}}$
(2.129)

where it has been assumed that the water density is almost constant: that is, the term ${\textstyle \delta _{u}{\rho _{w}}=-{\dfrac {\rho _{w}}{K_{w}}}\nabla \delta \mathbf {u} \approx 0}$ 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:

 ${\displaystyle \delta _{p_{w}}\left(\int _{\Omega _{t}}{\frac {\partial w_{j}}{\partial x_{i}}}\left(\sigma _{ij}^{\prime }+p_{w}\,\delta _{ij}\right)\mathrm {d} \Omega _{t}-\int _{\Omega _{t}}w_{j}\rho \,g_{j}\,\,\mathrm {d} \Omega _{t}\right)=}$ ${\displaystyle =\int _{\Omega _{t}}{\frac {\partial w_{j}}{\partial x_{j}}}\,\,\delta {p_{w}}\,\,\mathrm {d} \Omega _{t}}$
(2.130)

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

 ${\displaystyle \delta \mathbf {R} ^{mec}=\int _{\Omega _{n+1}}{\frac {\partial w_{j}}{\partial x_{i}}}{\Big (}\mathbb {D} _{ijkl}+\left(\sigma _{il}^{\prime }+p_{w}\delta _{il}\right)\delta _{jk}+p_{w}\left(\delta _{ij}\delta _{kl}-2I_{ijkl}^{4S}\right){\Big )}{\frac {\partial \delta u_{k}}{\partial x_{l}}}\,\mathrm {d} \Omega _{t}\,\,-}$ ${\displaystyle -\int _{\Omega _{t}}w_{j}\,g_{j}\rho _{w}\delta _{kl}{\frac {\partial \delta u_{k}}{\partial x_{l}}}\,\,\,\mathrm {d} \Omega _{t}-\int _{\gamma }w_{i}{\overline {t}}_{i}{\Big (}\delta \mathbf {u} _{,1}\cdot \left(\mathbf {g} _{2}\times \mathbf {n} \right)-\delta \mathbf {u} _{,2}\cdot \left(\mathbf {g} _{1}\times \mathbf {n} \right){\Big )}\mathrm {d} \gamma \,\,+}$ ${\displaystyle +\int _{\Omega _{t}}{\frac {\partial w_{j}}{\partial x_{j}}}\delta {p_{w}}\mathrm {d} \Omega _{t}}$
(2.131)

#### Linearization of the mass conservation equation

In order to calculate the linearization of the mass balance equation, let us first perform a push-back of Equation 2.121, so that all the quantities are referred to the initial configuration:

 ${\displaystyle \int _{\Omega _{t}}q\left({\frac {-1}{K_{w}}}\Delta {p}_{w}+\nabla \cdot \Delta \mathbf {u} \right){\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}-\Delta t\int _{\Omega _{t}}\nabla q\cdot \mathbf {v} ^{d}\,\,{\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}-\Delta t\int _{\Gamma _{g}}q\,\,{\overline {g}}\,\,{\frac {1}{J}}\,\,\mathrm {d} \Gamma _{t}=}$ ${\displaystyle =\int _{\Omega _{0}}q\left({\frac {-1}{K_{w}}}\Delta {p}_{w}+{\frac {\partial \Delta u_{i}}{\partial X_{j}}}F_{ji}^{-1}\right)\,\,\mathrm {d} \Omega _{0}-\Delta t\int _{\Omega _{0}}{\frac {\partial q}{\partial X_{j}}}F_{ji}^{-1}v_{i}^{d}\,\,\mathrm {d} \Omega _{0}-\Delta t\int _{\Gamma _{g}}q\,\,{\overline {g}}\,\,\mathrm {d} \Gamma _{0}=0}$

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

 ${\displaystyle \delta _{\mathbf {u} }\left(\int _{\Omega _{0}}q\left({\frac {-1}{K_{w}}}\Delta {p}_{w}+{\frac {\partial \Delta u_{i}}{\partial X_{j}}}F_{ji}^{-1}\right)\,\,\mathrm {d} \Omega _{0}-\Delta t\int _{\Omega _{0}}{\frac {\partial q}{\partial X_{j}}}F_{ji}^{-1}v_{i}^{d}\,\,\mathrm {d} \Omega _{0}-\Delta t\int _{\Gamma _{g}}q\,\,{\overline {g}}\,\,\mathrm {d} \Gamma _{0}\right)=}$ ${\displaystyle =\delta _{\mathbf {u} }\left(\int _{\Omega _{0}}q{\frac {\partial \Delta u_{i}}{\partial X_{j}}}F_{ji}^{-1}\,\,\mathrm {d} \Omega _{0}\right)-\delta _{\mathbf {u} }\left(\Delta t\int _{\Omega _{0}}{\frac {\partial q}{\partial X_{j}}}F_{ji}^{-1}v_{i}^{d}\,\,\mathrm {d} \Omega _{0}\right)}$

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:

 ${\displaystyle \delta _{\mathbf {u} }\left(\int _{\Omega _{0}}q{\frac {\partial \Delta u_{i}}{\partial X_{j}}}F_{ji}^{-1}\,\,\mathrm {d} \Omega _{0}\right)-\delta _{\mathbf {u} }\left(\Delta t\int _{\Omega _{0}}{\frac {\partial q}{\partial X_{j}}}F_{ji}^{-1}v_{i}^{d}\,\,\mathrm {d} \Omega _{0}\right)=}$ ${\displaystyle \delta _{\mathbf {u} }\left(\int _{\Omega _{0}}q{\frac {\partial \Delta u_{i}}{\partial X_{j}}}F_{ji}^{-1}\,\,\mathrm {d} \Omega _{0}\right)-\delta _{\mathbf {u} }\left(\Delta t\int _{\Omega _{0}}{\frac {\partial q}{\partial X_{j}}}F_{ji}^{-1}k_{jp}\left({\frac {\partial p_{w}}{\partial X_{q}}}F_{qp}^{-1}+\rho _{w}g_{p}\right)\,\,\mathrm {d} \Omega _{0}\right)}$
(2.132)

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

 ${\displaystyle \delta _{\mathbf {u} }\left(\int _{\Omega _{0}}q{\frac {\partial \Delta u_{i}}{\partial X_{j}}}F_{ji}^{-1}\,\,\mathrm {d} \Omega _{0}\right)=}$ ${\displaystyle \int _{\Omega _{0}}q{\frac {\partial \delta u_{i}}{\partial X_{j}}}F_{ji}^{-1}\,\,\mathrm {d} \Omega _{0}+\int _{\Omega _{0}}q{\frac {\partial \Delta u_{i}}{\partial X_{j}}}\delta _{\mathbf {u} }\left(F_{ji}^{-1}\right)\,\,\mathrm {d} \Omega _{0}}$
(2.133)

The inverse of a matrix derivative may be computed as: ${\textstyle \delta _{\mathbf {u} }\left(\mathbf {F} ^{-1}\right)=-\mathbf {F} ^{-1}\cdot \delta _{\mathbf {u} }\left(\mathbf {F} \right)\cdot \mathbf {F} ^{-1}=\mathbf {F} ^{-1}\cdot \nabla \delta \mathbf {u} }$. Then, introducing this result to Equation 2.133 and performing a push-forward:

 ${\displaystyle \int _{\Omega _{0}}q{\frac {\partial \delta u_{i}}{\partial X_{j}}}F_{ji}^{-1}\,\,\mathrm {d} \Omega _{0}+\int _{\Omega _{0}}q{\frac {\partial \Delta u_{i}}{\partial X_{j}}}F_{jp}^{-1}{\frac {\partial \delta u_{p}}{\partial x_{i}}}\,\,\mathrm {d} \Omega _{0}=\int _{\Omega _{t}}q\left(\delta _{ij}-{\frac {\partial \Delta u_{j}}{\partial x_{i}}}\right){\frac {1}{J}}{\frac {\partial \delta u_{i}}{\partial x_{j}}}\,\,\mathrm {d} \Omega _{t}}$
(2.134)

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

 ${\displaystyle -\delta _{\mathbf {u} }\left(\Delta t\int _{\Omega _{0}}{\frac {\partial q}{\partial X_{j}}}F_{ji}^{-1}k_{jp}\left({\frac {\partial p_{w}}{\partial X_{q}}}F_{qp}^{-1}+\rho _{w}g_{p}\right)\,\,\mathrm {d} \Omega _{0}\right)=}$ ${\displaystyle =-\Delta t\int _{\Omega _{0}}{\frac {\partial q}{\partial X_{j}}}\delta _{\mathbf {u} }\left(F_{ji}^{-1}\right)v_{j}^{d}\,\,\mathrm {d} \Omega _{0}-\Delta t\int _{\Omega _{0}}{\frac {\partial q}{\partial x_{i}}}k_{ip}{\frac {\partial p_{w}}{\partial X_{q}}}\delta _{\mathbf {u} }\left(F_{qp}^{-1}\right)\,\,\mathrm {d} \Omega _{0}}$
(2.135)

where, as commented previously, it has been assumed that the spatial permeability tensor, ${\textstyle \mathbf {k} }$, is constant; other assumptions made on the permeability tensor may induce different large displacement terms. So, the following expression may be obtained:

 ${\displaystyle \Delta t\int _{\Omega _{t}}{\frac {\partial q}{\partial x_{k}}}\left(\delta _{ik}v_{j}^{d}+k_{kj}{\frac {\partial p_{w}}{\partial x_{i}}}\right){\frac {1}{J}}{\frac {\partial \delta u_{i}}{\partial x_{j}}}\,\,\mathrm {d} \Omega _{t}}$
(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:

 ${\displaystyle \delta _{p_{w}}\left(\int _{\Omega _{t}}q\left({\frac {-1}{K_{w}}}\Delta {p}_{w}+\nabla \cdot \Delta \mathbf {u} \right){\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}-\Delta t\int _{\Omega _{t}}\nabla q\cdot \mathbf {v} ^{d}\,\,{\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}-\Delta t\int _{\Gamma _{g}}q\,\,{\overline {g}}\,\,{\frac {1}{J}}\,\,\mathrm {d} \Gamma _{t}\right)=}$ ${\displaystyle =\delta _{p_{w}}\left(\int _{\Omega _{t}}q{\frac {-1}{K_{w}}}\Delta {p}_{w}{\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}-\Delta t\int _{\Omega _{t}}\nabla q\cdot \mathbf {k} \cdot \left(\nabla p_{w}+\rho _{w}\mathbf {g} \right)\,\,{\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}\right)=}$ ${\displaystyle =\int _{\Omega _{t}}q{\frac {-1}{K_{w}}}{\frac {1}{J}}\,\,\delta {p}_{w}\,\,\mathrm {d} \Omega _{t}-\Delta t\int _{\Omega _{t}}\nabla q\cdot \mathbf {k} \cdot \nabla \delta p_{w}{\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}}$

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:

 ${\displaystyle \delta \mathbf {R} ^{mass}=\int _{\Omega _{t}}q\left(\delta _{ij}-{\frac {\partial \Delta u_{j}}{\partial x_{i}}}\right){\frac {1}{J}}{\frac {\partial \delta u_{i}}{\partial x_{j}}}\,\,\mathrm {d} \Omega _{t}+}$ ${\displaystyle +\Delta t\int _{\Omega _{t}}{\frac {\partial q}{\partial x_{k}}}\left(\delta _{ik}v_{j}^{d}+k_{kj}{\frac {\partial p_{w}}{\partial x_{i}}}\right){\frac {1}{J}}{\frac {\partial \delta u_{i}}{\partial x_{j}}}\,\,\mathrm {d} \Omega _{t}+}$ ${\displaystyle +\int _{\Omega _{t}}q{\frac {-1}{K_{w}}}{\frac {1}{J}}\,\,\delta {p}_{w}\,\,\mathrm {d} \Omega _{t}-\Delta t\int _{\Omega _{t}}\nabla q\cdot \mathbf {k} \cdot \nabla \delta p_{w}{\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}}$
(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.

#### Matrix form of the linearization

Finally, the matrix expression of the terms defined in Equation 2.120 may be obtained as:

 ${\displaystyle \mathbf {K} _{uu}=\int _{\Omega _{t}}\mathbf {B} ^{T}\cdot \left({\underline {C}}+p_{w}\,\left({\underline {1}}\otimes {\underline {1}}-2\,{\underline {\mathbf {I} }}^{4S}\right)\right)\cdot \mathbf {B} \,\,\mathrm {d} \Omega _{t}+}$ ${\displaystyle +\int _{\Omega _{t}}\mathbf {B} _{NL}^{T}\cdot {\underline {\underline {\boldsymbol {\sigma }}}}\cdot \mathbf {B} _{NL}\,\,\mathrm {d} \Omega _{t}\,\,-\int _{\Omega _{t}}\rho _{w}\,\,\mathbf {N} _{u}^{T}\cdot \mathbf {g} \otimes {\underline {1}}\cdot \mathbf {B} \,\,\mathrm {d} \Omega _{t}}$
(2.138)
 ${\displaystyle \mathbf {K} _{up_{w}}=\int _{\Omega _{t}}\mathbf {B} ^{T}\cdot {\underline {1}}\cdot \mathbf {N} \,\,\mathrm {d} \Omega _{t}=\mathbf {Q} }$
(2.139)

 ${\displaystyle \mathbf {K} _{p_{w}u}=\int _{\Omega _{t}}\mathbf {B} ^{T}\cdot {\underline {1}}\cdot \mathbf {N} \,\,{\frac {1}{J}}\,\,\mathrm {d} \Omega _{t}+\mathbf {K} _{p_{w}u}^{geo}}$
(2.140)
 ${\displaystyle \mathbf {K} _{p_{w}p_{w}}=-\int _{\Omega _{t}}{\dfrac {1}{K_{w}}}\mathbf {N} ^{T}\cdot \mathbf {N} \,\,\mathrm {d} \Omega _{t}-\Delta t\int _{\Omega _{t}}\left(\nabla \mathbf {N} \right)^{T}\cdot \mathbf {k} \cdot \nabla \mathbf {N} \,\,\mathrm {d} \Omega _{t}\,\,=}$ ${\displaystyle =-{\dfrac {1}{K_{w}}}\mathbf {M} -\Delta t\mathbf {H} }$
(2.141)

where ${\textstyle {\underline {\square }}}$ is the Voigt notation of the tensor ${\textstyle \square }$, and the definition of terms involved in the computation of the non-linear stiffness matrix, ${\textstyle \mathbf {B} _{NL}}$ and ${\textstyle {\underline {\underline {\boldsymbol {\sigma }}}}}$, may be found in Equations 2.80 and 2.81 for two-dimensional problems and, for three-dimensional cases, in [60]. Due to the complexity, the geometrical terms of the mass balance equation that appear in Equation 2.140, ${\textstyle \mathbf {K} _{p_{w}u}^{geo}}$, 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.

### 2.4.8 Stabilization of the hydro-mechanical problem

Undrained conditions in water-saturated soils with nearly incompressible constituents results in quasi-incompressible behavior. Therefore, as the problem approaches zero permeability and incompressilibity of the soil and water constituents, the system of the discretized equations describing the ${\textstyle \mathbf {u} -p_{w}}$ formulation has the similar structure to that found when using a mixed ${\textstyle \mathbf {u} -p}$ formulation of Solid Mechanics problems [66]. Depending on the interpolants used to discretize both nodal variables, the problem may become ill-posed from a mathematical point of view, which may result in non-uniqueness and mesh-dependence of the solution, the well-known 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 non-compliance of the Babuska-Brezzy conditions or the patch test due to an improper finite-dimensional space in the finite element discretization.

It is well known that in the hydro-mechanical problem, the pore pressure field tends to exhibit oscillations, which tend to increase when the time step is reduced [67]. As shown in the one-dimensional (Terzhagi's equation) analysis of the discrete equations, the use of implicit time-integration 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:

 ${\displaystyle \Delta t\leq {\frac {h^{2}}{6\,c_{v}}}}$
(2.142)

where ${\textstyle c_{v}}$ is the coefficient of consolidation and ${\textstyle h}$ the element size. By numerical examples, it can be shown that this limit also applies to two- and three-dimensional 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 equal-order, low-order 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:

 ${\displaystyle {\begin{cases}\mathbf {P} ({^{n+1}{\boldsymbol {\sigma }}}^{\prime })+\mathbf {Q} \cdot ^{n+1}{\tilde {\mathbf {p} }}_{w}=\mathbf {f} ^{ext}\\[8pt]\mathbf {Q} ^{\star T}\cdot \Delta {\tilde {\mathbf {u} }}-{\dfrac {1}{K_{w}}}\mathbf {M} \cdot \Delta {\tilde {\mathbf {p} }}_{w}-\Delta t\,\mathbf {H} \cdot ^{n+1}{\tilde {\mathbf {p} }}_{w}+{\underline {\mathbf {M} ^{s}\cdot \Delta {\tilde {\mathbf {p} }}_{w}}}=\Delta t\,\mathbf {f} ^{p_{w}}\\\end{cases}}}$
(2.143)

where ${\textstyle \mathbf {M} ^{s}}$ is the stabilization matrix, whose definition depends on the employed stabilization technique.

#### Polynomial pressure projection

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 ${\textstyle \mathbf {u} -p}$ 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:

1. A mixed equal order interpolation of the scalar and vector fields.
2. A ${\displaystyle L_{2}}$ projection of the scalar variables (volume or pressure variables).

The method is obtained by modifying the mixed variational equation (i.e. the pressure continuity equation) by using local ${\textstyle L_{2}}$ 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 equal-order approximations and leads to a stable variational formulation. Unlike other stabilization methods, the Polynomial Pressure Projection does not require the calculation of higher-order 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 ${\textstyle p\in L_{2}}$, the ${\textstyle L_{2}}$ projection operator ${\textstyle {\breve {p}}:L_{2}\rightarrow Q^{0}}$ is defined by

 ${\displaystyle \int _{\Omega _{0}}{\breve {q}}\,(p-{\breve {p}})\,d{\Omega _{0}}=0\;\;\;\forall {\breve {q}}\in Q^{0}}$
(2.144)

where ${\textstyle {\breve {p}}}$ is the best approximation of the pressure ${\textstyle p}$ in the space of polynomials of order ${\textstyle {\mathcal {O}}(Q^{0})}$.

 ${\displaystyle \int _{\Omega _{0}^{e}}(q-{\breve {q}})\,\alpha \,(p-{\breve {p}})\,d{\Omega _{0}}=\int _{\Omega _{0}^{e}}\alpha \,(q\,p-{\breve {q}}\,{\breve {p}})\,d{\Omega _{0}}}$
(2.145)

where ${\textstyle \alpha }$ is the stabilization parameter.

Finally, the discrete stabilization matrix, ${\textstyle \mathbf {M} ^{s}}$, may be expressed:

 ${\displaystyle \mathbf {M} ^{s}=\int _{\Omega _{t}}\alpha \,\mathbf {N} ^{T}\cdot \mathbf {N} {\dfrac {1}{J}}\;d\Omega _{t}-\int _{\Omega _{t}}\alpha \,{\breve {\mathbf {N} }}^{T}\cdot {\breve {\mathbf {N} }}{\dfrac {1}{J}}\;d\Omega _{t}}$
(2.146)

where ${\textstyle {\breve {\mathbf {N} }}}$ are the set of polynomials introduced in Equation 2.144; in the case of linear triangles, these local element polynomials are ${\textstyle {\breve {\mathbf {N} }}^{e}=[1/3,\;1/3,\;1/3]}$.

#### Estimation of the stabilization factor

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 time-step for implicit methods is used.

The rational of the method proposed by [74] is quite simple and is based on the solution of the one-dimensional 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:

• The nodal water pressure should be greater than the initial one at every time step.
• At a given time, the water pressure field should monotonically increase or decrease along the bar.

By using this rational, the authors obtained that the critical time step should be:

 ${\displaystyle \Delta t\geq {\frac {h^{2}}{6\,c_{v}}}}$
(2.147)

where ${\textstyle c_{v}={\dfrac {M\,k^{\prime }}{g\,\rho _{w}}}=M\,k}$ is the coefficient of consolidation, ${\textstyle M}$ is the constrained modulus and ${\textstyle k^{\prime }}$ is the permeability in velocity units, whereas ${\textstyle k=k^{\prime }\,g\,\rho _{w}}$

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: