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. 1 Introduction The discovery of fundamental mechanical effects in saturated deformable porous solids, and the formulation of the ﬁrst 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 ﬁrst to state that the pore liquid pressure does not have any inﬂuence 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 ﬁnally 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 ﬂuid 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 ﬂuid 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 ﬁnal 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 ﬂows through the porous skeleton according to Darcy’s law WithE,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, TS A historical review of theory of deformable porous media, ﬁnite element modeling and solution algorithms the stress, and uthe soil displacement, he introduced coefﬁcients H,Rand and arrived at the following relations TS= 2G EL+ 12(divu)I  I (1) = divu+1 R H (2) where ELis the well-known linearized Green strain tensor given by EL=1 2(grad u+ (grad u)T) The corresponding stresses must satisfy the equilibrium conditions of a stress ﬁeld (inertia forces and body forces are neglected): divTS=0 (3) The relation (1) along with the conditions (3) result in Gdiv grad u+ (+G) grad div u grad=0 (4) whereE (1+)(12)is a Lamé parameter. In order to have a complete set of equations for determining uand, one more equation is necessary. This equation can be obtained by introducing Darcy’s law governing the ﬂow of water in a porous medium: wF=kFgrad (5) where wFrepresents the volume of water ﬂowing per second and unit area through the faces of an elementary cube and kFdenotes the permeability coefﬁcient 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 ﬂow of ﬂuid 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. @ @t=divwF(6) Combining (2), (5) and (6), kFdiv grad= @(ELI) @t+1 R H@ @t(7) (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 ﬂuid and isotropy of the pore skeleton in order to render a set of equations applicable to an anisotropic porous solid containing a viscous compressible ﬂuid. The resulting set of equations for the speciﬁc case of isotropic porous solid had the same number (three) of coefﬁcients to be determined as the original equations (4)and (7). [8] reworked the set of equations with a reduced number (two) of coefﬁcients and further devised tests to determine those coefﬁcients. 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): div (2GEL+(divu)I pI) =0 (8) 2 A historical review of theory of deformable porous media, ﬁnite element modeling and solution algorithms @ @t+ div wF= 0 (9) wF=kFgradp (10) whereandare the increment in ﬂuid content and total stress respectively and are given by =1 Mp+ divu (divudivuF) (11) = 2GEL+(divu)I pI (12) where is the porosity, uFis the ﬂuid displacement and pis the pore ﬂuid pressure. The interested reader is refered to [10] for a better understanding of the ﬂuid 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 ﬂuid inside the pores with = 1being the fully saturated state and <1being the partially saturated state. The following three tests were devised to determine Mand . •The jacketed compressibility test: A specimen of the material is enclosed in a thin imper- meable jacket and then subjected to an external ﬂuid pressure p0. In this test the entire pressure of the ﬂuid is transmitted to the solid portions of the surfaces of the specimen. To insure constant internal ﬂuid pressure, the inside of the jacket may be made to communicate with the atmosphere through a tube. The dilatation eof the specimen is measured and a coefﬁcient of jacketed compressibility is determined by =e p0 •The unjacketed compressibility test: A sample of the material is immersed in a ﬂuid to which is applied a pressure p0. In this case the pressure acts both on the solid and ﬂuid portions of the surfaces of the specimen. When the ﬂuid pressure has penetrated the pores completely, the dilatation of the sample is then measured and an unjacketed compressibility coefﬁcientis determined by =e p0 •The unjacketed coefﬁcient of ﬂuid content test: A unit volume of porous material containing ﬂuid is placed within a closed chamber which has been ﬁlled with ﬂuid. Fluid is then injected into the chamber under pressure and the volume of injected ﬂuid is measured. The volume of ﬂuid injected per unit pressure will be the sum of the solid compressibility , the volume of ﬂuid which has entered the pores1, and a ﬁxed quantity depending upon the elastic properties of the chamber and the ﬂuid. The porous material is then removed from the chamber and its volume replaced by ﬂuid. The volume of ﬂuid injected per unit pressure is again measured and in this case will be the sum of the same ﬁxed quantity as in the previous measurement, and the ﬂuid compressibility crepresenting the new unit volume occupied by the ﬂuid. Therefore the difference between the volumes injected with and without the porous material in the chamber will be given by V=+ c If the ﬂuid compressibility cis then measured independently, the coefﬁcient of ﬂuid content may then be determined. 1Withp0refering to the pressure being applied to the chamber, the coefﬁcient of ﬂuid content was deﬁned as = p0 3 A historical review of theory of deformable porous media, ﬁnite element modeling and solution algorithms [8] expressed Mand in terms of,and for the case of a homogeneous, isotropic and elastically linear porous solid as follows M=1 +2  = 1  The increment in ﬂuid content was stated to be a work conjugate to increment in ﬂuid pressure p. Following the work of [11] and [12], the Biot modulus and Biot constant were established as M=1 c +  Ks(13) = 1K Ks(14) whereKandKsrepresent 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 ﬂow. [13] reworked the Biot thery in terms of parameters that were representative of the skeletal response at t= 0called the ‘undrained response’ and t!1 called the ‘drained response’. The undrained response was characterized by the ﬂuid 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 coefﬁcient 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 conﬁning pressure in undrained conditions as follows B=p j=0 (15) It is clear from the deﬁnition (15) that0B1. [15] arrived at the following relation for the Skempton coefﬁcient B=1 K1 Ks 1 K1 Ks+ (c1 Ks) The undrained parameters are obtained by setting the ﬂuid increment to zero ( = 0). From (11), that results in p= Mdivu which when substituted in (12) results in = 2GEL+ (+ 2M) divu I (16) from which the undrained bulk modulus Kucan be obtained as Ku=1 3tr() tr(EL)=2G 3++ 2M=K+ 2M (17) From (13), (14) and (17), the undrained bulk modulus can be written as Ku=K 1 + 2 c K + (  )(1 ) 4 A historical review of theory of deformable porous media, ﬁnite element modeling and solution algorithms It is easy to see that in case of a porous solid saturated with compressible pore ﬂuid i.e. = 1, c>0, the elastic skeleton is compressible under undrained conditions Ku=K 1 +1 c K tr(EL) = div u=1 3tr() Ku6= 0 (18) and in case of a porous solid saturated with incompressible pore ﬂuid i.e. = 1,c= 0, the elastic skeleton is incompressible under undrained conditions Ku=1 tr(EL) = div u=1 3tr() Ku= 0 (19) Further, from (16), it is easy to see that ij= 2GEL ij (ifi6=j) which is the same with the case of drained deformation. This means that shear modulus is the same in drained and undrained deformations. 1.1 From Navier-Stokes equations to Darcy’s law Although the Darcy’s law is the most commonly used equation for ﬂow in porous media, ﬂuid mechanics theory suggests that ﬂuid ﬂow is in fact governed by the Navier-Stokes equations of balance of momentum of the ﬂuid. 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 ﬁrst presented in [16] and later elucidated upon by [17]. [18] used the averaging technique to show that the Navier-Stokes equations for ﬂow in a deforming anisotropic porous medium reduced to the Darcy’s law in the limit of slow ﬂow. 2 Finite element modeling The ﬁrst successful attempt in applying the ﬁnite element method to consolidation problem was reportedly made by [19] for the case of elastic porous solid saturated with incompressible pore ﬂuid. [20] reworked the approach with the assumption of incompressibility of pore ﬂuid 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 coefﬁcient 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 ﬂuid 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. 5 A historical review of theory of deformable porous media, ﬁnite element modeling and solution algorithms 2.1 The Stokes’ problem under undrained conditions with incompressible ﬂuid 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 uandpis in fact the Stokes’ problem given by Gdiv grad ugradp=0 divu= 0 The work of [23] clearly suggested that the use of equal order ﬁnite element interpolation for p anduled to oscillatory solutions for pto the Stokes’ problem. The reason lied in the failure to satisfy the well-known Lady ˘zenskaja-Babu ˘ska-Brezzi (LBB) condition (see [25]) which poses a constraint on the ﬁnite element spaces used for pandufor the solution to be stable. [26] spawned a class of methods called ‘stabilized methods’ meant to circumvent the deﬁciencies of equal order interpolation for panduin 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. The patch test for ﬁnite elements was ﬁrst presented by the British engineer Bruce Irons at a conference on ‘Matrix Methods in Structural Mechanics’ held at Wright-Patterson Air Force Base, Ohio, USA on 26-28 October, 1965. The report can be found at http://contrails.iit.edu/items/show/8575 . [29] suggested the use of LBB stable Taylor-Hood elements (see [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 ﬂuid velocity and stress variables. The interested reader is refered to Chapter III of [32] for a description of various mixed ﬁnite element spaces. The reader shall also ﬁnd 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 ﬂow 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. 2.2 Mixed formulation for ﬂow with compressible ﬂuid constituent We have seen in the previous module that LBB stable mixed ﬁnite element spaces are a popular choice for displacement and pressure in the Stokes’ equations. But, we have seen in (18) that if the pore ﬂuid 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 ﬂow along with the Darcy law (10) for steady state conditions (@ @t= 0) in the presence of gravity as divwF= 0 kF1wF+ gradp=Fg whereFis the ﬂuid density and gis the gravity vector. The above system represents a saddle point problem for wFandp. Thus, although we can avoid the necessity of the use of mixed ﬁnite element spaces for displacement and pressure by maintaining a strictly positive pore ﬂuid compressibility, 6 A historical review of theory of deformable porous media, ﬁnite element modeling and solution algorithms the above argument suggests a need, instead, for mixed ﬁnite element spaces for the ﬂuid velocity and pressure under steady state conditions in the presence of gravity. 3 Solution algorithms The early works of [19] and [20] along with later renditions of [36 –38] are refered to as fully coupled schemes where the pressures, ﬂuid 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 –42] are refered to as staggered solution schemes in which an operator splitting strategy is used to split the coupled problem into well-posed ﬂow 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 ﬁnite 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 ﬂow and mechanics include the works of [45 –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 saddle- point 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 ﬁne 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 ﬁnite elements as three dimensional objects, thereby computing projection operators using intersection volumes. The micromechanical analyses for the case of anisotropic poroelasticity [49 –51] revealed that the modiﬁcation to the stress applied to the porous solid due to the presence of pore ﬂuid pressure is not hydrostatic, as it is in the case of isotropic poroelasticity [12]. Further, unlike in case of isotropic poroelasticity where the solid-ﬂuid 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 ﬂow in deformable porous media was performed in [53 –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 [57 –79]. 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 ﬂow effects without violating ﬁrst principles [80–85]. The use of staggered solution algorithms has found a lot of popularity amongst the coupled ﬂow and geomechanics community, and the most popular among the splitting techniques is the ﬁxed stress split strategy. It decouples the ﬂow and geomechanics equations by imposing a 7 A historical review of theory of deformable porous media, ﬁnite element modeling and solution algorithms Algorithm 1: Popular ﬁxed stress split algorithm for multiphase geomechanics /*Pre-processing start */ 1Initialize pore pressures; 2Solve geomechanics to equilibriate stresses; 3t 0; /*Pre-processing end */ 4while t<Tdo /*Loop over time steps */ 5 while Not converged do /*Loop over coupling iterations for the staggered solution algorithm with pre-defined convergence criterion. The criterion can be based off any metric that reflects the spatio-temporal evolution of the system */ 6 while Not converged do /*Loop over Newton iterations for flow */ 7 Solve multiphase ﬂow keeping rate of volumetric stress ﬁxed; 8 Update variables; 9 while Not converged do /*Loop over Newton iterations for geomechanics */ 10 Solve geomechanics; 11 Update variables; 12 t t+dt; constraint on the ﬂow solve, and then solves the ﬂow 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 –101]. While the work of [42] is appreciated as the benchmark for the ﬁxed 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 ﬁxed stress split technique is provided in [103]. The works of [104 –108] provide a linear algebra point of view to the solution of the system of equations using the decoupling technique.