A historical review of theory of deformable porous media, ﬁnite element modeling and solution algorithms

This work provides a historical perspective of the development of the theory of deformable porous media, along with its ﬁnite element modeling and algorithms to solve the resulting system of equations.


Introduction
The discovery of fundamental mechanical effects in saturated deformable porous solids, and the formulation of the first porous media theories, are mainly due to two distinguished professors at the Technische Hochschule of Vienna, namely Paul Fillunger and Karl von Terzaghi. [1] was the first to state that the pore liquid pressure does not have any influence upon the strength of the porous solid. Besides this comment, he remarked that the pore water pressure does not affect the material behaviour of the porous solid at all. [2] started the development of the idea of 'effective stress' within the framework of the treatment of the consolidation problem for clay layers. According to [3], the formalism of the idea of effective stress was finally given by [4]. The basic idea of the principle of effective stress is that the 'effective' stress responsible for causing soil deformation is the excess of the total stress over the pore fluid pressure. This implies that the linear momentum balance for poroelasticity can be obtained by simply replacing the total stress in the linear momentum balance for elasticity by the effective stress. In 1941, a Belgian physicist named Maurice Anthony Biot extended Terzaghi's one -dimensional theory to the three -dimensional case and further introduced a parameter representative of the degree of saturation of the fluid inside the pores of the solid (see [5]). He assumed the following properties of soils • isotropy of the material, • reversibility of stress-strain relations under final equilibrium conditions, • linearity of stress-strain relations, • small strains, • the water contained in the pores is incompressible, • the water may contain air bubbles, • the water flows through the porous skeleton according to Darcy's law With E, G, ν representing the Young's modulus, shear modulus, and Poisson's ratio respectively, θ the increment of water volume per unit volume of soil, σ the increment of water pressure, T S A historical review of theory of deformable porous media, finite element modeling and solution algorithms the stress, and u the soil displacement, he introduced coefficients H, R and α and arrived at the following relations where E L is the well-known linearized Green strain tensor given by The corresponding stresses must satisfy the equilibrium conditions of a stress field (inertia forces and body forces are neglected): The relation (1) along with the conditions (3) result in is a Lamé parameter. In order to have a complete set of equations for determining u and σ, one more equation is necessary. This equation can be obtained by introducing Darcy's law governing the flow of water in a porous medium: where w F represents the volume of water flowing per second and unit area through the faces of an elementary cube and k F denotes the permeability coefficient of the soil. [6] observed in tests with natural sand, the proportionality of the total volume of water running through the sand and the loss of pressure. The relationship he arrived at is called the 'Darcy's law' and is repeatedly used as the equation that governs the flow of fluid in porous media. Furthermore, under the assumption of the incompressibility of water, the rate of water content of an element of soil must be equal to the volume of water entering per second through the surface of the elements.
Combining (2), (5) and (6), (4) and (7) are the basic relations of Biot's theory of consolidation, where (4) describes the settlement of the soil, and (7) describes the change in the water pressure. Furthermore, it is important to note that these two relations are coupled. [7] relaxed the assumption of incompressibility of the pore fluid and isotropy of the pore skeleton in order to render a set of equations applicable to an anisotropic porous solid containing a viscous compressible fluid. The resulting set of equations for the specific case of isotropic porous solid had the same number (three) of coefficients to be determined as the original equations (4) and (7). [8] reworked the set of equations with a reduced number (two) of coefficients and further devised tests to determine those coefficients. The interested reader is refered to the review article of [9] for a brief description of these tests. The authors summarized the equations for the isotropic case (inertia forces and body forces are neglected): A historical review of theory of deformable porous media, finite element modeling and solution algorithms ∂ζ ∂t + div w F = 0 (9) where ζ and σ are the increment in fluid content and total stress respectively and are given by where φ is the porosity, u F is the fluid displacement and p is the pore fluid pressure. The interested reader is refered to [10] for a better understanding of the fluid mass increment. The variables M and α would be refered to as the Biot modulus and the Biot constant respectively. The Biot constant is representative of the degree of saturation of pore fluid inside the pores with α = 1 being the fully saturated state and α < 1 being the partially saturated state.
The following three tests were devised to determine M and α.
• The jacketed compressibility test: A specimen of the material is enclosed in a thin impermeable jacket and then subjected to an external fluid pressure p . In this test the entire pressure of the fluid is transmitted to the solid portions of the surfaces of the specimen. To insure constant internal fluid pressure, the inside of the jacket may be made to communicate with the atmosphere through a tube. The dilatation e of the specimen is measured and a coefficient of jacketed compressibility κ is determined by • The unjacketed compressibility test: A sample of the material is immersed in a fluid to which is applied a pressure p . In this case the pressure acts both on the solid and fluid portions of the surfaces of the specimen. When the fluid pressure has penetrated the pores completely, the dilatation of the sample is then measured and an unjacketed compressibility coefficient δ is determined by • The unjacketed coefficient of fluid content test: A unit volume of porous material containing fluid is placed within a closed chamber which has been filled with fluid. Fluid is then injected into the chamber under pressure and the volume of injected fluid is measured. The volume of fluid injected per unit pressure will be the sum of the solid compressibility δ, the volume of fluid γ which has entered the pores 1 , and a fixed quantity depending upon the elastic properties of the chamber and the fluid. The porous material is then removed from the chamber and its volume replaced by fluid. The volume of fluid injected per unit pressure is again measured and in this case will be the sum of the same fixed quantity as in the previous measurement, and the fluid compressibility c representing the new unit volume occupied by the fluid. Therefore the difference between the volumes injected with and without the porous material in the chamber will be given by If the fluid compressibility c is then measured independently, the coefficient of fluid content γ may then be determined.
A historical review of theory of deformable porous media, finite element modeling and solution algorithms [8] expressed M and α in terms of κ, δ and γ for the case of a homogeneous, isotropic and elastically linear porous solid as follows The increment in fluid content ξ was stated to be a work conjugate to increment in fluid pressure p. Following the work of [11] and [12], the Biot modulus and Biot constant were established as where K and K s represent the effective and grain bulk moduli respectively. We shall see in Chapter 3 that if the Biot modulus is assumed to be time-invariant, then (9) corresponds to the mass conservation equation for linearized slightly compressible single phase flow. [13] reworked the Biot thery in terms of parameters that were representative of the skeletal response at t = 0 called the 'undrained response' and t → ∞ called the 'drained response'. The undrained response was characterized by the fluid mass increment being zero and the drained response was characterized by the pore pressure being zero. The stress-strain relations under drained conditions are identical to the ones of non -porous media, provided they are expressed in terms of the effective stress. [14] introduced a parameter called Skempton coefficient that characterizes the instantaneous response (t = 0) of the porous skeleton. It is the ratio between the induced pore pressure and mean stress under a confining pressure in undrained conditions as follows It is clear from the definition (15) that 0 ≤ B ≤ 1. [15] arrived at the following relation for the Skempton coefficient Ks ) The undrained parameters are obtained by setting the fluid increment to zero (ξ = 0). From (11), that results in from which the undrained bulk modulus K u can be obtained as From (13), (14) and (17), the undrained bulk modulus can be written as It is easy to see that in case of a porous solid saturated with compressible pore fluid i.e. α = 1, c > 0, the elastic skeleton is compressible under undrained conditions and in case of a porous solid saturated with incompressible pore fluid i.e. α = 1, c = 0, the elastic skeleton is incompressible under undrained conditions Further, from (16), it is easy to see that which is the same with the case of drained deformation. This means that shear modulus is the same in drained and undrained deformations.

From Navier-Stokes equations to Darcy's law
Although the Darcy's law is the most commonly used equation for flow in porous media, fluid mechanics theory suggests that fluid flow is in fact governed by the Navier-Stokes equations of balance of momentum of the fluid. The key component in arriving at a connection between the Darcy's law and the Navier-Stokes equations is the 'averaging theorem'. The theorem is based upon the well-known Reynolds transport theorem and relates the volume average of the spatial derivative to the spatial derivative of the volume average of any quantity: scalar, vector or tensor.
The theorem was first presented in [16] and later elucidated upon by [17]. [18] used the averaging technique to show that the Navier-Stokes equations for flow in a deforming anisotropic porous medium reduced to the Darcy's law in the limit of slow flow.

Finite element modeling
The first successful attempt in applying the finite element method to consolidation problem was reportedly made by [19] for the case of elastic porous solid saturated with incompressible pore fluid. [20] reworked the approach with the assumption of incompressibility of pore fluid and solid grains being relaxed. They used four-noded quadrilateral interpolation functions for both the pore pressure and the displacement variables but with additional incompatible interpolation functions for the latter. [21] presented the idea of 'incompatible mode elements' to resolve the problem of shear locking in thin beams. Shear locking refers to the large, unphysical shear strains that arise in thin beams due to the use of standard four-noded linear quadrilateral plane stress elements. [22] established limits on the value of the time-marching coefficient for the unconditional stability of the fully discrete formulation of the equations of consolidation. [23] studied the response of a fully discrete formulation of the equations of two-dimensional consolidation of an elastic porous solid saturated with incompressible pore fluid with six-noded triangular elements and eight-noded quadrilateral elements. They reported that schemes that used equal order interpolation for pore pressure and displacement gave initial errors in pore pressure which do not dissipate in time and are oscillatory in space. [24] established lower bounds on the time step size below which the fully discrete formulation yields an inaccurate pore pressure distribution with violent oscillations.
A historical review of theory of deformable porous media, finite element modeling and solution algorithms

The Stokes' problem under undrained conditions with incompressible fluid and solid constituents
From (4) and (19), it is clear that the elasticity problem under undrained conditions for the case α = 1, c = 0, along with appropriate boundary conditions on u and p is in fact the Stokes' problem given by The work of [23] clearly suggested that the use of equal order finite element interpolation for p and u led to oscillatory solutions for p to the Stokes' problem. The reason lied in the failure to satisfy the well-known Ladyzenskaja-Babuska-Brezzi (LBB) condition (see [25]) which poses a constraint on the finite element spaces used for p and u for the solution to be stable. [26] spawned a class of methods called 'stabilized methods' meant to circumvent the deficiencies of equal order interpolation for p and u in obtaining a stable solution to the problem. [27] studied the problem in the context of a saddle-point systems and proposed an extension to the well-known patch test as an alternative to the mathematically rigorous LBB condition. The interested reader is refered to [28] for an understanding of saddle-point problems.  [30]) and proposed a post-processing technique that improves the order of convergence of pore pressure by one. [31] reformulated the problem as a least-squares minimization problem with lowest order Taylor-Hood spaces for pressure and displacement variables and compatible Raviart-Thomas spaces for fluid velocity and stress variables. The interested reader is refered to Chapter III of [32] for a description of various mixed finite element spaces. The reader shall also find an excellent rendition of the treatment of the Stokes' problem in Chapter VI. [33] presented a formulation that combats the nonphysical pressure oscillations using Raviart-Thomas-Nedelec spaces for the flow variables and a discontinuous Galerkin space for the displacement variable. Raviart-Thomas-Nedelec spaces are the three-dimensional analogues of the Raviart-Thomas spaces on rectangles. The interested reader is refered to [34] for an understanding of Discontinuous Galerkin methods. [35] presented a formulation that addresses the problem of oscillatory pore pressures using discontinuous Galerkin elements on areas with high pressure gradients and continuous Galerkin elements on the remaining domain.

Mixed formulation for flow with compressible fluid constituent
We have seen in the previous module that LBB stable mixed finite element spaces are a popular choice for displacement and pressure in the Stokes' equations. But, we have seen in (18) that if the pore fluid compressibility is strictly positive, the incompressibility constraint under undrained conditions is not operative and hence the Stokes' problem for displacement and pressure does not arise. Now, we rewrite the equation of mass conservation (9) for linearized slightly compressible single phase flow along with the Darcy law (10) for steady state conditions ( ∂ξ ∂t = 0) in the presence of gravity as div w F = 0 where ρ F is the fluid density and g is the gravity vector. The above system represents a saddle point problem for w F and p. Thus, although we can avoid the necessity of the use of mixed finite element spaces for displacement and pressure by maintaining a strictly positive pore fluid compressibility, the above argument suggests a need, instead, for mixed finite element spaces for the fluid velocity and pressure under steady state conditions in the presence of gravity.

Solution algorithms
The early works of [19] and [20] along with later renditions of [36][37][38] are refered to as fully coupled schemes where the pressures, fluid velocities and displacements are solved simultaneously. The fully coupled approach, although unconditionally stable, requires careful implementation with substantial local memory requirements and specialized linear solvers. On the other hand, the works of [39][40][41][42] are refered to as staggered solution schemes in which an operator splitting strategy is used to split the coupled problem into well-posed flow and mechanics subproblems which are then solved sequentially. In the mid 1970s, the efforts of Ted Belytschko at Northwestern University, Thomas Hughes at California Institute of Technology and Carlos Felippa at Lockheed Palo Alto Research Laboratories would spawn the class of partitioned solution algorithms to coupled dynamical problems designed to take advantage of increasing modularity in commercial finite element codes. An excellent review of these partitioned solution procedures is given in [43]. The interested reader is also refered to [44] for an understanding of the genesis of the idea of fractional steps. Carefully crafted splitting strategies would lend solution accuracies on par with the fully coupled approaches.
Previous attempts at incorporating non-matching grids for flow and mechanics include the works of [45][46][47][48]. [45] reformulated a staggered solution algorithm as a special case of a fully coupled scheme and implemented the algorithm on overlapping nonmatching rectilinear grids, but avoided three dimensional intersection calculations, instead evaluating the displacement-pressure coupling submatrices using a midpoint integration rule. [46] implemented a procedure in which a saddlepoint system with mortar spaces on nonmatching interfaces of a decomposed geomechanics domain is solved by applying a balancing Neumann-Neumann preconditioner. The procedure involves subdomain to mortar and mortar to subdomain projections, Lagrange multiplier solve and computationally expensive parallel subdomain solves at each time step. [47] implemented a procedure for nested grids in which coarse scale basis functions for the poromechanical solve are obtained in terms of fine scale basis functions by solving local equilibrium problems on each coarse scale poromechanical element. These coarse scale basis functions are then used to construct prolongation and restriction operators, which are then employed to construct a two-stage preconditioner for the coarse scale poromechanical solve. [48] used computational geometry to treat finite elements as three dimensional objects, thereby computing projection operators using intersection volumes. The micromechanical analyses for the case of anisotropic poroelasticity [49][50][51] revealed that the modification to the stress applied to the porous solid due to the presence of pore fluid pressure is not hydrostatic, as it is in the case of isotropic poroelasticity [12]. Further, unlike in case of isotropic poroelasticity where the solid-fluid coupling parameter is a scalar [3,5,11,12], the coupling parameter in case of anisotropic poroelasticity is a tensor [7,10,52].
Some of the earliest work on multiphase flow in deformable porous media was performed in [53][54][55][56]. The need for sophisticated algorithms to solve the coupled system of equations was then felt, which led to a lot of contributions to the understanding of solvers for such systems . Although the bulk of the aformentioned work is devoted to the implementation of nonlinear solvers, due diligence is also given to arriving at thermodynamically consistent set of equations that couple the deformation and multiphase flow effects without violating first principles [80][81][82][83][84][85].
The use of staggered solution algorithms has found a lot of popularity amongst the coupled flow and geomechanics community, and the most popular among the splitting techniques is the fixed stress split strategy. It decouples the flow and geomechanics equations by imposing a algorithms Algorithm 1: Popular fixed stress split algorithm for multiphase geomechanics / * Pre-processing start * / constraint on the flow solve, and then solves the flow problem followed by the geomechanics problem in repeated iterations (see Algorithm 1) until a certain convergence criterion is met at each time step [42,52,[86][87][88][89][90][91][92][93][94][95][96][97][98][99][100][101]. While the work of [42] is appreciated as the benchmark for the fixed stress split technique, there have been further important extensions to the algorithm like the multirate method [102] and the multiscale method [48]. An important piece of work in the realm of theoretical convergence analysis of the fixed stress split technique is provided in [103]. The works of [104][105][106][107][108] provide a linear algebra point of view to the solution of the system of equations using the decoupling technique.