## Abstract

The objective of this monograph is the derivation and implementation of a robust Finite Element formulation for the solution of solid-pore fluid coupled problems in multi-fractured porous media.

A coupled displacement-pore pressure FEM formulation for solving solid-pore fluid interaction problems is first introduced. The interaction between both components is governed by two equations: the balance of momentum for the mixture solid-fluid and the mass balance for the pore fluid.

Under nearly undrained-incompressible conditions, such formulation suffers from instability problems because of the violation of Babuska-Brezzi conditions. In order to work with elements of equal order interpolation for the displacement and pore pressure, the formulation is stabilized by means of the Finite Increment Calculus method (FIC). The FIC-stabilized formulation is tested against stable elements with a higher order interpolation for the displacement field in 2D and 3D examples.

Continuum damage mechanics is the basis of the crack growth strategy for the proposed fracture propagation technique. The strain softening models used for quasi-brittle materials favour spurious strain localization and ill-posedness of the boundary value problem if the damage variable only depends on the strain state at the point under consideration.

An integral-type non-local damage model associated to a characteristic length parameter is presented as a method to control the size of the fracture process zone and fully regularize the problem. Two examples are solved assessing the robustness of the model in front of changes in the spatial discretization.

Quasi-zero-thickness interface elements are formulated to represent discontinuities in the porous domain. A bilinear cohesive fracture model is used to describe its mechanical behaviour, and a formulation derived from the cubic law models the fluid flow through the crack.

Finally, a new methodology for the simulation of fracture propagation processes in saturated porous media is presented. The non-local damage model is used in conjunction with the interface elements to predict the degradation pattern of the domain and insert new fractures followed by remeshing. Fluid-driven fracture propagation examples in 2D and 3D are presented to illustrate the accuracy of the proposed technique.

## Acknowledgements

After three years of hard work I can say that most of the objectives set at the beginning have been achieved. The time spent in the research and elaboration of this monograph has been truly rewarding for the acquired knowledge and for all the gained friends. Indeed, this project would have been impossible without all the people around me that have helped me in some way.

First, I would like to express my sincere gratitude to Prof. Eugenio Oñate for offering me the opportunity of working at CIMNE in this monograph. As my supervisor he always gave me wise advice and helped me work better, but he also gave me freedom to think for myself and let me grow as a researcher, which I really appreciate. I also want to thank him along with Sonia Sagrista for encouraging me to go on a three-month secondment at Tsinghua University (Beijing) in the framework of the TCAINMAND project. I can say that it has been one of the best experiences in my life.

I specially want to thank all the officemates at CIMNE that have been close to me kindly clarifying all my doubts, and with whom I enjoyed several meals during these years: Salva Latorre, Pablo Agustín Becker, Ferran Arrufat, Merce López, Guillermo Casas, Lorenzo Gracia, Miguel Ángel Celigueta, Miquel Santasusana and Roberto Flores.

I am deeply grateful to Prof. Zhuo Zhuang as the supervisor of my secondment in Tsinghua University. Although I was very far from home, he always made me feel comfortable. I really appreciate all the interesting discussions we had every Thursday during the group meetings. In the same way, I would like to thank Prof. Zhanli Liu and all the officemates, specially: Tao Wang, Qinglei Zeng, Yue Gao, Liu Fengxian, Liyuan Wang, Ramiro Mena and Jan Rojek. I am also very thankful to my room mate in Tsinghua Jaime Santirso, and to Duan Wenjie.

I would also like to thank other colleagues at CIMNE that have helped me in various matters: Fernando Salazar, Enrique Escolano, Miguel Pasenau, Javi Gárate, Josep Maria Carbonell, Pooyan Dadvand, Riccardo Rossi, Jordi Cotela, Carlos Roig, Rubén Zorrilla and Vicente Mataix.

I also acknowledge the Agencia de Gestió d'Ajuts Universitaris i de Recerca (AGAUR) for the financial support.

Finally, this work is dedicated to my family: Lluís, Marta, Oriol, MªAngels and Gabriel. They are the most important people in my life and represent a fundamental and unique support.

# 1 Introduction

## 1.1 Motivation and presentation

Problems involving the flow of a fluid through a porous medium have always been the object of special interest in geomechanics. Many geomaterials contain pores and other cavities that are filled with fluid in saturated or unsaturated conditions. The pressure of such fluid and the deformation of the solid medium are closely related in a challenging problem that must be analysed in a coupled manner.

Further complexity is added when one takes into consideration fluid-driven fracture propagation in the solid domain. Indeed, depending on the fluid pressure distribution and the external efforts acting over the domain, not only pre-existing discontinuities can open and close, but they can also propagate and conform new fractures, introducing an additional degree of coupling.

From the classical contributions of Terzaghi [1] and Biot [2] important effort has been made in the last decades to develop competent numerical methods in poromechanics that allow analysing and understanding the complexity of the mechanisms involved in problems with a coupling between a solid skeleton and a fluid phase. Such techniques are becoming essential tools of for different fields and applications.

In hydrogeology, for instance, these numerical methods are important to represent the water flow through porous soil when characterizing aquifers performance (storage capacity, water transmissivity, water content, etc.), and also to study the dispersion and transport of dissolved contaminants [3].

For soil mechanics, these tools are specially interesting for solving problems in which the seepage flow and pressures interact with the dynamic behaviour of the solid skeleton (Figure 1a).

 (a) Seepage through earth dam. (b) Simplified scheme of an hydraulic fracturing process. Figure 1: Examples of problems with important interaction between a solid skeleton and a fluid phase.

An important field of application is in petroleum engineering. Here, the oil-gas-soil interaction takes the leading role along with the hydraulic fracturing (Figure 1) as a common technique to enhance reservoir permeability and well efficiency [4,5].

Other applications for numerical methods in poromechanics have recently emerged for a variety of fields, including the underground storage of carbon dioxide [6], the geothermal energy production [7], the analysis of interstitial flow in bone tissue [8], and even the study of fractures in epithelial cell sheets [9,10].

There are clearly different applications for this kind of problem, but in all of them one can always distinguish, at least, two distinct components: a solid skeleton and a fluid phase. The solid medium is usually composed by particles of different shapes and sizes enclosing little cavities that can be totally or partially filled by the fluid.

It seems natural that a proper analysis of these kind of problems should require taking into account both parts, i.e. solid and fluid, for the derivation of the governing equations and the definition of the boundary conditions.

Thereby, when solving this coupled problem, it will be necessary to capture the seepage of the fluid through the porous medium, and the deformations of the solid skeleton at any point of the domain. From the numerical perspective, one of the most widespread approaches is the displacement-pressure formulation, in which the main unknowns are the displacements of the solid skeleton ${\textstyle {\boldsymbol {u}}}$ and the pore pressure ${\textstyle p}$.

However, in the limit of nearly incompressible pore fluid and small permeability, coupled poromechanics formulations suffer from instability problems. Finite elements exhibit locking in the pressure field and spurious oscillations appear due to the violation of the so-called Babuska-Brezzi conditions [11,12]. The oscillations can be overcome by locally refining the mesh or by using shape functions of the displacement field one order higher than those of the pressure field. In practical applications this is, however, not feasible because of the increment in the computational cost. In this sense, stabilization methods have been found to be powerful tools to circumvent the Babuska-Brezzi conditions violation without compromising the efficiency of the computation.

Discontinuities in porous materials such as concrete, soil and rock have a noticeable influence in the mechanical and hydraulic behaviour of a multi-phase system. They represent a jump in the displacement field as well as in the stress and strain fields. The presence of these fracture planes modifies the global material response in terms of strength and stiffness, and introduces directional preferences in the deformation of the solid and in the fluid flow. It is thereby essential to properly describe the initiation and propagation of these discontinuities with a suitable fracture mechanics approach.

In the context of the Finite Element Method (FEM), modelling of discontinuities can be done by means of interface elements. These special elements, which can be used to define pre-existing or propagating cracks, act as joints between the discontinuities in the displacement field and can be used to represent the enhanced conductivity of the medium.

Thereby, in the modelling of the fractured domain, it is necessary to implement a constitutive model that allows determining the initiation and propagation of new failure zones, and a model characterizing the mechanical behaviour of the fractures. In this regard, the combination of a “smeared crack” model with a “discrete crack” approach permits working with pre-defined discontinuities and, at the same time, insert new joint elements when necessary.

Concerning the propagation of fractures, suitable regularization methods and remeshing techniques become essential to minimize the mesh dependency when determining the direction of the crack growth.

### Kratos framework

This work has been implemented in Kratos Multiphysics [13]. Kratos is an Open-Source framework for building parallel multi-disciplinary simulation software. Modularity, extensibility and High Performance Computing (HPC) are some of its main objectives. Kratos has BSD licence and is written in C++ with extensive Python interface.

It features a “core and applications approach”, where standard tools (databases, linear algebra, solvers, etc.) come as part of the “core” and are available as building blocks in the development of “applications”, which derive from the core and focus on the solution of particular problems of interest (see Figure 2). We note here that all the implementations of this monograph gave as a result the “Poromechanics Application” in Kratos.

 Figure 2: Kratos' “core and applications approach”.

## 1.2 Objectives

The main objective of this work is to advance in the development of finite element methods for the analysis of coupled solid-pore fluid problems in which the initiation and propagation of cracks represents a key issue and is strongly influenced by the interaction between the different media. To this end, the following goals were set for the three-years of work:

• Derivation and testing of a stable FEM formulation for coupled solid-pore fluid problems under fully saturated conditions for 2D and 3D analysis.
• Implementation and validation of a robust continuum damage model for the analysis of degradation of quasi-brittle materials.
• Development of special interface elements compatible with the implemented formulation, and valid for the definition of pre-existing fractures in the porous media as well as propagating cracks.
• Development of a crack tracking technique based on the analysis of the damage pattern around the crack tip.
• Implementation of an efficient technique to insert new interface elements in the domain and remesh the model every time the crack propagates.
• Application of the previous implementations for solving a problem of interest in engineering. An example of a problem where all the components are required is the analysis of hydraulic fracture processes.

## 1.3 State of the art

### 1.3.1 Coupled solid-pore fluid formulation

The first work regarding the coupling between a fluid and a porous medium can be attributed to Terzaghi, who proposed a one-dimensional theory of consolidation in 1943 [1]. The generalization to three-dimensional cases was presented by Biot who also described methods for determining the elastic coefficients of the theory [2,14,15]. Although these first works were restricted to linear elastic materials, they have been the basis for a lot of subsequent works in geophysics, soil and rock mechanics.

The transition between Biot's theoretical formulation and the first successful numerical results is found in Ghaboussi and Wilson's work, in which the assumption of incompressibility of fluid and solid particles was relaxed [16]. Zienkiewicz and Bettes [17] proposed a formulation based only on the displacement of solid and the pore pressure, which is suitable for the low-frequency situations such as earthquake problems.

De Boer and Kowalski [18] developed a plasticity theory for saturated porous media and, due to the increasing interest in non-linear applications, Zienkiewicz and co-workers expanded the theory to a generalized incremental form for non-linear material behaviour and large strain effects for liquefaction analysis of soil structures [19,20].

In order to overcome the violation of Babuska-Brezzi conditions, several stabilization techniques can be found in the literature in the context of fluid mechanics and solid mechanics, and have been also extended to poromechanics problems.

Some of these techniques are based on the pioneer work of Chorin [21] and use special time stepping schemes to generate stabilization terms [22,23,24].

Other types of approaches have also achieved good stabilized results by modifying the fluid mass balance equation with the addition of a characteristic term that needs to be calibrated [25,26,27,28]. Here one can also include the methods based on the concept of Polynomial-Pressure-Projections, in which the stabilizing terms use element-local projections of the pressure field to counteract the inherent instabilities [29,30,31,32].

There is an additional group of stabilization techniques that make use of the Finite Increment Calculus method (FIC) developed by Oñate and co-workers. In this approach, stabilization terms are derived by expanding the residual of the governing equations in a Taylor series. Previously stabilized FEM formulations for quasi and fully incompressible fluids and solids were based on the first order form of the FIC balance equation in space [33,34,35,36], while recent works are based on the second-order FIC form of the mass balance equation [37,38].

### 1.3.2 Fracture modelling

In the last decades various authors have developed numerical methods for the initiation and propagation of fractures in the context of porous media.

One approach that has obtained notable attention in the last years is the extended finite element method (XFEM) developed for crack analysis [39,40,41,42,43,44]. The XFEM handles discontinuities in an approximation space with jump functions, i.e. local enrichment functions, and so there is no need to represent cracks in a mesh. This considerably reduces the cost of remeshing during crack growth, but in return, the method demands a higher computational cost in terms of number of degrees of freedom and numerical integration.

The other obvious approach that has shown significant impact includes the group of numerical techniques purely based on the finite element method (FEM). In this group, numerous methods can be found in the literature, but probably two main subgroups can be distinguished: the “smeared crack” approaches, continuum based methods in which the influence of developing fractures is incorporated into the constitutive stress-strain law [45,46,47,48], and the “discrete crack” models, in which each single discontinuity is represented explicitly [49,50,51,52].

For the analysis of material degradation and crack initiation in porous media, most of the techniques developed for continuum damage mechanics are applicable.

The concept of damage was introduced by Kachanov in 1958 in the context of creep rupture in metals [53]. Rabotnov [54] gave the problem physical meaning by suggesting that the reduction of the sectional area was measured by means of the damage parameter.

The thermodynamic formalism involved in the irreversible process of damage was developed by Lemaitre and Chaboche [55], and other important contributions to our knowledge were given by: Mazars and Pijaudier-Cabot [56], Simo and Ju [57], Oller et al. [58], Oliver et al. [59,60], and Cervera et al. [61,62], to name but a few.

If the damage parameter depends only on the strain state at the point under consideration, numerical simulations exhibit a pathological mesh dependence and the energy consumed by the system tends to zero as the mesh is refined. Classical constitutive models require an extension in the form of a characteristic length to properly model the thickness of localized zones.

Such extension can be introduced by means of second gradient models [63], micro-polar [64], strain gradient [65], viscous [66] and non-local terms [67].

Non-local elasticity was developed by Eringen [68] who later extended it to non-local elasto-plasticity [69]. Subsequently, it was found that certain non-local formulations could act as efficient localization limiters with a regularizing effect on problems with strain localization [70].

Modelling of discontinuities between two bodies in contact has been performed via contact elements in solid mechanics [71,72,73], and the idea has also been used in fracture mechanics to model cohesive forces by means of cohesive elements [74,75].

Since Goodman et al. proposed the “zero-thickness” interface element to describe the mechanical behaviour of pre-existing joints in rock masses [76], many authors have developed strategies to adapt that element for the solution of fracture processes in solid-pore fluid coupled problems.

Similarly as in the formulation of the continuum solid-pore fluid mixture, two equilibrium equations usually govern the coupled behaviour of the fractures. One equation deals with the mechanical behaviour of the crack, whereas the other equation describes the balance of fluid mass within the fracture.

 (a) Single-noded zero-thickness interface element. (b) Double-noded zero-thickness interface element. (c) Triple-noded zero-thickness interface element. Figure 3: Two-dimensional zero-thickness interface elements.

Concerning the modelling of the fluid flow through the discontinuity, three different types of zero-thickness interface elements can be found in the literature: single, double and triple noded elements (Figure 3). The single-noded elements are the simplest ones and represent pipe elements that are superimposed onto the standard continuum element edges [77]. With this elements the hydraulic head is uniform across the discontinuity and thus only the longitudinal conductivity is taken into account. The triple-noded elements allow including the influence of a transversal conductivity through the discontinuity [78]. The two nodes at each side of the interface represent the potentials in the pore pressure between the two sides of the discontinuity. The third node, placed at the middle of the interface, stores the average potential of the longitudinal fluid through the fracture. Finally, the double-noded elements take into account both types of conductivity but the middle node variable is eliminated in terms of the variables of the external nodes. Ng and Small [79] used this double-noded zero-thickness interface element to model flow problems with pre-existing discontinuities but with no hydraulic potential drop between the two interface walls. Segura and Carol [80] introduced the transversal conductivity for zero-thickness elements to account for the exchange of fluid between the discontinuity and the porous material.

For analysing the mechanical behaviour of fractures, we find essentially two different approaches: those based on linear elastic fracture mechanics (LEFM), and those based on non-linear fracture mechanics (NLFM). LEFM techniques were first proposed to solve fracture propagation problems by means of remeshing without considering a fracture process zone (FPZ) before the crack tip. This approach is applicable in large structures where the size of the FPZ is negligible. However, for quasi-brittle analyses, the consideration of a non-linear fracture process zone where the energy is dissipated before it completely fails was found to be essential. In those cases the NLFM technique is usually applied and a softening law relates the cohesive stress to the crack opening in the FPZ. The first procedure based on the cohesive fracture model was originally introduced by Barenblatt [81,82] for brittle materials and by Dugdale [83] for plastic materials. Hillerborg et al. [84] developed the first fictitious crack model for Mode I fracture. It was extended later for the mixed mode fracture, from which Camacho and Ortiz [85] proposed a suitable fracture criterion that is widely used in the literature.

A crucial part in the modelling of fracture propagation is the criterion for determining the direction of the crack growth. Some criteria are based on the local evaluation of the stress field at the crack tip, e.g. the maximum circumferential stress criterion [86]. Others are based on the energy distribution at the fractured zone, such as the minimum strain energy density criterion [87] or the maximal strain energy release rate criterion [88]. There is also another group of crack growth criteria based on continuum damage mechanics [89].

Schrefler et al. [90] and Secchi et al. [91] implemented the double-noded interface element with cohesive fracture model, working together with an adaptive mesh refinement technique in porous materials. Also, Khoei et al. implemented a model of cohesive crack growth using an adaptive mesh refinement procedure for two-dimensional problems [92]. In that work, a bilinear cohesive law originally proposed by Espinosa and Zavattieri [93] is adopted, and a modified superconvergent patch recovery (SPR) technique is used for the error estimation.

More recently, Moës et al. presented an alternative technique combining fracture mechanics with level set procedures [94,95].

## 1.4 Organization

The book is organized as follows.

In Chapter 2 the displacement-pore pressure formulation is described, starting from the governing equations of the problem and the fully coupled system of equations to be solved. Following that, the instability produced under undrained-incompressible conditions is presented and the stable formulation based on the second-order FIC form of the mass balance equation in space is thoroughly detailed. Two academic examples are solved at the end of the chapter comparing the FIC-FEM elements of equal order interpolation to stable elements with a higher order interpolation for the displacement field.

Chapter 3 starts introducing the essential concepts on continuum damage mechanics. Details are given on the fundamental concepts of the isotropic damage theory, including the equivalent strain forms and damage evolution laws that have been implemented in this work. Next, the regularization techniques developed to overcome the problems associated to strain localization are illustrated. Special attention is given to the integral-type non-local damage model, highlighting the most relevant aspects of its numerical implementation. Two examples are shown at the end, where we assess the robustness of the implemented damage models against changes in the spatial discretization.

In Chapter 4 the formulation of the developed quasi-zero-thickness interface elements is detailed, explaining the mechanical behaviour of the fracture and the model governing the fluid flowing in it. After that, we present the proposed fracture propagation technique combining the non-local damage model with the interface elements. Finally, two plane-strain examples are solved to test the accuracy of the crack propagation methodology, and one additional case is included to show the performance of the generalized 3D formulation.

Chapter 5 summarizes the most relevant conclusions of this work, pointing out its main contributions. In the end, the new perspectives to future works are briefly commented.

## 1.5 Related publications

The following publications have emanated from the work carried out in this monograph. The list is ordered chronologically.

• I. de Pouplana and E. Oñate. Combination of a non-local damage model for quasi-brittle materials with a mesh-adaptive finite element technique. Finite Element in Analysis and Design, vol. 112, pp. 26-39, 2016. DOI: 10.1016/j.finel.2015.12.011
• I. de Pouplana and E. Oñate. A FIC-based stabilized mixed finite element method with equal order interpolation for solid-pore fluid interaction problems. International Journal for Numerical and Analytical Methods in Geomechanics, vol. 41, pp. 110-134, 2016. DOI: 10.1002/nag.2550
• I. de Pouplana, L. Gracia, F. Salazar and E. Oñate. Cracking of a concrete arch dam due to seasonal temperature variations. Theme A - ${\displaystyle 14^{th}}$ International Benchmark Workshop on the Numerical Analysis of Dams, 2017.
• L. Gracia, I. de Pouplana, F. Salazar and E. Oñate. Static and seismic analysis of an arch-gravity dam. Theme B - ${\displaystyle 14^{th}}$ International Benchmark Workshop on the Numerical Analysis of Dams, 2017.
• I. de Pouplana and E. Oñate. Finite element modelling of fracture propagation in saturated media using quasi-zero-thickness interface elements. Article sent to Computers and Geotechnics, 2017.
• E. Oñate and I. de Pouplana. A methodology for analysis of delamination in composites using joint elements. In preparation.

# 2 Coupled Solid-Pore Fluid Interaction Problem

## 2.1 Introduction

Sometimes two or more physical systems interact with each other in such a way that it is not possible to obtain a solution independent of any of the systems without considering all the systems at the same time. It is said that those systems belong to a coupled problem, with the coupling being stronger or weaker depending on the level of interaction between each system.

Our problem of interest belongs to the type of coupled problems in which the different domains of the physical systems overlap. In such cases, the coupling occurs through the governing differential equations describing the different physical phenomena. In particular, we are going to study solid-pore fluid interaction problems in which a unique mesh represents a porous medium.

The problem falls in the field of poromechanics [96], a branch of physics that studies porous materials whose mechanical behaviour is significantly influenced by the pore fluid.

A porous medium is composed of a matrix, containing solid components of different shape and size, and a porous space, that can be totally or partially filled by a fluid. Thereby, any infinitesimal volume of a porous medium can be treated as the superimposition of two material particles: the skeleton particle and the fluid particle (Figure 4).

 Figure 4: Infinitesimal volume of porous medium as a combination of two phases: a solid skeleton and a fluid medium.

A continuous description of such an heterogeneous medium requires taking into account the mechanical coupling between the solid skeleton and the fluid phase when formulating the governing equations of the problem.

In the last decades important effort has been made in the development of powerful numerical methods for poromechanics that allow analysing porous media considering the mechanical coupling between the solid and the fluid phase.

A one-dimensional theory of consolidation was first proposed by Terzaghi in 1943 [1] and then Biot extended it to three-dimensional cases [2,14,15]. Although these first works were restricted to linear elastic materials, they have been the basis for a lot of subsequent works in geophysics, soil and rock mechanics.

The first numerical solution of Biot's formulation was obtained by Ghaboussi and Wilson [16] and the work was further developed by Zienkiewicz et al. [97,17]. Later on, due to the increasing interest in non-linear applications, Zienkiewicz and co-workers expanded the theory to a generalized incremental form for non-linear materials and large deformation problems [19,20].

The mathematical formulation of solid skeleton and pore fluid interaction presented here is based on the model proposed by Zienkiewicz et al. [97]. The problem was originally formulated for fully saturated conditions in terms of the solid matrix displacement ${\textstyle u_{i}}$, the mean fluid velocity relative to the solid phase ${\textstyle w_{i}}$ and the fluid pore pressure ${\textstyle p}$, but in many geo-mechanical problems with no high-frequency phenomena involved, the fluid relative velocity ${\textstyle w_{i}}$ can be neglected and so the equations can be simplified to work with only two main variables: ${\textstyle u_{i}}$ and ${\textstyle p}$ [98].

In the limit of nearly incompressible pore fluid and small permeability, the coupled poromechanics formulations suffer from instability problems. Finite elements exhibit locking in the pressure field and spurious oscillations in the numerical solution for the pressure appear due to the violation of the so-called Babuska-Brezzi conditions [11,12]. The oscillations can be overcome by locally refining the mesh and by using shape functions of the displacement field one order higher than those of the pressure field. In practical applications this is, however, not the best approach because of the increment in the computational cost. In this sense, stabilization methods have been found to be powerful tools to circumvent the Babuska-Brezzi conditions violation without compromising the efficiency of the computation.

Several stabilization techniques have been investigated in the past years in the context of computational fluid dynamics [21,25,26,22,33,34] and (incompressible) solid mechanics [99,100,28,35,36], and have been also extended to poromechanics problems [27,23,24,101,29].

Although each stabilization approach has its differential aspects, they can be classified in three main categories.

The first category comprises those techniques in which special time stepping schemes are applied in order to generate stabilization terms. Probably the earliest work in this category is due to Chorin [21] who proposed a technique to deal with incompressible fluid problems which is now referred as the fractional step method or the operator-splitting method. Such a staggered time stepping algorithm has been found to provide stabilization in the steady-state when used in a convenient form [21,22,23,24].

The second type of techniques are more direct stabilization methods based on the perturbation of the fluid mass conservation equation. Instead of using special time-stepping algorithms that give additional terms in the steady-state approximation, the fluid balance equation is modified by adding a stabilization term multiplied by a parameter that needs to be carefully calibrated [25,26,27,28]. This group also includes the methods based on the concept of Polynomial-Pressure-Projections, in which the additional stabilizing terms use element-local projections of the pressure field to counteract the inherent instabilities [29,30,31,32].

The third category considered here is the Finite Increment Calculus method (FIC) [33,34,35,36,37,38], which is the approach also adopted in the present work. The FIC technique is based on expressing the equations of balance of mass and momentum in a space/time domain of finite size, and retaining higher order terms in the Taylor series expansion typically used for expressing the change in the transported variables within the balance domain. Apart from the standard terms of infinitesimal theory, the FIC form of the balance equations contains derivatives of the classical differential equations multiplied by characteristic distances in space and/or time.

Previously stabilized FIC-FEM formulations were based on the first-order form of the FIC balance equation in space and can be found in the literature for quasi and fully incompressible fluids and solids [33,34,35,36], and even for 1D and 2D poromechancis problems [102]. In this work, the more recent second-order FIC form of the mass balance equation [37,38] has been adopted as the basis for deriving the stabilized FIC-FEM formulation for 2D and 3D poromechanics problems. This new formulation takes advantage of the second order derivatives terms to provide a simpler procedure for obtaining a residual-based stable form of the mass balance equation suitable for finite element analysis. Proof of the good results of this formulation regarding convergence and mass conservation are given in [38].

This chapter is organized as follows. First, the finite element formulation is introduced, starting from the governing equations of the problem and the fully coupled system of equations to be solved. After that, the instability problem is stated and the stabilization by means of the second-order FIC form of the mass balance equation in space is thoroughly detailed. Finally, two examples are solved in order to test the implemented FIC-stabilized element and compare it to stable elements with a higher order interpolation for the displacement field.

## 2.2 Finite element formulation

 ${\displaystyle \sigma _{ij}'=\sigma _{ij}+\alpha p\delta _{ij}}$
(2.1)

where ${\textstyle \delta _{ij}}$ is the Kronecker delta and ${\textstyle \alpha }$ is Biot's coefficient [103]:

 ${\displaystyle \alpha =1-{\frac {K}{K_{s}}}\leq 1}$
(2.2)

with ${\textstyle K_{s}}$ being the bulk modulus of the solid phase and ${\textstyle K}$ the bulk modulus of the porous medium:

 ${\displaystyle K={\frac {E}{3\left(1-2\nu \right)}}}$
(2.3)

where ${\displaystyle E}$ is the Young's modulus and ${\textstyle \nu }$ is the Poisson's ratio.

In order to account for the coupling between the solid and fluid phases, the behaviour of a saturated porous medium is governed by the combination of two equations: the balance of momentum for the mixture solid-fluid, and the mass balance for the pore fluid.

Balance of momentum for the mixture solid-fluid

 ${\displaystyle {\frac {\partial \sigma _{ij}}{\partial x_{j}}}+\rho b_{i}=\rho {\ddot {u}}_{i}}$
(2.4)

where ${\textstyle b_{i}}$ is the body force per unit mass, ${\textstyle {\ddot {u}}_{i}}$ is the acceleration of the solid skeleton and ${\textstyle \rho }$ is the density of the solid-fluid mixture ${\textstyle \rho =\phi \rho _{f}+(1-\phi )\rho _{s}}$, being ${\textstyle \phi }$ the porosity of the soil, ${\textstyle \rho _{f}}$ the density of the fluid and ${\textstyle \rho _{s}}$ the density of the solid.

Let us introduce the general definition for the effective stress (2.1) into the balance of momentum equation (2.4). This gives

 ${\displaystyle {\frac {\partial }{\partial x_{j}}}\left(\sigma _{ij}'-\alpha p\delta _{ij}\right)+\rho b_{i}-\rho {\ddot {u}}_{i}=0}$
(2.5)

Mass balance for the pore fluid

 ${\displaystyle {\frac {\partial \zeta }{\partial t}}+{\frac {\partial q_{i}}{\partial x_{i}}}=0}$
(2.6)

where ${\textstyle \zeta }$ is the fluid mass content per unit volume and ${\textstyle q_{i}}$ represents the flow rate. It should be noted that Equation (2.6) is in a linearized form of the general mass balance expression as the fluid density variation effect has been ignored.

The mass balance equation is modified as follows. Let us first consider the constitutive equation presented in [2] relating the fluid mass content per unit volume ${\textstyle \zeta }$ with the volumetric strain of the solid skeleton ${\textstyle \epsilon }$ and the pore pressure ${\textstyle p}$

 ${\displaystyle \zeta =\alpha \epsilon +{\frac {p}{Q}}}$
(2.7)

where ${\textstyle Q}$ is a combined compressibility of the fluid-solid phases, also called Biot's modulus [103]

 ${\displaystyle {\frac {1}{Q}}={\frac {\alpha {-\phi }}{K_{s}}}+{\frac {\phi }{K_{f}}}}$
(2.8)

being ${\textstyle K_{f}}$ the bulk modulus of the fluid phase.

Let us now make use of the generalized form of Darcy's law

 ${\displaystyle q_{i}=-{\frac {1}{\mu }}k_{ij}\left({\frac {\partial p}{\partial x_{j}}}-\rho _{f}b_{j}\right)}$
(2.9)

where ${\textstyle \mu }$ is the dynamic viscosity of the fluid and ${\textstyle k_{ij}}$ is the intrinsic permeability matrix of the porous medium.

Taking into account the relations from equations (2.7) and (2.9), the mass balance equation (2.6) can be rewritten as

 ${\displaystyle \alpha {\dot {\epsilon }}+{\frac {\dot {p}}{Q}}+{\frac {\partial }{\partial x_{i}}}\left[-{\frac {1}{\mu }}k_{ij}\left({\frac {\partial p}{\partial x_{j}}}-\rho _{f}b_{j}\right)\right]=0}$
(2.10)

Equations (2.5) and (2.10) have to be supplemented by a constitutive law for the solid skeleton. In general, for any non-linear material we may consider an incremental form with a tangent constitutive tensor ${\textstyle D_{ijkl}}$ as

 ${\displaystyle d\sigma _{ij}'=D_{ijkl}\,d\varepsilon _{kl}}$
(2.11)

with

 ${\displaystyle \varepsilon _{ij}={\frac {1}{2}}\left({\frac {\partial u_{i}}{\partial x_{j}}}+{\frac {\partial u_{j}}{\partial x_{i}}}\right)}$
(2.12)

The boundary conditions for this problem are specified as:

 ${\displaystyle u_{i}={\tilde {u}}_{i}\qquad {\hbox{ on}}\quad \Gamma _{u}}$ (2.13) ${\displaystyle \left(\sigma _{ij}'-\alpha p\delta _{ij}\right)n_{j}={\tilde {t}}_{i}\qquad {\hbox{ on}}\quad \Gamma _{t}}$ (2.14)

where ${\textstyle {\tilde {u}}_{i}}$ and ${\textstyle {\tilde {t}}_{i}}$ are the prescribed displacements and surface tractions, respectively.

 ${\displaystyle p={\tilde {p}}\qquad {\hbox{ on}}\quad \Gamma _{p}}$ (2.15) ${\displaystyle \left[-{\frac {1}{\mu }}k_{ij}\left({\frac {\partial p}{\partial x_{j}}}-\rho _{f}b_{j}\right)\right]n_{i}={\tilde {q}}_{n}\qquad {\hbox{ on}}\quad \Gamma _{q}}$ (2.16)

where ${\textstyle {\tilde {p}}}$ and ${\textstyle {\tilde {q}}_{n}}$ are the prescribed pore pressure and the normal flow rate, respectively.

In order to express the resultant system of equations in a more compact manner, we rewrite equations (2.5) and (2.10) in matrix form using Voigt notation as follows.

• Balance of momentum

 ${\displaystyle {\boldsymbol {S}}^{T}\left({\boldsymbol {\sigma }}'-\alpha p{\boldsymbol {m}}\right)+\rho {\boldsymbol {b}}-\rho {\ddot {\boldsymbol {u}}}={\boldsymbol {0}}}$
(2.17)

• Mass balance

 ${\displaystyle \alpha {\boldsymbol {m}}^{T}{\boldsymbol {S}}{\dot {\boldsymbol {u}}}+{\frac {\dot {p}}{Q}}+{\boldsymbol {\nabla }}^{T}\left[-{\frac {1}{\mu }}{\boldsymbol {k}}\left({\boldsymbol {\nabla }}p-\rho _{f}{\boldsymbol {b}}\right)\right]=0}$
(2.18)

For a general 3D case, we have

 ${\displaystyle {\boldsymbol {S}}={\begin{bmatrix}{\frac {\partial }{\partial x}}&0&0\\0&{\frac {\partial }{\partial y}}&0\\0&0&{\frac {\partial }{\partial z}}\\{\frac {\partial }{\partial y}}&{\frac {\partial }{\partial x}}&0\\0&{\frac {\partial }{\partial z}}&{\frac {\partial }{\partial y}}\\{\frac {\partial }{\partial z}}&0&{\frac {\partial }{\partial x}}\\\end{bmatrix}}\quad {\hbox{;}}\quad {\boldsymbol {\nabla }}={\begin{bmatrix}{\frac {\partial }{\partial x}},{\frac {\partial }{\partial y}},{\frac {\partial }{\partial z}}\end{bmatrix}}^{T}}$
(2.19)
 ${\displaystyle {\boldsymbol {\sigma }}'={\begin{bmatrix}\sigma _{xx}',&\sigma _{yy}',&\sigma _{zz}',&\sigma _{xy}',&\sigma _{yz}',&\sigma _{zx}'\end{bmatrix}}^{T}}$
(2.20)
 ${\displaystyle {\boldsymbol {m}}={\begin{bmatrix}1,&1,&1,&0,&0,&0\end{bmatrix}}^{T}}$
(2.21)

The above system of partial differential equations (2.17) and (2.18) can be discretized by interpolating the displacement and pressure fields as: ${\textstyle {\boldsymbol {u}}={\boldsymbol {N}}_{u}{\bar {\boldsymbol {u}}}}$ and ${\textstyle p={\boldsymbol {N}}_{p}{\bar {\boldsymbol {p}}}}$ where ${\textstyle {\bar {(\cdot )}}}$ denotes nodal values and

 ${\displaystyle {\boldsymbol {N}}_{u}={\begin{bmatrix}{\hat {N}}_{1}&0&0&{\hat {N}}_{2}&0&0&...&{\hat {N}}_{\hat {n}}&0&0\\0&{\hat {N}}_{1}&0&0&{\hat {N}}_{2}&0&...&0&{\hat {N}}_{\hat {n}}&0\\0&0&{\hat {N}}_{1}&0&0&{\hat {N}}_{2}&...&0&0&{\hat {N}}_{\hat {n}}\end{bmatrix}}}$
(2.22)

 ${\displaystyle {\boldsymbol {N}}_{p}={\begin{bmatrix}N_{1},N_{2},...,N_{n}\end{bmatrix}}}$
(2.23)

with ${\textstyle {\hat {N}}_{i}}$ and ${\textstyle N_{i}}$ being, respectively, the shape functions of the displacement and pressure interpolations, which do not necessarily need to coincide.

Following the standard Galerkin technique [104,105], we left multiply Equation (2.17) by test functions ${\textstyle {\boldsymbol {N}}_{u}^{T}}$ and Equation (2.18) by test functions ${\textstyle {\boldsymbol {N}}_{p}^{T}}$ and integrate them over the spatial domain ${\textstyle \Omega }$ to obtain the following set of ordinary differential equations

• Balance of momentum

 ${\displaystyle {\boldsymbol {r}}_{u}:={\boldsymbol {M}}{\ddot {\bar {\boldsymbol {u}}}}+\int _{\Omega }{\boldsymbol {B}}^{T}{\boldsymbol {\sigma }}'\,d\Omega {-}{\boldsymbol {Q}}{\bar {\boldsymbol {p}}}-{\boldsymbol {f}}_{u}={\boldsymbol {0}}}$
(2.24)

• Mass balance
 ${\displaystyle {\boldsymbol {r}}_{p}:={\boldsymbol {Q}}^{T}{\dot {\bar {\boldsymbol {u}}}}+{\boldsymbol {C}}{\dot {\bar {\boldsymbol {p}}}}+{\boldsymbol {H}}{\bar {\boldsymbol {p}}}-{\boldsymbol {f}}_{p}={\boldsymbol {0}}}$
(2.25)

where ${\textstyle {\boldsymbol {r}}_{u}}$ and ${\textstyle {\boldsymbol {r}}_{p}}$ are the residual vectors for the momentum and the mass balance equation, and

 ${\displaystyle {\boldsymbol {M}}=\int _{\Omega }{\boldsymbol {N}}_{u}^{T}\rho {\boldsymbol {N}}_{u}\,d\Omega \quad {\hbox{;}}\quad {\boldsymbol {Q}}=\int _{\Omega }{\boldsymbol {B}}^{T}\alpha {\boldsymbol {m}}{\boldsymbol {N}}_{p}\,d\Omega }$
(2.26)

 ${\displaystyle {\boldsymbol {C}}=\int _{\Omega }{\boldsymbol {N}}_{p}^{T}{\frac {1}{Q}}{\boldsymbol {N}}_{p}\,d\Omega \quad {\hbox{;}}\quad {\boldsymbol {H}}=\int _{\Omega }\left({\boldsymbol {\nabla }}{\boldsymbol {N}}_{p}\right)^{T}{\frac {1}{\mu }}{\boldsymbol {k}}{\boldsymbol {\nabla }}{\boldsymbol {N}}_{p}\,d\Omega }$
(2.27)
 ${\displaystyle {\boldsymbol {f}}_{u}=\int _{\Omega }{\boldsymbol {N}}_{u}^{T}\rho {\boldsymbol {b}}\,d\Omega {+}\int _{\Gamma _{t}}{\boldsymbol {N}}_{u}^{T}{\tilde {\boldsymbol {t}}}\,d\Gamma }$
(2.28)
 ${\displaystyle {\boldsymbol {f}}_{p}=\int _{\Omega }\left({\boldsymbol {\nabla }}{\boldsymbol {N}}_{p}\right)^{T}{\frac {1}{\mu }}{\boldsymbol {k}}\rho _{f}{\boldsymbol {b}}\,d\Omega {-\int }_{\Gamma _{q}}{\boldsymbol {N}}_{p}^{T}{\tilde {q}}_{n}\,d\Gamma }$
(2.29)

with ${\textstyle {\boldsymbol {B}}={\boldsymbol {S}}{\boldsymbol {N}}_{u}}$.

The time derivatives of ${\textstyle {\boldsymbol {u}}}$ and ${\textstyle p}$ are approximated using the Generalized Newmark scheme. Thus, for a new time step ${\textstyle n+1}$, we use the GN22 scheme for displacements [106]:

 ${\displaystyle {\boldsymbol {u}}_{n+1}={\boldsymbol {u}}_{n}+\Delta t{\dot {\boldsymbol {u}}}_{n}+\left({\frac {1}{2}}-\beta \right)\Delta t^{2}{\ddot {\boldsymbol {u}}}_{n}+\beta \Delta t^{2}{\ddot {\boldsymbol {u}}}_{n+1}}$
(2.30.a)
 ${\displaystyle {\dot {\boldsymbol {u}}}_{n+1}={\dot {\boldsymbol {u}}}_{n}+(1-\gamma )\Delta t{\ddot {\boldsymbol {u}}}_{n}+\gamma \Delta t{\ddot {\boldsymbol {u}}}_{n+1}}$
(2.30.b)

and the GN11 scheme for pressure:

 ${\displaystyle p_{n+1}=p_{n}+\left(1-\theta \right)\Delta t{\dot {p}}_{n}+\theta \Delta t{\dot {p}}_{n+1}}$
(2.31)

which are unconditionally stable for ${\textstyle \beta \geq {\frac {1}{4}}}$, ${\textstyle \gamma \geq {\frac {1}{2}}}$ and ${\textstyle \theta \geq {\frac {1}{2}}}$ [107].

The residual vectors for the momentum and mass balance equations can be rewritten as

 ${\displaystyle \left({\boldsymbol {r}}_{u}\right)_{n+1}={\boldsymbol {M}}\left[{\frac {{\bar {\boldsymbol {u}}}_{n+1}-{\hat {\bar {\boldsymbol {u}}}}_{n}}{\beta \Delta t^{2}}}\right]+\left(\int _{\Omega }{\boldsymbol {B}}^{T}{\boldsymbol {\sigma }}'\,d\Omega \right)_{n+1}-{\boldsymbol {Q}}{\bar {\boldsymbol {p}}}_{n+1}-\left({\boldsymbol {f}}_{u}\right)_{n+1}}$
(2.32)

 ${\displaystyle \left({\boldsymbol {r}}_{p}\right)_{n+1}={\boldsymbol {Q}}^{T}\left[{\frac {\gamma {\bar {\boldsymbol {u}}}_{n+1}}{\beta \Delta t}}+{\hat {\dot {\bar {\boldsymbol {u}}}}}_{n}\right]+{\boldsymbol {C}}\left[{\frac {{\bar {\boldsymbol {p}}}_{n+1}}{\theta \Delta t}}+{\hat {\dot {\bar {\boldsymbol {p}}}}}_{n}\right]+{\boldsymbol {H}}{\bar {\boldsymbol {p}}}_{n+1}-\left({\boldsymbol {f}}_{p}\right)_{n+1}}$
(2.33)

where ${\textstyle {\bar {\boldsymbol {u}}}_{n+1}}$ and ${\textstyle {\bar {\boldsymbol {p}}}_{n+1}}$ are the nodal unknowns at time ${\textstyle n+1}$, and ${\textstyle {\hat {\bar {\boldsymbol {u}}}}_{n}}$, ${\textstyle {\hat {\dot {\bar {\boldsymbol {u}}}}}_{n}}$ and ${\textstyle {\hat {\dot {\bar {\boldsymbol {p}}}}}_{n}}$ stand for values that are computed from known parameters at time ${\textstyle n}$ as:

 ${\displaystyle {\hat {\bar {\boldsymbol {u}}}}_{n}={\bar {\boldsymbol {u}}}_{n}+\Delta t{\dot {\bar {\boldsymbol {u}}}}_{n}+\left({\frac {1}{2}}-\beta \right)\Delta t^{2}{\ddot {\bar {\boldsymbol {u}}}}_{n}}$
(2.34)
 ${\displaystyle {\hat {\dot {\bar {\boldsymbol {u}}}}}_{n}={\dot {\bar {\boldsymbol {u}}}}_{n}+\left(1-\gamma \right)\Delta t{\ddot {\bar {\boldsymbol {u}}}}_{n}-{\frac {\gamma }{\beta \Delta t}}{\hat {\bar {\boldsymbol {u}}}}_{n}}$
(2.35)
 ${\displaystyle {\hat {\dot {\bar {\boldsymbol {p}}}}}_{n}={\frac {\theta {-1}}{\theta }}{\dot {\bar {\boldsymbol {p}}}}_{n}-{\frac {1}{\theta \Delta t}}{\bar {\boldsymbol {p}}}_{n}}$
(2.36)

The Newton-Raphson method is used to solve the problem iteratively. Thus, at a time step ${\textstyle n+1}$ and the iteration ${\textstyle i+1}$ we have:

 ${\displaystyle {\begin{bmatrix}\displaystyle {\frac {\partial {\boldsymbol {r}}_{u}}{\partial {\bar {\boldsymbol {u}}}}}&\displaystyle {\frac {\partial {\boldsymbol {r}}_{u}}{\partial {\bar {\boldsymbol {p}}}}}\\\displaystyle {\frac {\partial {\boldsymbol {r}}_{p}}{\partial {\bar {\boldsymbol {u}}}}}&\displaystyle {\frac {\partial {\boldsymbol {r}}_{p}}{\partial {\bar {\boldsymbol {p}}}}}\end{bmatrix}}_{n+1}^{i}{\begin{bmatrix}\delta {\bar {\boldsymbol {u}}}\\\delta {\bar {\boldsymbol {p}}}\end{bmatrix}}_{n+1}^{i+1}=-{\begin{bmatrix}{\boldsymbol {r}}_{u}\\{\boldsymbol {r}}_{p}\end{bmatrix}}_{n+1}^{i}}$
(2.37)

 ${\displaystyle {\begin{bmatrix}\displaystyle {\frac {1}{\beta \Delta t^{2}}}{\boldsymbol {M}}+{\boldsymbol {K}}&-{\boldsymbol {Q}}\\\displaystyle {\frac {\gamma }{\beta \Delta t}}{\boldsymbol {Q}}^{T}&\displaystyle {\frac {1}{\theta \Delta t}}{\boldsymbol {C}}+{\boldsymbol {H}}\end{bmatrix}}_{n+1}^{i}{\begin{bmatrix}\delta {\bar {\boldsymbol {u}}}\\\delta {\bar {\boldsymbol {p}}}\end{bmatrix}}_{n+1}^{i+1}=-{\begin{bmatrix}{\boldsymbol {r}}_{u}\\{\boldsymbol {r}}_{p}\end{bmatrix}}_{n+1}^{i}}$
(2.38)

where ${\textstyle {\boldsymbol {K}}}$ is the tangent stiffness matrix:

 ${\displaystyle {\boldsymbol {K}}={\frac {\partial }{\partial {\bar {\boldsymbol {u}}}}}\int _{\Omega }{\boldsymbol {B}}^{T}{\boldsymbol {\sigma }}'\,d\Omega =\int _{\Omega }{\boldsymbol {B}}^{T}{\boldsymbol {D}}{\boldsymbol {B}}\,d\Omega }$
(2.39)

Note that the system of equations (2.38) can be made symmetric by multiplying Equation (2.33) by ${\textstyle -\beta \Delta t/\gamma }$.

## 2.3 Stabilization using Finite Increment Calculus (FIC)

### 2.3.1 Introduction to the FIC-stabilization procedure

For a quasi-static problem the term involving the mass matrix is omitted. Also, in the undrained-incompressible limit, i.e. when ${\textstyle {\boldsymbol {k}}\to {\boldsymbol {0}}}$ and ${\textstyle Q\to \infty }$, the matrices ${\textstyle {\boldsymbol {C}}}$ and ${\textstyle {\boldsymbol {H}}}$ vanish and the system of equations (2.38) becomes

 ${\displaystyle {\begin{bmatrix}{\boldsymbol {K}}&-{\boldsymbol {Q}}\\-{\boldsymbol {Q}}^{T}&{\boldsymbol {0}}\end{bmatrix}}_{n+1}^{i}{\begin{bmatrix}\delta {\bar {\boldsymbol {u}}}\\\delta {\bar {\boldsymbol {p}}}\end{bmatrix}}_{n+1}^{i+1}=-{\begin{bmatrix}{\boldsymbol {r}}_{u}\\{\boldsymbol {r}}_{p}^{*}\end{bmatrix}}_{n+1}^{i}}$
(2.40)

The resultant system matrix is almost identical to that frequently encountered in the solution of incompressible elasticity or incompressible fluid mechanics problems [36,26]. In such cases, the spaces used to approximate the displacement (or the velocity) and pressure fields must fulfill the Babuska-Brezzi conditions [11,12] or the Zienkiewicz-Taylor patch test [108,109] in order to avoid spurious oscillations and locking in the pressure field.

With the equations presented so far, these conditions can be accomplished using shape functions for the displacement field one order higher than those for the pressure field. Some examples of “stable” 2D elements are depicted in Figure 5.

 (a) Quadratic displacement and linear pressure. (b) Biquadratic displacement and bilinear pressure. Figure 5: Some 2D elements that fulfil the Babuska-Brezzi conditions.

However, low order finite elements with equal order interpolation for the displacements and the pressure are very attractive due to their simplicity and efficiency, and so stabilization techniques must be applied if one aims at solving large scale 3D computations.

The stabilization approach implemented in this work relies on the Finite Increment Calculus (FIC) method [33,34,35,36,37,38] and affects only the continuity equation with the balance of momentum remaining unchanged.

For the sake of clarity, let us redefine the mass balance equation (2.10) as:

 ${\displaystyle r_{p}:=\alpha {\dot {\epsilon }}+{\frac {\dot {p}}{Q}}+{\frac {\partial }{\partial x_{j}}}\left[-{\frac {1}{\mu }}k_{jk}\left({\frac {\partial p}{\partial x_{k}}}-\rho _{f}b_{k}\right)\right]=0}$
(2.41)

As stated in the introduction, this technique is based on the second-order FIC form of the mass balance equation in space for a quasi-incompressible fluid. This can be written as [38]:

 ${\displaystyle r_{p}+{\frac {h_{i}^{2}}{12}}{\frac {\partial ^{2}r_{p}}{\partial x_{i}^{2}}}=0\qquad {\hbox{in }}\Omega \qquad i=1,n_{s}}$
(2.42)

where ${\textstyle n_{s}}$ are the number of space dimensions (i.e. ${\textstyle n_{s}=3}$ for 3D problems).

In a 2D case Equation (2.42) is obtained by expressing the balance of mass in a rectangular domain of finite size with dimensions ${\textstyle h_{1}\times h_{2}}$, where ${\textstyle h_{i}}$ are arbitrary distances, and retaining up to third-order terms in the Taylor series expansions used for expressing the change of mass within the balance domain [37]. The derivation of Equation (2.42) for a 1D problem is shown in [38].

The FIC term in Equation (2.42) plays the role of a space stabilization term, in which ${\textstyle h_{i}}$ are space dimensions related to the characteristic element dimensions. Note that for ${\textstyle h_{i}\to 0}$ the standard form of the mass balance equation (2.41), as given by the infinitesimal theory, is recovered.

Equation (2.42) can be interpreted as a particular class of residual-based stabilization methods applied to the strong form of the mass conservation equation. This ensures the consistency of the stabilization method in the discretized problem.

### 2.3.2 Modified FIC-stabilized form of the mass balance equation

Let us first expand Equation (2.42) using Equation (2.41) as

 ${\displaystyle \alpha {\dot {\epsilon }}+{\frac {1}{Q}}{\dot {p}}+{\frac {\partial }{\partial x_{j}}}\left[-{\frac {1}{\mu }}k_{jk}\left({\frac {\partial p}{\partial x_{k}}}-\rho _{f}b_{k}\right)\right]+{\frac {h_{i}^{2}}{12}}{\frac {\partial ^{2}}{\partial x_{i}^{2}}}\left(\alpha {\dot {\epsilon }}\right)}$ ${\displaystyle +{\frac {h_{i}^{2}}{12}}{\frac {\partial ^{2}}{\partial x_{i}^{2}}}\left({\frac {1}{Q}}{\dot {p}}\right)+{\frac {h_{i}^{2}}{12}}{\frac {\partial ^{2}}{\partial x_{i}^{2}}}\left\{{\frac {\partial }{\partial x_{j}}}\left[-{\frac {1}{\mu }}k_{jk}\left({\frac {\partial p}{\partial x_{k}}}-\rho _{f}b_{k}\right)\right]\right\}=0}$
(2.43)

Next we split the stress tensor into its deviatoric and volumetric components as

 ${\displaystyle \sigma _{ij}=s_{ij}+\sigma \delta _{ij}}$
(2.44)

where ${\textstyle s_{ij}}$ is the deviatoric stress tensor and ${\textstyle \sigma }$ is the hydrostatic pressure defined as ${\textstyle \sigma =\sigma _{ii}/3}$.

In the same way, we split the strain tensor into its deviatoric and volumetric parts as

 ${\displaystyle \varepsilon _{ij}=e_{ij}+{\frac {1}{3}}\epsilon \,\delta _{ij}}$
(2.45)

with ${\textstyle e_{ij}}$ being the deviatoric strain tensor and ${\textstyle \epsilon =\varepsilon _{ii}}$ the volumetric dilation.

Now let us write the isotropic linear elastic constitutive equations:

 ${\displaystyle s_{ij}=2Ge_{ij}}$
(2.46)

where ${\textstyle G}$ is the shear modulus,

 ${\displaystyle \sigma =\sigma '{-\alpha }p}$
(2.47)

where ${\textstyle \sigma '=\sigma _{ii}'/3}$ is the mean effective stress and ${\textstyle \alpha }$ is the Biot's coefficient defined in Equation (2.2).

Combining Equations (2.44), (2.45), (2.46) and (2.47), we express the stress tensor as

 ${\displaystyle \sigma _{ij}=2G\varepsilon _{ij}-{\frac {2G}{3}}\epsilon \delta _{ij}+\sigma '\delta _{ij}-\alpha p\delta _{ij}}$
(2.48)

At this point, we substitute the above expression in the balance of momentum (2.4). This gives

 ${\displaystyle {\frac {\partial }{\partial x_{j}}}\left(2G\varepsilon _{ij}\right)-{\frac {\partial }{\partial x_{i}}}\left({\frac {2G}{3}}\epsilon \right)+{\frac {\partial \sigma '}{\partial x_{i}}}-{\frac {\partial }{\partial x_{i}}}\left(\alpha p\right)+\rho b_{i}-\rho {\ddot {u}}_{i}=0}$
(2.49)

From this point forward, in the derivation of equations (2.43), (2.49) and in the following, we neglect the space and time changes of ${\textstyle \alpha }$, ${\textstyle Q}$, ${\textstyle G}$ and ${\textstyle \rho }$ in the derivatives.

If we now apply the divergence operator to both sides of Equation (2.49), differentiate it with respect to time, and rearrange the terms we can obtain

 ${\displaystyle {\frac {\partial ^{2}{\dot {\epsilon }}}{\partial x_{i}^{2}}}={\frac {3}{2G}}\left[2G{\frac {\partial ^{2}{\dot {\varepsilon }}_{ij}}{\partial x_{i}\partial x_{j}}}+{\frac {\partial ^{2}{\dot {\sigma }}'}{\partial x_{i}^{2}}}-\alpha {\frac {\partial ^{2}{\dot {p}}}{\partial x_{i}^{2}}}+\rho {\frac {\partial {\dot {b}}_{i}}{\partial x_{i}}}-\rho {\frac {D}{Dt}}\left({\frac {\partial {\ddot {u}}_{i}}{\partial x_{i}}}\right)\right]}$
(2.50)

In the previous equation the term ${\textstyle \rho {\frac {\partial {\dot {b}}_{i}}{\partial x_{i}}}}$ can be neglected and the term involving the third derivative of the displacements with respect to time can be obtained from the mass balance equation. Thus if we differentiate Equation (2.10) twice with respect to time assuming that the time derivative of the body force per unit mass is negligible, and we take into account the definition of the volumetric strain ${\textstyle \epsilon ={\frac {\partial u_{i}}{\partial x_{i}}}}$, we obtain

 ${\displaystyle {\frac {D}{Dt}}\left({\frac {\partial {\ddot {u}}_{i}}{\partial x_{i}}}\right)={\frac {1}{\alpha }}\left[{\frac {\partial }{\partial x_{j}}}\left({\frac {1}{\mu }}k_{jk}{\frac {\partial {\ddot {p}}}{\partial x_{k}}}\right)-{\frac {1}{Q}}{\frac {D{\ddot {p}}}{Dt}}\right]}$
(2.51)

Introducing (2.51) into (2.50) gives

 ${\displaystyle {\frac {\partial ^{2}{\dot {\epsilon }}}{\partial x_{i}^{2}}}={\frac {3}{2G}}\left\{2G{\frac {\partial ^{2}{\dot {\varepsilon }}_{ij}}{\partial x_{i}\partial x_{j}}}+{\frac {\partial ^{2}{\dot {\sigma }}'}{\partial x_{i}^{2}}}-\alpha {\frac {\partial ^{2}{\dot {p}}}{\partial x_{i}^{2}}}\right.}$ ${\displaystyle \left.-{\frac {\rho }{\alpha }}\left[{\frac {\partial }{\partial x_{j}}}\left({\frac {1}{\mu }}k_{jk}{\frac {\partial {\ddot {p}}}{\partial x_{k}}}\right)-{\frac {1}{Q}}{\frac {D{\ddot {p}}}{Dt}}\right]\right\}}$
(2.52)

At this point, we can already substitute the above relation in Equation (2.43) giving

 ${\displaystyle \alpha {\dot {\epsilon }}+{\frac {1}{Q}}{\dot {p}}+{\frac {\partial }{\partial x_{j}}}\left[-{\frac {1}{\mu }}k_{jk}\left({\frac {\partial p}{\partial x_{k}}}-\rho _{f}b_{k}\right)\right]+{\frac {h_{i}^{2}\alpha }{8G}}\left\{2G{\frac {\partial ^{2}{\dot {\varepsilon }}_{ij}}{\partial x_{i}\partial x_{j}}}\right.}$ ${\displaystyle \left.+{\frac {\partial ^{2}{\dot {\sigma }}'}{\partial x_{i}^{2}}}-\alpha {\frac {\partial ^{2}{\dot {p}}}{\partial x_{i}^{2}}}-{\frac {\rho }{\alpha }}\left[{\frac {\partial }{\partial x_{j}}}\left({\frac {1}{\mu }}k_{jk}{\frac {\partial {\ddot {p}}}{\partial x_{k}}}\right)-{\frac {1}{Q}}{\frac {D{\ddot {p}}}{Dt}}\right]\right\}+{\frac {h_{i}^{2}}{12}}{\frac {\partial ^{2}}{\partial x_{i}^{2}}}\left({\frac {1}{Q}}{\dot {p}}\right)}$ ${\displaystyle +{\frac {h_{i}^{2}}{12}}{\frac {\partial ^{2}}{\partial x_{i}^{2}}}\left\{{\frac {\partial }{\partial x_{j}}}\left[-{\frac {1}{\mu }}k_{jk}\left({\frac {\partial p}{\partial x_{k}}}-\rho _{f}b_{k}\right)\right]\right\}=0}$
(2.53)

The above equation can be simplified if we take into account that the last term involves the fourth order spatial derivative of the pore pressure and the third order spatial derivative of the body force per unit mass. In practical problems these terms are either zero (for constant body forces) or much smaller than the others. Hence, these terms will be omitted for the rest of this work. Furthermore, numerical results have shown that the terms ${\textstyle {\frac {\partial }{\partial x_{j}}}\left({\frac {1}{\mu }}k_{jk}{\frac {\partial {\ddot {p}}}{\partial x_{k}}}\right)}$ and ${\textstyle {\frac {1}{Q}}{\frac {D{\ddot {p}}}{Dt}}}$ can be neglected without loss of accuracy. Thereby, Equation (2.53) is written in the simpler form as

 ${\displaystyle \alpha {\dot {\epsilon }}+{\frac {1}{Q}}{\dot {p}}+{\frac {\partial }{\partial x_{j}}}\left[-{\frac {1}{\mu }}k_{jk}\left({\frac {\partial p}{\partial x_{k}}}-\rho _{f}b_{k}\right)\right]}$ ${\displaystyle +{\frac {h_{i}^{2}\alpha }{8G}}\left\{2G{\frac {\partial ^{2}{\dot {\varepsilon }}_{ij}}{\partial x_{i}\partial x_{j}}}+{\frac {\partial ^{2}{\dot {\sigma }}'}{\partial x_{i}^{2}}}-\alpha {\frac {\partial ^{2}{\dot {p}}}{\partial x_{i}^{2}}}\right\}+{\frac {h_{i}^{2}}{12}}{\frac {\partial ^{2}}{\partial x_{i}^{2}}}\left({\frac {1}{Q}}{\dot {p}}\right)=0}$
(2.54)

In the following we will assume ${\textstyle h_{i}=h}$ where ${\textstyle h}$ is a characteristic length related to a typical average dimension of each element in the mesh. After rearranging terms, Equation (2.54) can be rewritten as

 ${\displaystyle \alpha {\dot {\epsilon }}+{\frac {1}{Q}}{\dot {p}}+{\frac {\partial }{\partial x_{i}}}\left[-{\frac {1}{\mu }}k_{ij}\left({\frac {\partial p}{\partial x_{j}}}-\rho _{f}b_{j}\right)\right]}$ ${\displaystyle +\tau {\frac {\partial }{\partial x_{i}}}\left[2G{\frac {\partial {\dot {\varepsilon }}_{ij}}{\partial x_{j}}}+{\frac {\partial {\dot {\sigma }}'}{\partial x_{i}}}-\left(\alpha {-}{\frac {2G}{3\alpha Q}}\right){\frac {\partial {\dot {p}}}{\partial x_{i}}}\right]=0}$
(2.55)

where ${\textstyle \tau }$ is a stabilization parameter given by

 ${\displaystyle \tau ={\frac {h^{2}\alpha }{8G}}}$
(2.56)

The form of the stabilization parameter in Equation (2.56) resembles that typically used in other stabilized procedures. We note that this term has naturally formed from the FIC formulation.

In the examples solved in this work we have chosen ${\textstyle h=l^{e}}$, where ${\textstyle l^{e}}$ is a characteristic element length that, for 2D problems, is taken as the diameter of a circle with the area of the element, while for 3D problems, it is the diameter of a sphere with the volume of the element. In essence

 ${\displaystyle l^{e}={\sqrt {\frac {4A^{e}}{\pi }}}\quad {\hbox{in 2D}}}$
(2.57)

 ${\displaystyle l^{e}={\sqrt {=}}3\quad {\hbox{in 3D}}}$
(2.58)

where ${\textstyle A^{e}}$ and ${\textstyle V^{e}}$ represent the area and the volume of the element, respectively.

The presence of the characteristic element length ${\textstyle l^{e}}$ in Equation (2.56) helps us control the diffusion of the stabilized solution, avoiding over-diffusive numerical results provided that fine enough discretizations are used. This fact is clearly shown in the example of Section 2.4.2.

### 2.3.3 Variational form of the FIC-stabilized mass balance equation

The variational expression of the FIC-stabilized mass balance equation is obtained by multiplying Equation (2.55) by arbitrary test functions ${\textstyle v}$ and integrating over the domain ${\textstyle \Omega }$ to give

 ${\displaystyle \int _{\Omega }v\,\alpha {\dot {\epsilon }}\,d\Omega {+}\int _{\Omega }v\,{\frac {1}{Q}}{\dot {p}}\,d\Omega {+}\int _{\Omega }v\,{\frac {\partial }{\partial x_{i}}}\left[-{\frac {1}{\mu }}k_{ij}\left({\frac {\partial p}{\partial x_{j}}}-\rho _{f}b_{j}\right)\right]\,d\Omega }$ ${\displaystyle +\int _{\Omega }v\,\tau {\frac {\partial }{\partial x_{i}}}\left[2G{\frac {\partial {\dot {\varepsilon }}_{ij}}{\partial x_{j}}}+{\frac {\partial {\dot {\sigma }}'}{\partial x_{i}}}-\left(\alpha {-}{\frac {2G}{3\alpha Q}}\right){\frac {\partial {\dot {p}}}{\partial x_{i}}}\right]\,d\Omega =0}$
(2.59)

Integrating by parts the last two integrals of Equation (2.59) and applying the divergence theorem yields

 ${\displaystyle \int _{\Omega }v\,\alpha {\dot {\epsilon }}\,d\Omega {+}\int _{\Omega }v\,{\frac {1}{Q}}{\dot {p}}\,d\Omega -\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}\left[-{\frac {1}{\mu }}k_{ij}\left({\frac {\partial p}{\partial x_{j}}}-\rho _{f}b_{j}\right)\right]\,d\Omega }$ ${\displaystyle +\int _{\Gamma }v\left[-{\frac {1}{\mu }}k_{ij}\left({\frac {\partial p}{\partial x_{j}}}-\rho _{f}b_{j}\right)\right]n_{i}\,d\Gamma }$ ${\displaystyle -\int _{\Omega }\tau {\frac {\partial v}{\partial x_{i}}}\left[2G{\frac {\partial {\dot {\varepsilon }}_{ij}}{\partial x_{j}}}+{\frac {\partial {\dot {\sigma }}'}{\partial x_{i}}}-\left(\alpha {-}{\frac {2G}{3\alpha Q}}\right){\frac {\partial {\dot {p}}}{\partial x_{i}}}\right]\,d\Omega }$ ${\displaystyle +\int _{\Gamma }v\,\tau \left[2G{\frac {\partial {\dot {\varepsilon }}_{ij}}{\partial x_{j}}}+{\frac {\partial {\dot {\sigma }}'}{\partial x_{i}}}-\left(\alpha {-}{\frac {2G}{3\alpha Q}}\right){\frac {\partial {\dot {p}}}{\partial x_{i}}}\right]n_{i}\,d\Gamma =0}$
(2.60)

where ${\textstyle n_{i}}$ are the components of the unit normal vector to the external boundary ${\textstyle \Gamma }$ of ${\textstyle \Omega }$.

Introducing the boundary condition (2.16) and rearranging terms we have

 ${\displaystyle \int _{\Omega }v\,\alpha {\dot {\epsilon }}\,d\Omega {+}\int _{\Omega }v\,{\frac {1}{Q}}{\dot {p}}\,d\Omega {+}\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}{\frac {1}{\mu }}k_{ij}{\frac {\partial p}{\partial x_{j}}}\,d\Omega {-}\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}{\frac {1}{\mu }}k_{ij}\rho _{f}b_{j}\,d\Omega }$ ${\displaystyle +\int _{\Gamma }v\,{\tilde {q}}_{n}\,d\Gamma {-}\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}\tau 2G{\frac {\partial {\dot {\varepsilon }}_{ij}}{\partial x_{j}}}\,d\Omega {-}\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}\tau {\frac {\partial {\dot {\sigma }}'}{\partial x_{i}}}\,d\Omega }$ ${\displaystyle +\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}\tau \left(\alpha {-}{\frac {2G}{3\alpha Q}}\right){\frac {\partial {\dot {p}}}{\partial x_{i}}}\,d\Omega {+}\int _{\Gamma }v\,\tau \left(2G{\frac {\partial {\dot {\varepsilon }}_{ij}}{\partial x_{j}}}n_{i}+{\frac {\partial {\dot {\sigma }}'}{\partial x_{i}}}n_{i}-\alpha {\frac {\partial {\dot {p}}}{\partial x_{i}}}n_{i}\right)\,d\Gamma }$ ${\displaystyle +\int _{\Gamma }v\,\tau {\frac {2G}{3\alpha Q}}{\frac {\partial {\dot {p}}}{\partial x_{i}}}n_{i}\,d\Gamma =0}$
(2.61)

Using the stress decomposition in (2.48), we can rewrite the boundary condition (2.14) as

 ${\displaystyle 2G\varepsilon _{ij}n_{j}-{\frac {2G}{3}}\epsilon n_{i}+\sigma 'n_{i}-\alpha pn_{i}={\tilde {t}}_{i}}$
(2.62)

Next if we apply the divergence operator to both sides of Equation (2.62), differentiate it with respect to time and rearrange terms we have

 ${\displaystyle 2G{\frac {\partial {\dot {\varepsilon }}_{ij}}{\partial x_{j}}}n_{i}+{\frac {\partial {\dot {\sigma }}'}{\partial x_{i}}}n_{i}-\alpha {\frac {\partial {\dot {p}}}{\partial x_{i}}}n_{i}={\frac {2G}{3}}{\frac {\partial {\dot {\epsilon }}}{\partial x_{i}}}n_{i}+{\frac {\partial {\dot {\tilde {t}}}_{i}}{\partial x_{i}}}}$
(2.63)

In all typical problems the divergence of the traction vector is zero and so it is the term ${\textstyle {\frac {\partial {\dot {\tilde {t}}}_{i}}{\partial x_{i}}}}$. Moreover, numerical tests have shown that the results are not affected by the term ${\textstyle {\frac {2G}{3}}{\frac {\partial {\dot {\epsilon }}}{\partial x_{i}}}n_{i}}$. Consequently, the stabilized mass balance equation is rewritten as

 ${\displaystyle \int _{\Omega }v\,\alpha {\dot {\epsilon }}\,d\Omega {+}\int _{\Omega }v\,{\frac {1}{Q}}{\dot {p}}\,d\Omega {+}\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}{\frac {1}{\mu }}k_{ij}{\frac {\partial p}{\partial x_{j}}}\,d\Omega }$ ${\displaystyle -\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}{\frac {1}{\mu }}k_{ij}\rho _{f}b_{j}\,d\Omega +\int _{\Gamma }v\,{\tilde {q}}_{n}\,d\Gamma {-}\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}\tau 2G{\frac {\partial {\dot {\varepsilon }}_{ij}}{\partial x_{j}}}\,d\Omega }$ ${\displaystyle -\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}\tau {\frac {\partial {\dot {\sigma }}'}{\partial x_{i}}}\,d\Omega +\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}\tau \left(\alpha {-}{\frac {2G}{3\alpha Q}}\right){\frac {\partial {\dot {p}}}{\partial x_{i}}}\,d\Omega }$ ${\displaystyle +\int _{\Gamma }v\,\tau {\frac {2G}{3\alpha Q}}{\frac {\partial {\dot {p}}}{\partial n}}\,d\Gamma =0}$
(2.64)

where ${\textstyle {\frac {\partial {\dot {p}}}{\partial n}}}$ is the derivative of ${\textstyle {\dot {p}}}$ in the direction of the normal to the external boundary. This derivative can be approximated as shown in Figure 6.

 Figure 6: Computation of the term ${\displaystyle {\frac {\partial {\dot {p}}}{\partial n}}}$ at the side ${\displaystyle ij}$ of a triangle adjacent to the external boundary ${\displaystyle \Gamma }$.

Using the previous argument, the stabilized FIC form for the mass balance equation (2.64) finally becomes

 ${\displaystyle \int _{\Omega }v\,\alpha {\dot {\epsilon }}\,d\Omega +\int _{\Omega }v\,{\frac {1}{Q}}{\dot {p}}\,d\Omega +\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}{\frac {1}{\mu }}k_{ij}{\frac {\partial p}{\partial x_{j}}}\,d\Omega }$ ${\displaystyle -\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}{\frac {1}{\mu }}k_{ij}\rho _{f}b_{j}\,d\Omega +\int _{\Gamma }v\,{\tilde {q}}_{n}\,d\Gamma -\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}\tau 2G{\frac {\partial {\dot {\varepsilon }}_{ij}}{\partial x_{j}}}\,d\Omega }$ ${\displaystyle -\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}\tau {\frac {\partial {\dot {\sigma }}'}{\partial x_{i}}}\,d\Omega +\int _{\Omega }{\frac {\partial v}{\partial x_{i}}}\tau \left(\alpha -{\frac {2G}{3\alpha Q}}\right){\frac {\partial {\dot {p}}}{\partial x_{i}}}\,d\Omega }$ ${\displaystyle -\int _{\Gamma }v\,{\frac {\tau }{3h_{n}}}{\frac {4G}{3\alpha Q}}{\dot {p}}\,d\Gamma }$
(2.65)

where ${\textstyle h_{n}}$ is an arbitrary distance in the normal direction to the boundary. In our work ${\textstyle h_{n}}$ has been taken as the characteristic length ${\textstyle l^{e}}$ of the element adjacent to the external boundary.

### 2.3.4 Discretized form of the momentum and FIC-stabilized mass balance equations

Combining the discretized form of Equation (2.65) along with the residual obtained from the original balance of momentum equation (2.24), the system of equations to be solved reads

• Balance of momentum
 ${\displaystyle {\boldsymbol {r}}_{u}={\boldsymbol {M}}{\ddot {\bar {\boldsymbol {u}}}}+\int _{\Omega }{\boldsymbol {B}}^{T}{\boldsymbol {\sigma }}'\,d\Omega {-}{\boldsymbol {Q}}{\bar {\boldsymbol {p}}}-{\boldsymbol {f}}_{u}={\boldsymbol {0}}}$
(2.66)
• Mass balance

 ${\displaystyle {\hat {\boldsymbol {r}}}_{p}=\left({\boldsymbol {Q}}^{T}-{\boldsymbol {R}}\right){\dot {\bar {\boldsymbol {u}}}}-{\boldsymbol {l}}+\left({\boldsymbol {C}}+{\boldsymbol {T}}-{\boldsymbol {T}}_{b}\right){\dot {\bar {\boldsymbol {p}}}}+{\boldsymbol {H}}{\bar {\boldsymbol {p}}}-{\boldsymbol {f}}_{p}={\boldsymbol {0}}}$
(2.67)

with the stabilizing terms

 ${\displaystyle {\boldsymbol {l}}=\int _{\Omega }\left({\boldsymbol {\nabla }}{\boldsymbol {N}}_{p}\right)^{T}{\frac {\tau }{3}}{\hat {\boldsymbol {S}}}^{T}{\dot {\boldsymbol {\sigma }}}'\,d\Omega \quad {\hbox{;}}\quad {\boldsymbol {R}}=\int _{\Omega }\left({\boldsymbol {\nabla }}{\boldsymbol {N}}_{p}\right)^{T}\tau 2G{\boldsymbol {S}}^{T}{\boldsymbol {I}}_{v}{\boldsymbol {B}}\,d\Omega }$
(2.68)
 ${\displaystyle {\boldsymbol {T}}=\int _{\Omega }\left({\boldsymbol {\nabla }}{\boldsymbol {N}}_{p}\right)^{T}\tau \left(\alpha -{\frac {2G}{3\alpha Q}}\right)\left({\boldsymbol {\nabla }}{\boldsymbol {N}}_{p}\right)\,d\Omega }$
(2.69)
 ${\displaystyle {\boldsymbol {T}}_{b}=\int _{\Gamma }{\boldsymbol {N}}_{p}^{T}{\frac {\tau }{h_{n}}}{\frac {4G}{3\alpha Q}}{\boldsymbol {N}}_{p}\,d\Gamma }$
(2.70)

where ${\textstyle {\hat {\boldsymbol {S}}}={\boldsymbol {m}}{\boldsymbol {\nabla }}^{T}}$ and the matrix ${\textstyle {\boldsymbol {I}}_{v}}$ appears due to the difference between the strain tensor and the strain vector expressed in Voigt notation. For a general 3D case

 ${\displaystyle {\boldsymbol {I}}_{v}={\frac {1}{2}}{\begin{bmatrix}2&&&&&0\\&2&&&&\\&&2&&&\\&&&1&&\\&&&&1&\\0&&&&&1\\\end{bmatrix}}}$
(2.71)

Note that since the test functions ${\textstyle v}$ are zero on the Dirichlet boundaries ${\textstyle \Gamma _{p}}$, matrix ${\textstyle {\boldsymbol {T}}_{b}}$ must be only computed on those external boundaries where no pressure Dirichlet conditions are imposed.

Using the stabilized residual vector (2.67), the system of equations for a linear elastic material, quasi-static and undrained-incompressible limit case becomes

 ${\displaystyle {\begin{bmatrix}{\boldsymbol {K}}&-{\boldsymbol {Q}}\\\displaystyle {\frac {\gamma }{\beta \Delta t}}\left[{\boldsymbol {Q}}^{T}-\left({\boldsymbol {R}}+{\boldsymbol {L}}\right)\right]&\displaystyle {\frac {1}{\theta \Delta t}}{\boldsymbol {T}}^{*}\end{bmatrix}}_{n+1}^{i}{\begin{bmatrix}\delta {\bar {\boldsymbol {u}}}\\\delta {\bar {\boldsymbol {p}}}\end{bmatrix}}_{n+1}^{i+1}=-{\begin{bmatrix}{\boldsymbol {r}}_{u}\\{\hat {\boldsymbol {r}}}_{p}\end{bmatrix}}_{n+1}^{i}}$
(2.72)

with

 ${\displaystyle {\boldsymbol {L}}=\int _{\Omega }\left({\boldsymbol {\nabla }}{\boldsymbol {N}}_{p}\right)^{T}{\frac {\tau }{3}}{\hat {\boldsymbol {S}}}^{T}{\boldsymbol {D}}{\boldsymbol {B}}\,d\Omega }$
(2.73)
 ${\displaystyle {\boldsymbol {T}}^{*}=\int _{\Omega }\left({\boldsymbol {\nabla }}{\boldsymbol {N}}_{p}\right)^{T}\tau \alpha \left({\boldsymbol {\nabla }}{\boldsymbol {N}}_{p}\right)\,d\Omega }$
(2.74)

Note that the diagonal of the iteration matrix is now different from zero. Matrices ${\textstyle {\boldsymbol {R}}}$ and ${\textstyle {\boldsymbol {L}}}$ emerging from the stabilized FIC formulation are essential for ensuring the consistency of the residual-based stabilization method used.

The consistency of residual does not directly affect the diffusion of the numerical solution, but it helps improving the condition number of the system matrix and facilitates the convergence of the Newton-Raphson scheme.

We highlight that the FIC-FEM formulation presented is applicable to equal order interpolation approximations of the displacements and the pressure of any degree. In this work we have implemented the FIC-FEM formulation using 3-noded triangles and 4-noded quadrilaterals for 2D problems, and 4-noded tetrahedra and 8-noded hexahedra for 3D problems, with equal order interpolation for the displacements and pressure.

## 2.4 Examples

This section includes two academic examples designed to assess the performance of the stabilized FIC-FEM formulation in situations near the undrained-incompressible limit. The first example is an elastic soil column subjected to a surface load. The second problem is an elastic soil foundation also subjected to a surface load.

The soil column problem is analysed in a 2D framework under plane strain conditions, whereas the soil foundation case is tackled as a 3D problem. The problems have been solved with simple elements with equal order interpolation for the displacements and the pore pressure, using the direct (non-stabilized) formulation and the FIC-stabilized formulation presented in the previous sections. Stable mixed elements with higher order interpolation for the displacements have been used to validate the FIC-stabilized element.

In both problems, the porous medium is under totally saturated conditions with isotropic permeability. The effect of gravity is not considered.

This example consists on a ${\textstyle 1\times 30}$ ${\textstyle m}$ column of saturated soil subjected to a surface loading that lies on a rigid rock bed. The geometry and boundary conditions of the problem are shown in Figure 7.

 Figure 7: Elastic soil column subjected to surface loading. Geometry and boundary conditions. Dimensions in ${\displaystyle m}$.

The material properties of the soil column are summarized in Table 1. It must be noticed that the indices 1 and 2 in the solid bulk modulus ${\textstyle K_{s}}$, the fluid bulk modulus ${\textstyle K_{f}}$ and the intrinsic permeability ${\textstyle k}$ denote the different cases that have been considered here. Index 1 corresponds to the nearly undrained-incompressible limit, and index 2 corresponds to a more relaxed condition.

Table 1. Elastic soil column problem. Material properties.
 Property Value Young modulus (${\textstyle E}$) ${\textstyle 2.5\cdot 10^{7}\;N/m^{2}}$ Poisson's ratio (${\textstyle \nu }$) 0.2 Solid density (${\textstyle \rho _{s}}$) ${\textstyle 2000\;Kg/m^{3}}$ Fluid density (${\textstyle \rho _{f}}$) ${\textstyle 1000\;Kg/m^{3}}$ Porosity (${\textstyle \phi }$) 0.3 Dynamic viscosity (${\textstyle \mu }$) ${\textstyle 0.001\;N/m^{2}\cdot s}$ Solid bulk modulus 1 (${\textstyle K_{s,1}}$) ${\textstyle 1.5\cdot 10^{17}\;N/m^{2}}$ Fluid bulk modulus 1 (${\textstyle K_{f,1}}$) ${\textstyle 3\cdot 10^{14}\;N/m^{2}}$ Intrinsic Permeability 1 (${\textstyle k_{1}}$) ${\textstyle 1\cdot 10^{-14}\;m^{2}}$ Solid bulk modulus 2 (${\textstyle K_{s,2}}$) ${\textstyle 1.5\cdot 10^{12}\;N/m^{2}}$ Fluid bulk modulus 2 (${\textstyle K_{f,2}}$) ${\textstyle 3\cdot 10^{9}\;N/m^{2}}$ Intrinsic Permeability 2 (${\textstyle k_{2}}$) ${\textstyle 1\cdot 10^{-11}\;m^{2}}$

This problem is solved in a 2D configuration under plane strain conditions. The geometry is discretized with a structured mesh of 20 quadrilateral elements. We have solved the problem using 4-noded quadrilateral elements with bilinear shape functions for both the pressure and the displacements (Q4P4), and 9-noded quadrilateral elements with biquadratic shape functions for the displacements and bilinear shape functions for the pressure (Q9P4).

The main purpose of this first case is to verify that the FIC-stabilized formulation captures properly the pressure distribution along the soil column under undrained-incompressible conditions.

 Figure 8: Surface step loading applied in the elastic soil column problem.

Figure 9 shows the normalized pore pressure along the normalized height of the soil column for nearly undrained conditions at a time ${\textstyle t=2}$ ${\textstyle s}$.

Understanding the results of Figure 9 requires taking into account the fact that gravity is not considered. Indeed, since the only force acting over the soil column is the surface load, the pressure applied by the surface load is completely transferred to the pore pressure throughout the column height and thus the ratio ${\textstyle p/q}$ should equal 1.

The results obtained with the non-stabilized formulation and the 4-noded quadrilaterals (hereafter called Q4P4-Direct) are not able to properly capture the pressure distribution along the column. The element locks and the pressure oscillates with arbitrary values. In contrast, by using the Q4P4 element stabilized with the FIC approach presented in this work (here called Q4P4-FIC), the correct pressure distribution is obtained. Similar good results are obtained with the Q9P4 element, as expected.

Comparing the graphs in Figures 9a and 9b one can see the effect of the compressibility of the materials on the sought solution. Certainly, as we decrease the value of Biot's modulus ${\textstyle Q}$, or equivalently, the bulk modulus of the solid and fluid phases ${\textstyle K_{s}}$ and ${\textstyle K_{f}}$, the amplitude of the pressure oscillations obtained with the Q4P4-Direct element decreases near the rock bed. On the other hand, the results obtained with the Q4P4-FIC and Q9P4 elements remain unaltered.

 (a) Quasi-incompressible limit (${\displaystyle Q=10^{15}N/m^{2}}$). (b) Relaxed incompressibility (${\displaystyle Q=10^{10}N/m^{2}}$). Figure 9: Normalized pore pressure along the soil column (${\displaystyle k=1\cdot 10^{-14}\;m^{2}}$, ${\displaystyle \Delta t=0.02}$ ${\displaystyle s}$, ${\displaystyle t=2}$ ${\displaystyle s}$).

In this second case the objective is to analyse the evolution of the pore pressure in time, and assess the effect of the soil permeability on the dissipation of the pore pressure.

 Figure 10: Surface cyclic loading applied in the elastic soil column problem.

To that purpose, Figure 11 shows the temporal evolution of the pore pressure at a node located 1.5 m below the surface.

Like in the previous case, the pressure obtained here should reflect the load applied on the surface. From the results in Figure 11, it is clear that the three tested elements show sinusoidal pressure evolutions in time. However, from the nearly undrained case of Figure 11a, it appears that while the Q4P4-FIC element shows a pressure that correctly oscillates from 500 to 1500 ${\textstyle N/m^{2}}$, the Q4P4-Direct element presents pressure values that vary from 1000 to 3000 ${\textstyle N/m^{2}}$, making no sense. One must notice that although the Q9P4 element shows higher pressures than the stabilized element, after refining the mesh up to 60 elements, the Q9P4 results are in good agreement with the solution obtained with 20 Q4P4-FIC elements (Figure 12). This evidences the good accuracy of the stabilized FIC-FEM Q4P4 element.

 (a) Quasi-undrained limit (${\displaystyle k=10^{-14}m^{2}}$). (b) Partially drained (${\displaystyle k=10^{-11}m^{2}}$). Figure 11: Evolution of the pore pressure with time at a node near the surface (${\displaystyle Q=1\cdot 10^{15}N/m^{2}}$, ${\displaystyle \Delta t=0.05}$ ${\displaystyle s}$).

Looking at the graph depicted in Figure 11b one can see that, by increasing the intrinsic permeability ${\textstyle k}$, the problem is no longer as bad conditioned as before, and so the gap in the amplitude of the three elements diminishes. Moreover, a greater permeability implies that the water can flow towards the surface due to the deformation of the column. Thereby, a dissipation in the pore pressure occurs as the computation advances in time.

Finally, it is also interesting to notice that the pressure becomes negative after the soil has drained a certain amount of water. This is easily understood if one thinks about the elastic behaviour considered for the porous medium. Since the applied load is oscillating, the deformation of the elastic soil experiments a similar behaviour. In consequence, when the soil recovers its initial position as the load reduces, a suction appears that makes the pore pressure negative at some points near the surface.

 Figure 12: Time evolution of the pore pressure. Comparison between 20 Q4P4-FIC elements and 60 stable Q9P4 elements (${\displaystyle Q=1\cdot 10^{15}N/m^{2}}$, ${\displaystyle k=10^{-14}m^{2}}$, ${\displaystyle \Delta t=0.05}$ ${\displaystyle s}$).

This second example consists on a pillar foundation on a saturated soil stratum lying over a rigid rock bed. The geometry and boundary conditions are sketched in Figure 13.

 Figure 13: Elastic soil foundation subjected to surface loading. Geometry and boundary conditions. Dimensions in ${\displaystyle m}$.

The material properties for the soil are the same as for the previous example (Table 2).

Table 2. Elastic soil foundation problem. Material properties.
 Property Value Young modulus (${\textstyle E}$) ${\textstyle 2.5\cdot 10^{7}\;N/m^{2}}$ Poisson's ratio (${\textstyle \nu }$) 0.2 Solid density (${\textstyle \rho _{s}}$) ${\textstyle 2000\;Kg/m^{3}}$ Fluid density (${\textstyle \rho _{f}}$) ${\textstyle 1000\;Kg/m^{3}}$ Porosity (${\textstyle \phi }$) 0.3 Dynamic viscosity (${\textstyle \mu }$) ${\textstyle 0.001\;N/m^{2}\cdot s}$ Solid bulk modulus 1 (${\textstyle K_{s,1}}$) ${\textstyle 1.5\cdot 10^{17}\;N/m^{2}}$ Fluid bulk modulus 1 (${\textstyle K_{f,1}}$) ${\textstyle 3\cdot 10^{14}\;N/m^{2}}$ Intrinsic Permeability 1 (${\textstyle k_{1}}$) ${\textstyle 1\cdot 10^{-14}\;m^{2}}$ Solid bulk modulus 2 (${\textstyle K_{s,2}}$) ${\textstyle 1.5\cdot 10^{12}\;N/m^{2}}$ Fluid bulk modulus 2 (${\textstyle K_{f,2}}$) ${\textstyle 3\cdot 10^{9}\;N/m^{2}}$ Intrinsic Permeability 2 (${\textstyle k_{2}}$) ${\textstyle 1\cdot 10^{-11}\;m^{2}}$

In this case the soil is subjected to a surface step loading of ${\textstyle 1\cdot 10^{4}N/m^{2}}$ applied in a linear increment during ${\textstyle 0.1}$ ${\textstyle s}$ (Figure 14).

 Figure 14: Surface step loading applied in the elastic soil foundation problem.

As stated before, this problem is solved by means of a 3D analysis. We use 4-noded tetrahedral elements with linear interpolation for the displacements and the pressure in the non-stabilized form (T4P4-Direct) and the FIC-stabilized formulation (T4P4-FIC). Ten-noded stable tetrahedral elements with quadratic shape functions for the displacement field and linear shape functions for the pressure field (T10P4) have also been used in the analysis for comparison purposes.

This example has two objectives. First, it should help analysing the effect of the spatial discretization on the solution obtained with the FIC-stabilized formulation. Second, it may allow us to assess whether the three element types considered here converge to the expected solution when the mesh is refined.

In order to do so, two different unstructured spatial discretizations have been used: a coarse uniform mesh (Figure 15a) and a refined non-uniform mesh (Figure 15).

 (a) Coarse uniform mesh: 9985 elements. (b) Refined non-uniform mesh: 13532 elements. Figure 15: Spatial discretizations used for the elastic soil foundation problem.

In Figure 16 we represent the evolution of the maximum pore pressure with time under nearly undrained-incompressible conditions for both the coarse and refined meshes using ${\textstyle \Delta t=0.02}$ ${\textstyle s}$. The impermeability of the case makes that the pore pressure remains constant after the first ${\textstyle 0.1}$ ${\textstyle s}$ of application of the load for the two stable elements, the T10P4 and the T4P4-FIC. We see, though, a different behaviour in the unstable T4P4-Direct element, showing a certain dissipation for the finer mesh.

Comparing the T10P4 and the T4P4-FIC elements we notice that, while the first one surpasses the threshold of ${\textstyle 10000}$ ${\textstyle N/m^{2}}$ for both meshes, the latter is always below the value of the applied load. In any case, both elements approach the expected solution as the mesh is refined, which indicates a consistent response.

 (a) Coarse mesh. (b) Refined mesh. Figure 16: Time evolution of the maximum pore pressure under undrained-incompressible conditions (${\displaystyle Q=1\cdot 10^{15}\;N/m^{2}}$, ${\displaystyle k=1\cdot 10^{-14}\;m^{2}}$, ${\displaystyle \Delta t=0.02}$ ${\displaystyle s}$).

If we look at the contour lines of the pore pressure distribution obtained at time ${\textstyle t=1}$ ${\textstyle s}$ in Figure 17, we can understand the abnormal behaviour of the T4P4-Direct element. Indeed, the direct mixed formulation with equal order interpolation elements locks under undrained-incompressible conditions which leads to the incoherent pressure distribution obtained.

 (a) T4P4-Direct elem. in coarse mesh (${\displaystyle p_{max}=16282N/m^{2}}$). (b) T4P4-Direct elem. in refined mesh (${\displaystyle p_{max}=15851N/m^{2}}$). (c) T10P4 elem. in coarse mesh (${\displaystyle p_{max}=15197N/m^{2}}$). (d) T10P4 elem. in refined mesh (${\displaystyle p_{max}=12283N/m^{2}}$). (e) T4P4-FIC elem. in coarse mesh (${\displaystyle p_{max}=4390.8N/m^{2}}$). (f) T4P4-FIC elem. in refined mesh (${\displaystyle p_{max}=9609.7N/m^{2}}$). Figure 17: Pore pressure distribution under undrained-incompressible conditions (${\displaystyle Q=1\cdot 10^{15}\;N/m^{2}}$, ${\displaystyle k=1\cdot 10^{-14}\;m^{2}}$, ${\displaystyle \Delta t=0.02}$ ${\displaystyle s}$, ${\displaystyle t=1}$ ${\displaystyle s}$).

Next we evaluate the capabilities of the FIC stabilized formulation to reproduce a solution where stabilization is not needed. For this purpose the previous problem has been solved for partially drained and compressible conditions up to ${\textstyle t=5}$ s (Figure 18). Results are again shown at ${\textstyle t=1}$ s for ${\textstyle \Delta t=0.02}$ s (Figure 19). Under these relaxed conditions, the equal order interpolation mixed element with the direct formulation is able to capture acceptable pressure distributions as expected (Figures 19a and 19b).

The higher permeability of this second case also permits capturing the so-called Mandel-Cryer effect, i.e. the increase of pore pressure due to the application of the load is immediate, but the dissipation due to the outflow of the pore fluid is delayed by the permeability of the porous medium [110]. Figures 18a and 18b show that, after the first 0.1 s of loading, the fluid starts draining due to the consolidation of the soil.

In order to properly understand the low peak pore pressure obtained with the T4P4-FIC element in Figure 18a, it is useful to recall the definition of the stabilization parameter in Equation (4.56) and observe carefully the mesh in Figure 15a. It can be noticed that the area of application of the load is covered by only two triangular faces. This implies that the characteristic element length around that particular zone is relatively large, and so it is the stabilization parameter. As a result, one obtains an over-diffusive solution that leads to an underestimation of the maximum pore pressure.

However, it is interesting to see that the larger peak pore pressure obtained with the T4P4-Direct and the T10P4 elements, caused by an initial locking of the pressure field, converges to the solution obtained with the T4P4-FIC element after the soil consolidates.

On the other hand, the results for the maximum pressure obtained with the refined non-uniform mesh (Figures 18b, 19b, 19d and 19f) converge all to an almost identical value. This evidences that the FIC stabilization does not alter negatively the solution obtained with the original non-stabilized mixed formulation.

 (a) Coarse mesh. (b) Refined mesh. Figure 18: Time evolution of the maximum pore pressure for partially compressible and drained conditions (${\displaystyle Q=1\cdot 10^{10}\;N/m^{2}}$, ${\displaystyle k=1\cdot 10^{-11}\;m^{2}}$, ${\displaystyle \Delta t=0.02}$ ${\displaystyle s}$).
 (a) T4P4-Direct elem. in coarse mesh (${\displaystyle p_{max}=7268.1N/m^{2}}$). (b) T4P4-Direct elem. in refined mesh (${\displaystyle p_{max}=8104.1N/m^{2}}$). (c) T10P4 elem. in coarse mesh (${\displaystyle p_{max}=10427N/m^{2}}$). (d) T10P4 elem. in refined mesh (${\displaystyle p_{max}=8173.8N/m^{2}}$). (e) T4P4-FIC elem. in coarse mesh (${\displaystyle p_{max}=4089N/m^{2}}$). (f) T4P4-FIC elem. in refined mesh (${\displaystyle p_{max}=8086.8N/m^{2}}$). Figure 19: Pore pressure distribution for partially compressible and drained conditions (${\displaystyle Q=1\cdot 10^{10}\;N/m^{2}}$, ${\displaystyle k=1\cdot 10^{-11}\;m^{2}}$, ${\displaystyle \Delta t=0.02}$ ${\displaystyle s}$, ${\displaystyle t=1}$ ${\displaystyle s}$).

# 3 Continuum Damage Mechanics

## 3.1 Introduction

Many engineering materials subjected to unfavourable mechanical and environmental conditions undergo micro-structural changes that decrease their strength. Since these changes impair the mechanical properties of these materials, the term damage is generally used.

This effect is particularly relevant and difficult to predict in brittle or quasi-brittle materials such as concrete, rocks, mortar or other geo-materials. For instance, concrete is a highly heterogeneous, anisotropic, brittle material with a very complex non-linear mechanical behaviour due to the occurrence of the localization of deformation. This localization of deformation can appear as cracks, if cohesive properties are dominant, or as shear zones, if frictional properties prevail. As a result of strain localization, material softening takes place and a significant reduction of the material stiffness occurs during cyclic loading. A good understanding of the mechanism of the formation of such localization is of crucial importance because it acts as a precursor of fracture and failure.

Continuum damage mechanics is a branch of continuum mechanics that describes the progressive loss of material integrity due to the propagation and coalescence of micro-cracks, micro-voids, and similar defects (Figure 20). These changes in the micro-structure lead to an irreversible material degradation, characterized by a loss of stiffness that can be observed on the macro-scale.

 Figure 20: Generic damaged section.

The main objective of standard continuum damage mechanics is to propose a continuum-mechanics based framework that allows to characterize, represent and model at the macroscopic scale the effects of distributed defects and their growth on the material behaviour.

The term “continuum damage mechanics” was first used by Hult [111], but the concept of damage was introduced by Kachanov in 1958 in the context of creep rupture [53]. In that work Kachanov introduced the concept of effective stress, and by using continuum damage he solved problems related to creep in metals. Rabotnov [54] gave the problem physical meaning by suggesting that the reduction of the sectional area was measured by means of the damage parameter.

The thermodynamic formalism involved in the irreversible process of damage was developed by Lemaitre and Chaboche [55]. Other important contributions on damage mechanics include: Mazars and Pijaudier-Cabot [112,56], Simo and Ju [57], Oller et al. [58], Oliver et al. [59,60], and Cervera et al. [61,62,113].

If the damage parameter depends only on the strain state at the point under consideration, and non enriched kinematics are adopted to regularize the problem, numerical simulations exhibit a pathological mesh dependence and the energy consumed by the fracture process tends to zero as the mesh is refined. This is the typical behaviour of the so-called local damage models, which are not able to properly describe both the thickness of localization and the distance between damaged zones. They suffer from mesh sensitivity (for size and alignment) and produce unreliable results. Strains concentrate in one element wide zones and the computed force-displacement curves are mesh-dependent. The reason behind these misbehaviours is that the differential equations of motion change their type, from elliptic to hyperbolic in static problems, and the boundary value problem becomes ill-posed [114].

Classical constitutive models require an extension in the form of a characteristic length to properly model the thickness of localized zones. Such extension can be done by means of second gradient models [63], micro-polar [64], strain gradient [65], viscous [66] and non-local terms [67]. In this monograph we have worked with the latter approach using a weighted spatial averaging of the internal variables. In this manner the kinematic and equilibrium equations remain standard, and the notions of stress and strain keep their usual meaning.

The first non-local models of this type, proposed in the 1960s, aimed at improving the description of elastic wave dispersions in crystals. Non-local elasticity was further developed by Eringen [68] who later extended it to non-local elasto-plasticity [69]. Subsequently, it was found that certain non-local formulations could act as efficient localization limiters with a regularizing effect on problems with strain localization [70].

The chapter is organized as follows. First, the basic concepts on continuum damage mechanics are introduced. Details are given on the essential components of the isotropic damage theory, including the equivalent strain forms and damage evolution laws implemented for this work. Following that, the regularization techniques developed to overcome the difficulties associated to strain localization are illustrated. Special attention is given to the integral-type non-local damage model, pointing out the most relevant aspects of its numerical implementation. Two examples are shown at the end, in order to assess the robustness of the implemented damage models against changes in the spatial discretization.

## 3.2 Isotropic damage theory

Damage models work with certain internal variables that characterize the density and orientation of micro-defects. To introduce its concepts and understand the fundamental theory of damage mechanics, let us explain the stress evolution in a simple uniaxial tensile case (Figure 21).

 Figure 21: Idealized material for the description of the uniaxial damage theory.

For a better understanding, the material is idealized as a bundle of fibres parallel to the loading direction (Figure 21). Initially, all the fibres respond elastically, and the stress is carried by the total cross section of all fibres ${\textstyle A}$ (Figure 22.1). As the applied strain is increased some fibres start breaking (Figure 22.2). Each fibre is assumed to be perfectly brittle, meaning that the stress in the fibre drops down to zero immediately after a critical strain level is reached. However, since the critical strain can differ from one fibre to another, the effective area ${\textstyle A_{e}}$ (the area of unbroken fibres that can still carry stress) decreases gradually from ${\textstyle A_{e}=A}$ to ${\textstyle A_{e}=0}$. Of course, when the applied force diminishes (Figure 23.2), the affected fibres remain broken and so the damaged area of the specimen is irrecoverable.

 Figure 22: Scheme of a uniaxial damage model through a monotonic loading process.
 Figure 23: Scheme of a uniaxial damage model through a non-monotonic loading process.

It is important to make a distinction between the nominal stress ${\textstyle \sigma }$, defined as the force per unit initial area of the cross section, and the effective stress ${\textstyle \sigma _{e}}$, defined as the force per unit effective area1. The nominal stress enters the Cauchy equations of equilibrium on the macroscopic level, while the effective stress is the “actual” stress acting in the material micro-structure. Since the external applied force is the same regardless of the definition of the stress, i.e. ${\textstyle F=\sigma A=\sigma _{e}A_{e}}$, we can write:

 ${\displaystyle \sigma ={\frac {A_{e}}{A}}\sigma _{e}}$
(3.1)

The ratio of the effective area to the total area, ${\textstyle A_{e}/A}$, is a scalar characterizing the integrity of the material. Within the classical approach, a very simple measure of the damage amplitude in a given plane is obtained by measuring the area of the intersection of all defects with that plane. Thereby, we can define the damage variable as:

 ${\displaystyle \mathbb {d} =1-{\frac {A_{e}}{A}}={\frac {A-A_{e}}{A}}={\frac {A_{d}}{A}}}$
(3.2)

where ${\textstyle A_{d}}$ is the damaged part of the area. An undamaged material is characterized then by ${\textstyle A_{e}=A}$, i.e. by ${\textstyle \mathbb {d} =0}$. Due to propagation of micro-defects and their coalescence, the damage variable grows and at the late stages of degradation process it approaches asymptotically the limit value ${\textstyle \mathbb {d} =1}$, corresponding to a complete damaged material with effective area reduced to zero.

In the simplest version of the presented scheme, each fibre is supposed to remain linear elastic up to the strain level at which it breaks. Consequently, the 1D effective stress ${\textstyle \sigma _{e}}$ is related to the elastic strain of the material ${\textstyle \varepsilon }$ by the uniaxial Hooke's law:

 ${\displaystyle \sigma _{e}=E\varepsilon }$
(3.3)

where ${\displaystyle E}$ is the elastic modulus of the undamaged material. Combining (2.1), (3.2) and (3.3), it is straightforward to see that the constitutive law for the nominal stress ${\textstyle \sigma }$ takes the form:

 ${\displaystyle \sigma =(1-\mathbb {d} )E\varepsilon }$
(3.4)

For the uniaxial model formulation, Equation (3.4) must be complemented with the damage evolution law, which can be characterized by the dependence between the damage variable and the applied strain:

 ${\displaystyle {\begin{array}{cc}\mathbb {d} =g(\varepsilon )&\quad 0\leq \mathbb {d} \leq 1\end{array}}}$
(3.5)

Function ${\textstyle g}$ affects the shape of the stress-strain diagram and can be directly identified from a uniaxial test.

The evolution of the effective stress, damage variable, and nominal stress in a material that remains elastic up to the peak stress is shown in Figure 22. This description is valid only for monotonic loading by an increasing applied strain ${\textstyle \varepsilon }$.

When the material is first stretched up to a certain strain level ${\textstyle \varepsilon _{1}}$ that induces damage ${\textstyle \mathbb {d} _{1}=g(\varepsilon _{1})}$ and then the strain decreases (Figure 23), the damaged area remains constant and the material responds as an elastic material with a reduced Young's modulus ${\textstyle E_{1}=(1-\mathbb {d} _{1})E}$. This means that, during unloading and partial reloading, the damage variable in (5.4) must be evaluated from the largest previously reached strain and not from the current strain ${\textstyle \varepsilon }$. This means that it is crucial to introduce an internal state variable ${\textstyle r}$ characterizing the maximum strain level reached in the previous history of the material up to a given time ${\textstyle t}$:

 ${\displaystyle r(t)=\max \left\{r_{y},\max _{\tau \leq t}\varepsilon (\tau )\right\}}$
(3.6)

The above expression implies that ${\textstyle r(t)\geq r_{y}}$, where ${\textstyle r_{y}}$ is the so-called damage threshold, a material parameter that represents the value of strain at which damage starts. In this formula, ${\textstyle t}$ is not necessarily the physical time (it can be any monotonically increasing parameter controlling the loading process). The damage evolution law (3.5) is then rewritten in the form:

 ${\displaystyle \mathbb {d} =g(r){\hbox{ with }}\left\{{\begin{array}{cc}g(r)=0&{\hbox{ if }}r=r_{y}\\0r_{y}\end{array}}\right.}$
(3.7)

which remains valid not only during monotonic loading but also during unloading and reloading. The evolution of the effective stress, damage variable, and nominal stress in a non-monotonic test is shown in Figure 23. Note that, upon a complete removal of the applied force, the strain returns to zero (due to elasticity of the unbroken fibres), i.e. the pure damage model does not take into account any permanent strains. Nevertheless, the material state is different from the initial virgin state because, according to (3.6) and (3.7), when the state variable ${\textstyle r}$ becomes greater than ${\textstyle r_{y}}$, the damage variable will not be zero again and thus the stiffness and strength mobilized in a new tensile loading process will be smaller than their initial values. Therefore, we can say that the loading history is always reflected by the value of the internal state variable ${\textstyle r}$.

Instead of defining the variable ${\textstyle r}$ like in (3.6), we can introduce a loading function ${\textstyle f(\varepsilon ,r)=\varepsilon {-}r}$ and postulate the loading-unloading conditions in the Kuhn-Tucker form:

 ${\displaystyle {\begin{array}{ccc}f\leq {0};&{\dot {r}}\geq {0};&{\dot {r}}f=0\end{array}}}$
(3.8)

The first condition indicates that ${\textstyle r}$ can never be smaller than ${\textstyle \varepsilon }$, while the second condition means that ${\textstyle r}$ cannot decrease. Finally, according to the third condition, ${\textstyle r}$ can grow only if the current values of ${\textstyle \varepsilon }$ and ${\textstyle r}$ are equal.

At this point, we can already formulate the extension of the uniaxial damage theory to general multiaxial stress states.

In this work we have chosen the simple isotropic damage model with a unique scalar variable. Isotropic damage models are based on the simplifying assumption that the stiffness degradation is isotropic, i.e. the stiffness moduli corresponding to different directions decrease proportionally, independently of the direction of loading. Since an isotropic elastic material is characterized by two independent elastic constrains, a general isotropic damage model should deal with two damage variables. The model with a single variable makes use of the additional assumption that the relative reduction of all the stiffness coefficients is the same, in other words, that the Poisson's ratio is not affected by damage. In this regard, the stress-strain law in Voigt notation is postulated as:

 ${\displaystyle {\boldsymbol {\sigma }}=(1-\mathbb {d} ){\boldsymbol {D}}_{e}{\boldsymbol {\varepsilon }}}$
(3.9)

where ${\textstyle {\boldsymbol {D}}_{e}}$ is the elastic constitutive matrix of the intact material. One can clearly notice that (3.9) is a generalization of (3.4). Also, Equation (3.9) can alternatively be written as:

 ${\displaystyle {\boldsymbol {\sigma }}=(1-\mathbb {d} ){\boldsymbol {\sigma }}_{e}}$
(3.10)

which is the multidimensional generalization of (3.1), and where ${\textstyle {\boldsymbol {\sigma }}_{e}}$ is the effective stress vector defined as:

 ${\displaystyle {\boldsymbol {\sigma _{e}}}={\boldsymbol {D}}_{e}{\boldsymbol {\varepsilon }}}$
(3.11)

Similarly as in the uniaxial case, we introduce a loading function ${\textstyle f}$ specifying the elastic domain and the states at which damage grows. The loading function now depends on the strain vector ${\textstyle {\boldsymbol {\varepsilon }}}$, and on a variable ${\textstyle r}$ that controls the evolution of the elastic domain. Physically. ${\textstyle r}$ represents a scalar measure of the largest strain level ever reached in the history of the material. States for which ${\textstyle f({\boldsymbol {\varepsilon }},r)<0}$ are supposed to be below the current damage threshold. Damage can grow only if the current state reaches the boundary of the elastic domain, i.e. when ${\textstyle f({\boldsymbol {\varepsilon }},r)=0}$. Essentially, we can postulate the damage criterion for a multiaxial isotropic damage model with:

 ${\displaystyle f({\boldsymbol {\varepsilon }},r)=\varepsilon _{eq}({\boldsymbol {\varepsilon }})-r}$
(3.12)

${\displaystyle \varepsilon _{eq}}$ in Eq (3.12) is the equivalent strain, a scalar magnitude of the strain level, and ${\textstyle r}$ is the largest value of equivalent strain calculated in the previous deformation history of the material up to its current state. In this regard, (3.6) can now be generalized to:

 ${\displaystyle r(t)=\max \left\{r_{y},\max _{\tau \leq t}\varepsilon _{eq}(\tau )\right\}}$
(3.13)

An important advantage of an isotropic damage model is that the stress evaluation algorithm is usually explicit and there is no need to use an iterative solution for non-linear equations.

Thereby, for a particular strain increment, the corresponding stress is obtained by computing the current value of equivalent strain, updating the maximum previously reached equivalent strain and the damage variable, and reducing the effective stress according to (3.10). In essence, one must follow the scheme of Table 3.

 Step ${\displaystyle n+1}$ 1. Evaluate effective stress: ${\textstyle {\boldsymbol {\sigma }}_{e,n+1}={\boldsymbol {D}}_{e}{\boldsymbol {\varepsilon }}_{n+1}}$ 2. Compute the new equivalent strain: ${\textstyle \varepsilon _{eq,n+1}}$ 3. Update ${\textstyle r}$ with (3.13): If ${\textstyle \varepsilon _{eq,n+1}>r_{n}\Rightarrow r_{n+1}=\varepsilon _{eq,n+1}}$ 4. Update damage variable: ${\textstyle \mathbb {d} _{n+1}=g(r_{n+1})}$ 5. Compute nominal stress: ${\textstyle {\boldsymbol {\sigma }}_{n+1}=(1-\mathbb {d} _{n+1}){\boldsymbol {\sigma }}_{e,n+1}}$

(1) Note that this concept is different from Terzaghi's effective stress ${\displaystyle \sigma '}$ defined in Eq (2.1).

### 3.2.1 Elastic-damage tangent constitutive tensor

The material non-linearity introduced by the damage model requires solving the system of equations by means of incremental-iterative techniques such as the Newton-Raphson's method. If one aims at achieving reasonable convergence rates, it is crucial to properly compute the tangent constitutive tensor.

To that purpose, it is essential to understand that the global stiffness of the structure changes through any damage process. In 1D problems one can easily distinguish the secant Young's modulus ${\textstyle E_{sec}}$ from the tangent one ${\textstyle E_{tan}}$ and compute them as:

 ${\displaystyle E_{sec}={\frac {\sigma }{\varepsilon }}\quad ;\quad E_{tan}={\frac {d\sigma }{d\varepsilon }}}$
(3.14)

Figure 24 shows the stress-strain relation on a unidimensional damage process. The branch on the left side of the curve corresponds to the elastic loading case. The stress increases proportionally to the strain growth, and goes back following the same path when the strain diminishes. The secant and tangent elastic modulus coincide through all this stage. After a certain strain value is reached, i.e. ${\textstyle \varepsilon {>}r_{y}}$, the relation between the stress and the strain becomes non-linear as the material starts degrading. In this situation, the secant and tangent Young's moduli diverge and the internal state variable starts increasing (${\textstyle dr>0}$). Finally, since no plastic deformation is considered, the strain is fully recovered upon unloading following a branch with a lower stiffness than the original.

 Figure 24: Unidimensional stress-strain diagram throughout a damage process.

Understanding the previous concepts is not so straightforward in a general 3D case and so it is convenient to derive a general expression for the elastic-damage tangent constitutive tensor ${\textstyle {\boldsymbol {D}}_{tan}\equiv {\boldsymbol {D}}}$.

Let us remind first that the tangent constitutive tensor defines the tangent stiffness matrix ${\textstyle {\boldsymbol {K}}}$ (2.39) and must satisfy the following relation:

 ${\displaystyle d{\boldsymbol {\sigma }}={\boldsymbol {D}}d{\boldsymbol {\varepsilon }}}$
(3.15)

Considering the stress vector defined in (3.9) as a function of the strain vector and the damage variable, i.e. ${\textstyle {\boldsymbol {\sigma }}={\boldsymbol {\sigma }}({\boldsymbol {\varepsilon }},\mathbb {d} )}$, we compute its total derivative as:

 ${\displaystyle d{\boldsymbol {\sigma }}={\frac {\partial {\boldsymbol {\sigma }}}{\partial {\boldsymbol {\varepsilon }}}}d{\boldsymbol {\varepsilon }}+{\frac {\partial {\boldsymbol {\sigma }}}{\partial \mathbb {d} }}d\mathbb {d} =\left(1-\mathbb {d} \right){\boldsymbol {D}}_{e}d{\boldsymbol {\varepsilon }}-{\boldsymbol {D}}_{e}{\boldsymbol {\varepsilon }}d\mathbb {d} =\left(1-\mathbb {d} \right){\boldsymbol {D}}_{e}d{\boldsymbol {\varepsilon }}-{\boldsymbol {\sigma }}_{e}d\mathbb {d} }$
(3.16)

Like in the 1D case (Figure 24), we distinguish two possible situations, bearing in mind that the Kuhn-Tucker conditions hold (Eq (3.8)):

• Elastic loading or unloading process (${\textstyle f({\boldsymbol {\varepsilon }},r)<0\;\&\;dr=0}$)

Whenever the structure is under an elastic loading or unloading process, there is no change in the internal historical variable ${\textstyle r}$ (3.13). As a result, the damage variable ${\textstyle \mathbb {d} =g(r)}$ also remains constant (${\textstyle d\mathbb {d} =0}$) and so the expression in (3.16) results:

 ${\displaystyle d{\boldsymbol {\sigma }}=\left(1-\mathbb {d} \right){\boldsymbol {D}}_{e}d{\boldsymbol {\varepsilon }}}$
(3.17)

Consequently, in this case the tangent constitutive tensor coincides with the secant constitutive tensor:

 ${\displaystyle {\boldsymbol {D}}={\boldsymbol {D}}_{sec}=\left(1-\mathbb {d} \right){\boldsymbol {D}}_{e}}$
(3.18)
• Loading process with growing damage (${\textstyle f({\boldsymbol {\varepsilon }},r)=0\;\&\;dr>0}$)

If the structure is under a loading process with the equivalent strain exceeding the damage threshold ${\textstyle r_{y}}$, the internal state variable increases with time, and so the damage variable (${\textstyle d\mathbb {d} >0}$).

From the definitions in (3.7), (3.12) and (3.13), we compute the total derivative of the damage variable as:

 ${\displaystyle d\mathbb {d} ={\frac {dg}{dr}}dr={\frac {dg}{dr}}d\varepsilon _{eq}={\frac {dg}{dr}}{\frac {d\varepsilon _{eq}}{d{\boldsymbol {\varepsilon }}}}\cdot d{\boldsymbol {\varepsilon }}}$
(3.19)

Where the derivative of the equivalent strain with respect to the strain vector ${\textstyle d\varepsilon _{eq}/d{\boldsymbol {\varepsilon }}}$ is estimated by means of the perturbation method [115]. This method is actually convenient as it allows us to compute the previous derivative regardless of the form of the equivalent strain:

 ${\displaystyle {\frac {d\varepsilon _{eq}}{d\varepsilon _{i}}}={\frac {\varepsilon _{eq}\left({\boldsymbol {\varepsilon }}+{\boldsymbol {\Delta \varepsilon }}_{i}\right)-\varepsilon _{eq}\left({\boldsymbol {\varepsilon }}-{\boldsymbol {\Delta \varepsilon }}_{i}\right)}{2\Delta \varepsilon _{i}}}}$
(3.20)

with ${\textstyle {\boldsymbol {\Delta \varepsilon }}_{i}}$ being the perturbation vector of the ${\textstyle i^{th}}$ strain component.

 ${\displaystyle {\boldsymbol {\Delta \varepsilon }}_{i}={\begin{bmatrix}0&\dots &\Delta \varepsilon _{i}&\dots &0\end{bmatrix}}^{T}}$
(3.21)

Now, introducing Equation (3.19) into (3.16) and rearranging terms gives

 ${\displaystyle d{\boldsymbol {\sigma }}=\left[\left(1-\mathbb {d} \right){\boldsymbol {D}}_{e}-{\frac {dg}{dr}}{\boldsymbol {\sigma }}_{e}\times {\frac {d\varepsilon _{eq}}{d{\boldsymbol {\varepsilon }}}}\right]d{\boldsymbol {\varepsilon }}}$
(3.22)

And so the expression defining the tangent constitutive tensor reads

 ${\displaystyle {\boldsymbol {D}}=\left(1-\mathbb {d} \right){\boldsymbol {D}}_{e}-{\frac {dg}{dr}}{\boldsymbol {\sigma }}_{e}\times {\frac {d\varepsilon _{eq}}{d{\boldsymbol {\varepsilon }}}}={\boldsymbol {D}}_{sec}-{\frac {dg}{dr}}{\boldsymbol {\sigma }}_{e}\times {\frac {d\varepsilon _{eq}}{d{\boldsymbol {\varepsilon }}}}}$
(3.23)

The typical scheme for the numerical implementation of the tangent constitutive tensor is summarized in Table (4).

 Step ${\displaystyle n+1}$ 1. Compute secant constitutive tensor: ${\displaystyle {\boldsymbol {D}}_{sec,n+1}=(1-\mathbb {d} _{n+1}){\boldsymbol {D}}_{e}}$ If ${\textstyle \varepsilon _{eq,n+1}\geq r_{n}}$ 2. Compute damage function derivative: ${\textstyle (dg/dr)_{n+1}}$ 3. Calculate equivalent strain derivative: ${\textstyle (d\varepsilon _{eq}/d{\boldsymbol {\varepsilon }})_{n+1}}$ 4. Compute tangent constitutive tensor: ${\displaystyle {\boldsymbol {D}}_{n+1}={\boldsymbol {D}}_{sec,n+1}-({\frac {dg}{dr}})_{n+1}{\boldsymbol {\sigma }}_{e,n+1}\times ({\frac {d\varepsilon _{eq}}{d{\boldsymbol {\varepsilon }}}})_{n+1}}$

### 3.2.2 Equivalent strain

To some extent, the equivalent strain presented in (3.12) plays a role similar to the yield function in plasticity, because it directly affects the shape of the elastic domain (Figure 25).

 Figure 25: 3D elastic domain for a generic equivalent strain.

There are numerous forms for the equivalent strain in the literature, but the simplest choice is to define it as the Euclidean norm of the strain vector:

 ${\displaystyle \varepsilon _{eq}=\|{\boldsymbol {\varepsilon }}\|={\sqrt {{\boldsymbol {\varepsilon }}\cdot {\boldsymbol {\varepsilon }}}}={\sqrt {{\boldsymbol {\varepsilon }}^{T}{\boldsymbol {\varepsilon }}}}}$
(3.24)

or as the energy norm:

 ${\displaystyle \varepsilon _{eq}={\sqrt {\frac {{\boldsymbol {\varepsilon }}^{T}{\boldsymbol {D}}_{e}{\boldsymbol {\varepsilon }}}{E}}}}$
(3.25)

The two norms of ${\textstyle {\boldsymbol {\varepsilon }}}$ introduced above, lead to symmetric elastic domain in tension and compression (Figure 26a). Nevertheless, several materials (rocks, concrete, ceramics, etc.) often show a non symmetric damage surface, i.e., the yield value in compression can be several times the value in tension. In order to overcome this limitation, it is necessary to adjust the definition of the equivalent strain.

 (a) Symmetric energy norm. (b) Adapted Simo & Ju norm. (c) Mazars norm. (d) Modified von Mises norm. Figure 26: Damage surfaces in the 2D principal stress space.

In this regard, one may work wit a definition of the equivalent strain based on the model proposed by Simo and Ju [57], using the energy norm of the strain and modifying it to take into account the different degradation rate in tension and compression of concrete-like materials (Figure 26b). Such definition takes the form:

 ${\displaystyle \varepsilon _{eq}=\left(\theta +{\frac {1-\theta }{\kappa }}\right){\sqrt {{\boldsymbol {\varepsilon }}^{T}{\boldsymbol {D}}_{e}{\boldsymbol {\varepsilon }}}}}$
(3.26)

where the variable ${\textstyle \theta }$ is a weighting factor computed from the eigenvalues of the effective stress tensor:

 ${\displaystyle \theta ={\frac {\sum _{i=1}^{3}\langle \sigma _{i}\rangle }{\sum _{i=1}^{3}|\sigma _{i}|}}}$
(3.27)

with ${\textstyle \langle .\rangle }$ denoting the “positive part” operator or McAuley brackets. For scalars ${\textstyle \langle x\rangle =\max(0,x)}$, i.e.,

 ${\displaystyle \langle x\rangle =\left\{{\begin{array}{ll}x&{\hbox{ if }}x>0\\0&{\hbox{ if }}x\leq {\hbox{ 0}}\end{array}}\right.}$
(3.28)

The parameter ${\textstyle \kappa }$ in (3.26) is defined as the ratio between the compression elastic limit and the tension elastic limit, i.e.

 ${\displaystyle \kappa ={\frac {\sigma _{y_{c}}}{\sigma _{y_{t}}}}}$
(3.29)

Eq (3.26) is not the only possible form of equivalent strain valid for quasi-brittle materials. Micro-cracks in concrete grow mainly if the material is stretched, and so it is natural to take into account only strains that are positive (tensile strains) and neglect those that are negative (compressive strains). This leads to the so-called Mazars definition of equivalent strains [116]:

 ${\displaystyle \varepsilon _{eq}={\sqrt {\sum _{i=1}^{3}\langle \varepsilon _{i}\rangle ^{2}}}}$
(3.30)

where ${\textstyle \varepsilon _{i}}$ are the principal values of the strain tensor.

Alternatively, the modified von Mises definition [117], reads:

 ${\displaystyle \varepsilon _{eq}={\frac {\kappa {-1}}{2\kappa (1-2\nu )}}I_{1}+{\frac {1}{2\kappa }}{\sqrt {\left({\frac {\kappa {-1}}{1-2\nu }}I_{1}\right)^{2}+{\frac {12\kappa }{(1+\nu )^{2}}}J_{2}}}}$
(3.31)

where ${\textstyle I_{1}}$ is the first invariant of the strain tensor and ${\textstyle J_{2}}$ is the second invariant of the deviatoric strain tensor. Given a generic symmetric strain tensor:

 ${\displaystyle [\varepsilon ]={\begin{bmatrix}\varepsilon _{xx}&\varepsilon _{xy}&\varepsilon _{xz}\\\varepsilon _{xy}&\varepsilon _{yy}&\varepsilon _{yz}\\\varepsilon _{xz}&\varepsilon _{yz}&\varepsilon _{zz}\end{bmatrix}}}$
(3.32)

The first invariant ${\textstyle I_{1}}$ is the trace of the strain tensor:

 ${\displaystyle I_{1}=tr([\varepsilon ])=\varepsilon _{xx}+\varepsilon _{yy}+\varepsilon _{zz}}$
(3.33)

Also, one can always decompose the strain tensor into its volumetric and deviatoric parts ${\textstyle [\varepsilon ]=[\varepsilon _{v}]+[\varepsilon _{d}]}$:

 ${\displaystyle [\varepsilon _{v}]={\begin{bmatrix}{\frac {I_{1}}{3}}&0&0\\0&{\frac {I_{1}}{3}}&0\\0&0&{\frac {I_{1}}{3}}\end{bmatrix}}}$
(3.34)
 ${\displaystyle [\varepsilon _{d}]={\begin{bmatrix}\varepsilon _{xx}-{\frac {I_{1}}{3}}&\varepsilon _{xy}&\varepsilon _{xz}\\\varepsilon _{xy}&\varepsilon _{yy}-{\frac {I_{1}}{3}}&\varepsilon _{yz}\\\varepsilon _{xz}&\varepsilon _{yz}&\varepsilon _{zz}-{\frac {I_{1}}{3}}\end{bmatrix}}}$
(3.35)

From the deviatoric strain tensor ${\textstyle [\varepsilon _{d}]}$ we can calculate ${\textstyle J_{1}}$ and ${\textstyle J_{2}}$:

 ${\displaystyle J_{1}=tr([\varepsilon _{d}])=\varepsilon _{xx}-{\frac {I_{1}}{3}}+\varepsilon _{yy}-{\frac {I_{1}}{3}}+\varepsilon _{zz}-{\frac {I_{1}}{3}}=I_{1}-3{\frac {I_{1}}{3}}=0}$
(3.36)
 ${\displaystyle J_{2}={\frac {1}{2}}([\varepsilon _{d}]:[\varepsilon _{d}]-J_{1}^{2})={\frac {1}{2}}([\varepsilon _{d}]:[\varepsilon _{d}])}$ ${\displaystyle ={\frac {1}{3}}\left[\varepsilon _{xx}^{2}+\varepsilon _{yy}^{2}+\varepsilon _{zz}^{2}-\varepsilon _{xx}\varepsilon _{yy}-\varepsilon _{xx}\varepsilon _{zz}-\varepsilon _{yy}\varepsilon _{zz}\right]+\varepsilon _{xy}^{2}+\varepsilon _{xz}^{2}+\varepsilon _{yz}^{2}}$
(3.37)

The Mazars form of equivalent strain (3.30) and the modified von Mises equation (3.31) lead to different damage surfaces as shown in Figures 26c and 26d. However, in both cases the elastic limit in tension ${\textstyle \sigma _{y_{t}}}$ is considerably lower than the elastic limit in compression ${\textstyle \sigma _{y_{c}}}$, and thus they are convenient for modelling damage progression of quasi-brittle materials such as rocks or concrete.

### 3.2.3 Damage evolution law

Once we have defined the energy norm in the strain space or equivalent strain ${\textstyle \varepsilon _{eq}}$, let us introduce the law ${\textstyle g(r)}$ governing the evolution of the damage variable (3.7).

There are various damage governing laws that can be effectively used to model damage growth in geo-materials. One option, especially suited for simplified analyses, is the linear softening law:

 ${\displaystyle g(r)={\frac {r_{max}}{r_{max}-r_{y}}}\left(1-{\frac {r_{y}}{r}}\right)}$
(3.38)

Such a relation limits the range of the internal state variable ${\textstyle r}$ between the damage threshold and a maximum admissible value ${\textstyle r_{max}}$, and leads to a linear softening branch in the stress-strain curve (Figure 27a).

 (a) Linear softening law. (b) Polynomial softening law. (c) Exponential softening law. Figure 27: Generic unidimensional stress-strain curves for different softening laws.

Two typical choices to describe the evolution of degradation in quasi-brittle materials are the polynomial law [118]:

 ${\displaystyle g(r)=1-{\frac {1}{1+S(r-r_{y})+R(r-r_{y})^{2}}}}$
(3.39)

and the exponential softening law [112]:

 ${\displaystyle g(r)=1-{\frac {r_{y}(1-R)}{r}}-R\exp\{-S(r-r_{y})\}}$
(3.40)

In the above expressions, the parameter ${\textstyle R}$ is associated to the residual strength of the material, whereas the parameter ${\textstyle S}$${\displaystyle S}$softening slope controls the slope of the softening branch after the peak of the stress-strain curve.

The polynomial equation (3.39) and the exponential softening law (3.40) mainly differ from the linear law (3.38) in that the material preserves some residual strength in the post-peak branch, as shown in the 1D stress-strain curves of Figures 27b and 27c.

An alternative exponential softening model was proposed in [59]:

 ${\displaystyle g(r)=1-{\frac {r_{y}}{r}}\exp \left\{A_{f}\left(1-{\frac {r}{r_{y}}}\right)\right\}}$
(3.41)

The parameter ${\textstyle A_{f}}$ is obtained from the following expression:

 ${\displaystyle A_{f}=\left({\frac {G_{f}}{l_{f}(r_{y})^{2}}}-{\frac {1}{2}}\right)^{-1}\geq 0}$
(3.42)

where ${\textstyle G_{f}}$${\displaystyle G_{f}}$ is the specific fracture energy per unit area, ${\textstyle l_{f}}$ is the characteristic length for the fractured domain, usually taken as the characteristic length of the finite elements, and ${\textstyle r_{y}}$ is the aforementioned damage threshold which, for quasi-brittle materials, can be estimated from the tensile strength ${\textstyle \sigma _{y_{t}}}$:

 ${\displaystyle r_{y}={\frac {\sigma _{y_{t}}}{\sqrt {E}}}}$
(3.43)

Another popular damage evolution law specifically designed for concrete was proposed by Mazars [116,112]. He divided the damage variable in two components: ${\textstyle \mathbb {d} _{c}}$ and ${\textstyle \mathbb {d} _{t}}$, which are computed from the equivalent strain of Eq (3.30), but using two different damage laws: ${\textstyle g_{c}}$ and ${\textstyle g_{t}}$. Function ${\textstyle g_{c}}$ is characterized by the uniaxial compressive test, while ${\textstyle g_{t}}$ corresponds to the uniaxial tensile test. Both functions were proposed in the same exponential form of (3.40):

 ${\displaystyle g_{c}(r)=1-{\frac {r_{y}(1-R_{c})}{r}}-R_{c}\exp\{-S_{c}(r-r_{y})\}}$
(3.44)
 ${\displaystyle g_{t}(r)=1-{\frac {r_{y}(1-R_{t})}{r}}-R_{t}\exp\{-S_{t}(r-r_{y})\}}$
(3.45)

The damage evolution law ${\textstyle g(r)}$ results from the combination of the two parts:

 ${\displaystyle g(r)=\alpha _{t}^{\beta }g_{t}(r)+(1-\alpha _{t})^{\beta }g_{c}(r)}$
(3.46)

The coefficient ${\textstyle \alpha _{t}}$ ranges from ${\textstyle 0}$ to ${\textstyle 1}$ and takes into account the character of the stress state. It is evaluated from:

 ${\displaystyle \alpha _{t}=\sum _{i=1}^{3}{\frac {\varepsilon _{ti}\langle \varepsilon _{i}\rangle }{\varepsilon _{eq}^{2}}}}$
(3.47)

where ${\textstyle \varepsilon _{ti}}$ are the principal strains due to positive effective stresses (or tensile stresses).

The parameter ${\textstyle \beta }$ in (3.46) was equal to ${\textstyle 1}$ in the original version of the model [116]. For values higher than ${\textstyle 1}$, ${\textstyle \beta }$ allows to slow down the evolution of damage under shear loading (when principal stresses have different sign).

The definition of ${\textstyle \varepsilon _{ti}}$ tells us that if all principal stresses are negative then ${\textstyle \alpha _{t}=0}$ and ${\textstyle \mathbb {d} =\mathbb {d} _{c}=g_{c}(r)}$, corresponding to a “purely compressive” stress state. On the other hand if all principal stresses are positive, i.e. in a “purely tensile” stress states, we have ${\textstyle \alpha _{t}=1}$ and ${\textstyle \mathbb {d} =\mathbb {d} _{t}=g_{t}(r)}$.

## 3.3 Local and Non-local Damage Models

In the last section we have presented the fundamental theory of the isotropic damage models with a unique scalar variable ${\textstyle \mathbb {d} }$. Although these models are relatively simple, they are often used to model the failure of concrete, rocks and other quasi-brittle materials that show important strain localization. If the damage parameter depends only on the strain state at the point under consideration, the numerical simulations exhibit a pathological mesh dependence, and physically unrealistic results are obtained.

This is the typical behaviour of the so-called local damage models, which are not able to properly describe both the thickness of localization and distance between them. They suffer from mesh sensitivity (for size and alignment) and produce unreliable results. The strains concentrate in one element wide zones and the computed force-displacement curves are mesh-dependent. The reason is that differential equations of motion change their type (from elliptic to hyperbolic in static problems) and the rate boundary value problem becomes ill-posed [114].

Classical constitutive models require an extension in the form of a characteristic length to properly model the thickness of localized zones. Such extension can be done with micro-polar [64,119], strain gradient [65,120,121], viscous [66,122,123] and non-local terms [124,70,67,125]. In this work we have developed the latter approach.

First of all, we present the problems derived from the strain localization phenomenon, using a simple uniaxial case as example. Following that, the basic concepts of the implemented non-local damage model are introduced.

### 3.3.1 Strain localization phenomenon

The idea of modelling damaged concrete and other quasi-brittle materials as strain-softening continua, was not immediately accepted by all the scientific community. Actually, most of the early analyses were not truly objective and, upon mesh refinement, their results would not converge to a robust solution. Let us explain the nature of the problem by means of a unidimensional example.

Consider a straight bar with a constant cross section ${\textstyle A}$ and a total length ${\textstyle L_{0}}$ under uniaxial tension (Figure 28). The material is assumed to obey a simple stress-strain law with linear elasticity up to the peak stress ${\textstyle \sigma _{y}}$, followed by linear softening (Figure 29). If the bar is loaded in tension by an applied displacement ${\textstyle u}$ at its right extreme, the response remains linear elastic up to ${\textstyle u_{y}=L_{0}\varepsilon _{y}}$, instant at which the force transmitted by the bar (reaction at the support) attains its maximum value ${\textstyle F_{y}=\sigma _{y}S}$. After that, the resistance of the bar starts decreasing until the strain reaches ${\textstyle \varepsilon _{f}}$ and the transmitted stress completely disappears.

 Figure 28: Bar under uniaxial tension.
 Figure 29: Stress-strain diagram of the linear softening law.

The equilibrium equations at a section in the middle imply that the axial forces are constant along the bar and so the stress profile must remain also uniform (Figure 30).

 Figure 30: Axial force acting along the bar.

However, at a given stress level ${\textstyle \sigma _{1}}$ between ${\textstyle 0}$ and ${\textstyle \sigma _{y}}$, there are actually two values of strain, ${\textstyle \varepsilon _{1,e}}$ and ${\textstyle \varepsilon _{1,s}}$, that satisfy the constitutive equation (Figure 31). This is quite straightforward if one considers that any cross section can either be damaged, or undamaged. Indeed, an undamaged section is on the elastic branch with ${\textstyle \sigma _{1}=E\varepsilon _{1,e}}$, whereas a damaged one falls in the softening branch with ${\textstyle \sigma _{1}=(1-\mathbb {d} _{1})E\varepsilon _{1,s}}$.

 Figure 31: Possible strain values corresponding to the same stress level.

Thereby, the strain profile along the bar does not have to be necessarily uniform. In fact, any piecewise constant strain distribution that jumps between the two possible strain values represents a valid solution.

 Figure 32: Piecewise constant strain profile along the bar.

In Figure 32 we have denoted by ${\textstyle L_{s}}$ the cumulative length of the softening regions and by ${\textstyle L_{e}=L_{0}-L_{s}}$ the cumulative length of the elastic regions. When the stress is completely relaxed, the strain in the elastic region is ${\textstyle \varepsilon _{e}=0}$ and the strain in the softening branch is ${\textstyle \varepsilon _{s}=\varepsilon _{f}}$. Thus the total elongation of the bar in this case is ${\textstyle u_{f}=L_{e}\varepsilon _{e}+L_{s}\varepsilon _{s}=L_{s}\varepsilon _{f}}$.