A topology-motivated mixed finite element method for dynamic response of porous media

A topology-motivated mixed finite element method for dynamic response of porous media

Z. Lotfian and M. V. Sivaselvan

In this paper, we propose a numerical method for computing solutions to Biot’s fully dynamic model of incompressible saturated porous media [1]. Our spatial discretization scheme is based on the three-field formulation (--) and the coupling of a lowest order Raviart-Thomas mixed element [2] for fluid variable fields (, ) and a nodal Galerkin finite element for skeleton variable field (). These mixed spaces are constructed based on the natural topology of the variables; hence, are physically compatible and able to exactly model the kind of continuity which is expected. The method automatically satisfies the well known LBB (inf-sup) stability condition and avoids locking that usually occurs in the numerical computations in the incompressible limit and very low hydraulic conductivity. In contrast to the majority of approaches, our three-field formulation can fully capture dynamic behavior of porous media even in high frequency loading phenomena with considerable fluid acceleration such as liquefaction and biomechanics of porous tissues under rapid external loading. Moreover, we address the importance of consistent initial conditions for poroelasticity equations with the incompressibility constraint, which represent a system of differential algebraic equations. The energy balance equation is derived for the full porous medium and used to assess the stability and accuracy of our time integration. To highlight the capabilities of our method, a variety of numerical studies are provided including verification with analytical and boundary element solutions, wave propagation analyses, hydraulic conductivity effects on damping and frequency content, energy balance analyses, mass lumping considerations, effects of mesh pattern and size, and stability analyses. We also explain some discrepancies commonly found in dynamic poroelasticity results in the literature.

A topology-motivated mixed finite element method for dynamic response of porous media

Z. Lotfian and M. V. Sivaselvan

Department of Civil, Structural and Environmental Engineering, University at Buffalo, Buffalo, NY 14221, U.S.A.

Received …

KEY WORDS:    dynamic poroelasticity, Raviart-Thomas, mixed finite element, LBB condition, locking

11footnotetext: Correspondence to: zahrasad@buffalo.edu

1 Introduction

Dynamics of porous media are of interest in a variety of fields including mechanics of saturated and partially saturated soil, porous biological materials and petroleum reservoirs. The dynamic response of saturated soil under earthquake and the subsequent soil liquefaction process is subject of many studies [3, 4, 5]. Additionally, the ability to accurately model the behavior of rock/soil and fluids in petroleum reservoirs is essential for maximizing oil recovery and assessing the performance of processes such as wellbore stability and subsidence. For this also, coupled modeling of geomechanics and fluid flow is required [6, 7]. Similarly in biomechanics, mechanical response of hydrated soft tissues such as brain, skin, and articular cartilage [8, 9, 10, 11, 12], bones [13], and even individual cells [14, 15] involves coupled response of a skeleton and pore fluid. Furthermore, all the results of dynamic poromechanics can be applied to dynamic thermomechanics based on the mathematical analogy between the formulations [16, 17].

In this paper, we propose a numerical method for dynamic poroelasticity, starting from Biot’s formulation [1]. The porous medium is considered as a solid skeleton saturated with fluid. While the skeleton itself is deformable and elastic, the material that makes up the skeleton and the fluid are modeled as incompressible. This is reasonable for many of the physical problems described above. We also restrict the formulation to linearized kinematics. Alternate formulations based on the Theory of Porous Media [18] have been shown to be equivalent in this realm of incompressible solid and fluid materials and linear kinematics [19]. Our focus is on fully dynamic problems, i.e., inertia effects in both the fluid and the skeleton have been taken into account. We note that this is in contrast to most previous work, where fluid acceleration terms are fully or partially neglected [3, 5, 11, 20].

Some of the specific features of our approach are as follows: (a) In contrast to much of the literature, we use a full three-field formulation in which the skeleton displacement (), pore pressure () and pore-fluid relative velocity () fields are all approximated independently and calculated directly. As a result (see section 3.1), the formulation is applicable to problems with considerable fluid acceleration, for instance in liquefaction analysis [4] and biomechanics of porous tissues under rapid external loading [12, 10]. Moreover, greater accuracy is attained compared to formulations where some fields are computed by post processing techniques as reported in [21, 22]. (b) We employ a topology-motivated spatial discretization based on the lowest order Raviart-Thomas () mixed element [2] to approximate the flow variables and a linear/quadratic standard Galerkin finite element for the skeleton displacement. This mixed approach automatically enforces the precise kind of continuity (local mass conservation) that is required. Additionally, in the incompressible limit (when material is modeled as incompressible and under low hydraulic conductivity) and the rigid skeleton limit, the proposed method satisfies the well known LBB (inf-sup) condition [23, 24], avoiding associated locking effects and the need for any additional stabilization techniques. (c) Taking a differential algebraic equations (DAE) point of view, we demonstrate the importance of consistent initial conditions for poroelasticity equations with the incompressibility constraint. (d) The energy balance equation for the whole porous medium is also derived and used to confirm the stability and accuracy of the Newmark time integration method with particular emphasis on consistent initial condition.

In many problems, such as soil mechanics and biomechanics of soft tissues, the fluid and skeleton material could be reasonably modeled as being incompressible. However when developing a finite element formulation, the incompressibility condition imposes a certain constraint (see section 3.2), which in turn gives rise to stability concerns. The discrete form of the Ladyženskaja-Babuška-Brezzi (LBB) condition provides the necessary and sufficient criterion for a stable finite element formulation. To satisfy this condition, several approaches including stabilized methods [25, 26, 27, 28, 6, 29, 30], discontinuous Galerkin [31, 32, 33, 34, 35, 36, 37], special elements such as Taylor-Hood elements [38] or MINI elements [39], and mixed finite elements [2, 40, 41] have been proposed in the literature. As we discuss in section 3.2, in a finite element formulation for poroelasticity, the LBB condition must be satisfied in two limiting cases — low hydraulic conductivity and rigid skeleton [6, 42]. In this paper, we develop a finite element formulation using a two fold approach to satisfy the LBB condition — (i) We discretize the pressure () and velocity () fields of the pore-fluid flow using Raviart-Thomas (RT) elements [2]. These elements are topology-inspired, associating velocity degrees of freedom (DOF) with element edges and pressure DOF with the element interior, and result in automatically satisfying the LBB condition in the rigid skeleton limit. (ii) We discretize the skeleton displacement field () using standard Galerkin finite elements in such a way that together with the fluid pressure field, the LBB condition in the low hydraulic conductivity limit is satisfied in the same manner as in incompressible elasticity problems [23, 43]. In the context of RT elements, we note other topology-inspired finite elements [44, 45, 39].

Most finite elements of the three-field category for porous media are obtained using standard nodal Galerkin discretization of all three fields [46, 22, 47, 10, 48, 49, 50, 51, 4, 52]. However, these finite elements usually fail to fulfill the LBB condition. As described above, in this paper, we exploit a mixed finite element approach to deal with the challenges that conventional node-based Galerkin methods are not able to handle easily. The main advantage of mixed elements is their ability to exactly model the kind of continuity which is expected i.e., continuity of the normal velocity across element edges which results in the desirable element-wise (local) mass conservation. Mixed elements form a structure that imitates, at the discrete level, all the operators that exist at the continuous level PDEs (such as divergence and gradient). In addition, boundary conditions can be applied more naturally. For these reasons, mixed elements such as RT elements [2], Nédélec elements [40] or a variety of others summarized in [41] are widely used in solid and fluid mechanics problems of incompressible elasticity and fluid flow. However, their application in coupled problems (simultaneous solution) has been limited. Phillips and Wheeler [53, 54] and Lipnikov [55] developed a formulation based on a combination of the lowest order Raviart-Thomas () mixed finite element for the and fields and the standard nodal linear Galerkin finite element for the field for quasistatic poroelasticity problems. We call this element . They elaborated on the solvability and error estimates of this finite element for quasistatic consolidation problems. In this paper, we extend the application of this element to dynamic problems. In doing so, we are able to effectively capture wave propagation effects (section 6.1). We also observe that careful consideration of the finite element mesh is required to obtain the proper relationship between hydraulic conductivity and the damping of the mechanical response (section 6.2). This relationship has not been captured correctly by many past results in the literature. Furthermore, we find that while stability in terms of the LBB condition is achieved for typical ranges of hydraulic conductivity, in the limiting case of very small hydraulic conductivity, the element exhibits instability in the form of checkerboard pattern in the field. For reasons described in section 6.3.2, this effect (also referred to as locking in this context), is more pronounced in dynamic problems than previously observed in quasistatic problems [56]. To resolve this issue, we propose a element which employs a quadratic Galerkin finite element to estimate the space of . The element is shown to pass the numerical inf-sup test and therefore the stability criterion (see section 6.3.2).

This paper is organized as follows. First we summarize Biot’s mathematical model for porous media in section 2. In section 3, we motivate and present our mixed finite element formulation; show how the poroelasticity equations reduce to two classical models in the limiting case, and discuss the associated stability considerations; we develop our spatial discretization scheme based upon the three-field formulation (--) and describe the proposed finite element formulation using RT elements for the - fields and nodal Galerikin element for the field. In section 4, we discuss time integration, in particular addressing the need for consistent initial conditions; we also present the energy balance equation for the full porous medium, which we later use to assess the stability and accuracy of time integration. We present the details of implementation in section 5, including aspects involving unconventional edge degrees of freedom. To demonstrate the performance of this approach, we consider detailed numerical examples in section 6. These include verification with analytical and boundary element solutions, wave propagation analysis, energy balance analyses, mass lumping considerations, effects of mesh pattern and size, hydraulic conductivity effects on damping and frequency, spurious pressure modes and locking, and stability analysis by means of inf-sup tests. We also comment on some discrepancies found in the results in the literature studies.

2 Governing equations

We describe Biot’s equations [1] under the following conditions: (i) fluid and skeleton materials are incompressible; (ii) there is a single fluid phase; (iii) skeleton has a linear elastic behavior; (iv) a linear form of Darcy’s law (i.e,. constant hydraulic conductivity) is used; (v) a linearized kinematics formulation is applied which also results in constant porosity; (vi) the additional apparent massalso called the added mass; in porous medium, when the fluid (skeleton) is moving, a force must be exerted on the skeleton (fluid) to prevent an average displacement of the latter which is defined through a mass coupling coefficient. used in Biot’s equations is neglected; (vii) to focus on the essential features of our formulation, we do not consider fluid sources or sinks, body forces or gravity effects, inclusion of which is straightforward. Accordingly, the linearized governing equations read [57]


These represent the momentum balance of the porous medium, the dynamic extension of Darcy’s law, and the mass balance equation (continuity equation) respectively. The three basic fields are skeleton displacement, , pore fluid pressure, , and fluid Darcy velocity, , as shown in Figure (a)a. Darcy velocity is defined as the relative rate of discharge per unit area of the medium, , where refers to the fluid absolute velocity. , , and dot denote the gradient operator, the divergence operator, and the derivative with respect to time, , respectively. is the fluid density, while is the average density of the medium, with and being solid phase density and porosity respectively. , , and represent the identity matrix, acceleration due to gravity, and the hydraulic conductivity (its associated matrix is assumed to be isotropic).

Figure 3: (a) Saturated porous medium with primary variables. (b) Domain corresponds to the skeleton boundaries is shown on the left while at the same time, corresponds to the fluid boundaries is shown on the right. In each figure boundary is divided into two disjoint subsets which cover the entire boundary of the domain; i.e., and .

The concept of effective stress is used as a principle that controls the constitutive behavior of a porous medium: where the total stress is composed of effective stress of skeleton and the pore pressure . Here, denotes the Biot-Willis coefficient which is commonly taken as 1 for soil material. For an isotropic elastic skeleton, the effective stress is given by the usual constitutive equation


where and are the Lamé’s first and second constant (shear modulus) of the skeleton under drained condition, respectively. is the strain.

Equation (1) is a system of coupled partial differential equations with , , and as the three unknown fields on the domain . The following are appropriate boundary conditions


where and are the prescribed skeleton displacement, traction, pressure and volumetric flux, and is the outward unit normal to the boundary . We note that for each of the flow and skeleton, the boundary is decomposed into two disjoint subsets which cover the entire boundary as shown in Figure (b)b.

3 Finite element formulation

The coupled nature of the poroelastic governing equations (1) precludes analytical solutions, except in some special cases like those presented in the seminal paper of Rice and Cleary [58]. In most cases, a numerical approximation must be computed. For this purpose, typically a finite element approach is used for spatial discretization, resulting in a system of differential algebraic equations (DAE) that can be solved by a time integration scheme. In this section, we discuss spatial discretization, specifically our approach of combining RT elements for the - fields with standard nodal Galerkin elements for the field. First we motivate our choice of a three-field formulation. Then we discuss the stability concerns that arise. Finally, we present the finite element formulation.

3.1 Three-field formulation

Zienkiewicz and Shiomi [3] first discussed various possible finite element formulations for Biot’s dynamic poroelastcity model with the notion of incompressibility. These can be classified into three-field (--) formulations and two-field (- or -) formulations. Zienkiewicz and Shiomi [3] did not present an implementation of a three-field formulation. A majority of current implementations are based on two-field formulations obtained by eliminating either the fluid Darcy velocity , resulting in a - formulation [59, 60, 61, 62, 20, 11, 5, 37] or the pressure, leading to a - formulation [63, 64, 3, 65, 66, 9]. These methods are efficient in terms of the computational cost as they reduce the number of DOFs. However, we pursue here a three-field formulation motivated by the following considerations.

Generality – When the fluid and the solid material making up the skeleton are modeled as incompressible, the pressure term does not appear in the mass conservation equation (1). Consequently, the fluid Darcy velocity, , cannot be directly eliminated from the Darcy’s law (1) [3, 22]. Therefore, to obtain a - formulation, either (a) penalty procedures must be used in the mass conservation equation [63, 64, 3, 65], or (b) the fluid relative acceleration, , (and hence inertia) must be neglected. The latter is a simplification, first introduced by Zienkiewicz and coworkers [59, 3], has been the most common means of obtaining a - formulation [59, 60, 61, 62, 20, 11, 5, 37]. This approach has been found to be typically valid for the slow and medium frequency level phenomena [59, 60]. However, in a variety of problems, it is desirable to take into account all acceleration terms; which is only possible by exploiting the three-field formulation, without using a penalty formulation. For instance using verification and validation procedures, Jeremić and coworkers [4, 52, 67] show that a three-field formulation is required to realistically model the saturated soil-foundation-structure interaction in earthquake as well as the liquefaction and cyclic mobility phenomena. They conclude that a full-field method can appropriately model the physical damping in the dynamic poroelasticity problem during seismic events by directly taking into account the fluid and skeleton interaction through Darcy’s equation. Accordingly, we could model and compute the velocity proportional damping energy as shown in equation (22) whereas in other studies an artificial damping has to be used [5]. Another area of interest is soft biological tissues that are often modeled mechanically as poroelastic media, and experience very rapid external forces [14]. This situation creates significant inertia forces in the fluid phase that is not negligible and requires use of a full formulation [12]. An improved performance of the full-field mixed finite element formulation over the two-field formulation is also reported in [10, 11] for large deformation simulations and in [68] for modeling of bones subjected to Osteoporosis disease. Presence of the fluid velocity field in an explicit fashion will also allow the formulation to be more readily extended for coupled transport phenomena.

Accuracy – The - formulation amounts to solving for pressure and skeleton displacement as primary variable fields. Typically post-processing techniques are required to compute the velocity field from the gradient of the pressure field, and is therefore of lower accuracy [21, 22]. Furthermore, both pressure and flux boundary conditions can be applied directly with a three-field formulation, but require special treatment if a two-field formulation is used [69].

3.2 Nature of stability considerations

To understand the nature of stability issues that are relevant in a three-field finite element formulation, we describe two limiting conditions in terms of hydraulic conductivity and rigidity of skeleton similar to [6, 42]. Both of these limiting cases give rise to saddle point problems. Without loss of generality, for these limiting cases we only consider the static version of the poroelasticity problem where all the inertia terms are zeroThe dynamic version affects only on ellipticity condition..

Case (i) – Low hydraulic conductivity: In the impermeable limit , the static version of (1) reduces to


where is the Laplacian operator. Note is replaced by as they are equivalent if the initial conditions satisfy. We recognize that these are the equations of the classical incompressible elasticity problem. Thus in the impermeable limit, there is no flow and the skeleton becomes incompressible but deformable.

Case (ii) – Rigid skeleton: The other limiting case is when the skeleton matrix become rigid (which implies ) and results in


These are the classical Darcy’s flow equations. We note that this second limiting case applies only to a three-field formulation.

Both these limit cases are saddle-point problems [70]. The coefficient matrix associated with the discrete form of the saddle point problem has a zero-block structure on lower diagonal (see equation (30) for example) that creates stability challenges. There are two stability criteria that provide the necessary and sufficient condition for well-posedness of such problems, namely the ellipticity condition and the LBB condition (also called inf-sup condition). We refer to [70, 71] for a thorough description on this subject. The discrete version of the ellipticity condition is satisfied for any choice of finite element with positive definite mass and stiffness matrices upon application of boundary conditions. The LBB stability condition is more important and challenging. Often, standard finite element approaches presented in the literature fail to satisfy this condition, leading to instability manifested as checkerboard patterns in the pressure field poor convergence behavior. We provide a comparison of stable and unstable finite element results in section 6.3.

Since the poroelasticity model degenerates to the saddle point problems (4) and (5), a finite element formulation must satisfy the stability conditions associated with both these limiting cases. Thus the poroelasticity problem has the additional requirement that both the - pair and the - pair must be able to satisfy the stability conditions simultaneously. We show in section 6.3.1 that the proposed element qualifies.

3.3 Finite element discretization

We start from the weak-form representation of (1)


where , , and are the skeleton displacement, fluid relative velocity, and pressure test functions respectively. We recognize now that the functions , , and must belong to the spaces


respectively. Here is the space of square integrable functions on , is the usual Sobolev space in which each function and its first order derivative belong to , and represents the space of vector-valued functions that belong to and whose divergences belong to as well. Further elaboration on these spaces can be found in [71] and also in [2, 55, 53]. Similarly, the test functions of , , and belong to the same spaces as their corresponding variable fields, except that instead of and a zero is placed. We call these test spaces , , and .

As discussed above, with the goal of satisfying the appropriate stability conditions, we simultaneously apply the lowest order Raviart-Thomas mixed finite element () [2] for flow variable fields and the standard linear/quadratic Galerkin finite elements for skeleton variable fields. Let the domain be partitioned with into non-overlapping triangular elements. We denote by and number of nodes and edges. Also position vector is .

Skeleton displacement discretization: We discretize the skeleton displacement field using standard linear/quadratic shape functions over triangles, ,


We note that the components of vector are skeleton displacement DOFs, , that are displacements at the element nodes. The skeleton and fluid DOFs are illustrated in Figure 4. We call the combination of linear Galerkin finite element with finite element, , and the combination of quadratic Galerkin finite element with finite element, . We propose the latter finite element to resolve the stability issues that the linear element creates in the limiting case of very small hydraulic conductivity as elaborated in section 6.3.

Fluid velocity discretization:


Based on element, Darcy velocity DOFs (Figure 4) are values of the normal velocity at the midpoints of edges of element (in fact, normal velocity is constant along each edge), defined as where is fluid velocity vector defined on edge . Also denotes the outward unit normal vector to edge of element , while is the unit normal vector to edge chosen with a global fixed orientation. This global convention is such that points from element with smaller index into its adjacent element with larger index (see Figure (a)a). Therefor, based on this convention, edges on the boundaries always have an outward global unit normal. In equation (10) is the position vector of the vertex node opposite the edge of element and is a sign convention equal to if points outward and otherwise is . The effect of this sign convention can be seen in Figure(b)b where the velocity interpolation function, , for the common edge () of the two adjacent elements is illustrated. Due to the convention, the element with smaller index has and the one with larger index has . Moreover, shown in Figure (c)c is the normal component of velocity along the common edge of two adjacent elements illustrating the continuity of such normal component. Also and refer to the area of element and the length of edge , respectively.

Pressure discretization: Pressure DOF is set to be constant elemental pressure, thus its interpolation function is a polynomial of degree 0 that takes a constant value of 1 at element and 0 at all other elements,


Fluid Darcy velocity and pressure DOFs are collected in vector and , respectively and are shown in Figure 4.

Figure 4: Fluid and skeleton DOFs in the coupled Galerkin/ elements – (a) linear Galerkin/ element, referred as , (b) quadratic Galerkin/ element referred as . Both elements are constructed on the natural topology of the variable fields, i.e., skeleton displacement fields are defined on the nodes using a standard Galerkin approximation, while fluid variable fields are defined on the element edges (velocity) and volumes (pressure) using element. We note that Darcy’s velocity DOFs are the values of the normal velocity at the midpoints of edges of elements and the pressure DOFs are set to be the constant elemental pressure.
Figure 8: (a) Two enumerated adjacent elements and that share the same edge . In this figure, and are the outward unit normal vector to edge of elements and , respectively, while is the unit normal vector to edge chosen with a global fixed orientation. This is such that is equal to the outward unit normal vector of element with smaller index, i.e., (and hence points into ). Also based on this convention, the element with smaller index has whereas the one with larger index has . (b) velocity interpolation function for the common edge of two adjacent elements, . (c) normal velocity along the common edge of two adjacent elements. Note the continuity of this component of velocity across element edge which results in element-wise mass conservation.

Given the described finite element approximation, we formulate a discrete version of the weak problem (6) by substituting approximation functions of in it.


The solution of the discretized problem must satisfies the integral equations for all test functions . In equation (11) denotes the Kronecker delta in the vector format, is the linearized deformation-displacement matrix, and is the elastic constitutive matrix under the drained condition. Note that we dropped the dependence of functions to and due to clarity. The test functions and the shape functions coincide in this formulations. Moreover, equation (11) represent a consistent variational formulation in which the same shape functions are applied for higher time derivatives. An alternative formulation is discussed in section 6.1.5. The prescribed skeleton displacement and fluid flux are imposed in a strong way, i.e., as essential boundary conditions, while the prescribed traction () and pressure () are entered the variational form as natural boundary conditions. Also, the constitutive equation (2) is used in strong form as substitute in (1). The compact form of (11) is the following semi-discrete expression


where the following terms are defined


We reiterate that such topology-based finite element discretization is physically compatible, and as a result all the above constant matrices (at the discrete level) imitates their corresponding operators at the continuous level in (1); for instance is the discrete version of divergence operator whose column represents the element where all entries are zeros except the rows representing the free edges belongs to that element. In those entries depending on the flux convention, is placed. Further elaboration on building these matrices is provided in Section 5. Also using the described finite element, specification of pressure and flux boundary conditions are straightforward.

4 Time Discretization

The resulting system of equations (12) not only consist of ordinary differential equations (ODE), but also algebraic equations. The numerical time integration of a system of Differential Algebraic Equations (DAE) is more complex than the solution of an ODE. Therefore, we utilize the Newmark constant average time integration scheme with a special attention to select a consistent initial condition (Section 4.1). Based on the Newmark method, governing equations of (12) is considered at


in which the following assumptions are made to approximate the values at each time step


where and

using approximations of (15) in the discretized equations of (14) results in the final time discretized version


where the following definitions are made


Rearranging and reverting Equation (16) to matrix form yields


represent a zero matrix of appropriate dimension. Above system of equations is symmetric with diagonal blocks and both symmetric positive definite.

4.1 Consistent initial condition

Following the three-field spatial approximation of the problem, the incompressibility condition, which is embedded in the mass conservation equation (12), forms an algebraic constraint. The presence of this constraint, makes (12) represent a DAE system in time rather than an ODE in time. This fact is also emphasized by Diebels and Ehlers [47]. We note the absence of the algebraic variable (pressure) from the constraint (12). We also note that not only the same matrix block appears in (18) and (18), but also the Schur complement of the upper left block matrix

is a nonsingular matrix with a bounded inverse. As a result, (18) yields to a Hessenberg index-2 system of DAEs as indicated by Ascher and Petzold [72] and has the form of an augmented Lagrangian equation. Consistency of initial conditions is one of the key aspects in DAEs§§§An alternative approach to deal with incompressible poroelasticity system of DAEs is to utilize the energy preserving method proposed by Bauchau [73]. His approach satisfies the algebraic constraint at the end point (if written in terms of displacement), while imposes the algebraic variable at the midpoint; hence, it does not require an initial condition on the algebraic variable. This method reduces to the constant average Newmark method over unconstrained problems and it preserve the energy of the system by eliminating of the work done by algebraic variable. The result of our Numerical Examples in Section 6 using the constant average Newmark method and Bauchau’s method perfectly match.; the initial values for variables in a DAE must satisfy not only the original equations in the system but also their differentials with respect to time [74, 75]. Accordingly, to select consistent initial conditions for acceleration terms and pressure, we must satisfy the original equations in the system as well as the differential of the mass conservation equation with respect to time, at time zero which reads


Solving this equation, is obtained. Note that in deriving 19, fluid and solid phases are assumed to be initially at rest without any deformation , which is an intuitive set of initial conditions that also satisfy the mass balance equation at time zero. We observe that an inconsistent set of initial conditions results in instability of pressure and acceleration time histories illustrated by nonphysical oscillations as numerical artifacts.

4.2 Equation of energy

In this section, we derive the work done by various forces acting on the system during a time step. We show that the work done by pressure (Lagrange multiplier) vanishes, as expected for incompressible media, so the total mechanical energy of the system is preserved.

The total incremental energy of the system is obtained by multiplying the discretized equations of motion (12) at the midpoint (i.e., average of the equations at time and ) by their respective displacement increments, and then adding them together. The resulting increment of the total energy is


where . Also , , , , , and denote the increments of kinetic energy of skeleton, kinetic energy of fluid, damping energy, strain energy of skeleton, work done by constraint force, and input energy respectively. Each of these terms are defined as


Where we employ the relationships of (15) for the displacement increments of solid and fluid. The statement of energy balance (20) says that the increment of the input energy must be equal to the increment of the kinetic energy, strain energy and damping forces energy. We note that the increment of the work done by constraint force is zero, , since the terms inside the parentheses in (21) are the discretized constraint equations that are set to be zero. Moreover, since we satisfy the constraint at initial time, we get as well.

Using (21), the energy components at each time step is recovered as


The damping and the input energy are time history dependent.

5 Implementation

The algorithm for implementation of the proposed numerical method for dynamic analysis of porous media is summarized in Procedure 1.

1:Form the coefficient matrices, , , , , , using Procedure (2)
2:Select a time increment (for an accurate result use Equation (27))
3:Calculate , , , using Equation (17)
4:Initialize , , to satisfy (12)
5:Initialize a consistent to satisfy (12) and (19) an example is given in Section 4.1
6:for  do
7:     Given , , and , calculate , , and using (13).
8:     Update and using Equation (17)
9:     Solve (18) for , , and
10:     Calculate (if required evaluate accelerations as well) using (15)
11:     For the energy balance analysis: Calculate , and using (22) Calculate and using (21) Calculate, and
12:end for
Procedure 1 Algorithm to solve dynamic poroelasticity

Procedure 2 gives an outline of the construction of coefficients defined in (13) for a two dimensional problem. Section 3.3 provides more details on the notations used here. and denote the number of free nodes and the number of free edges. In building these coefficient matrices, in addition to the connectivity (element to node) and nodal coordinates matrices used in the typical finite element problems, some additional data structures such as element to edge and node to edge are required. An example of the RT element implementation is elaborately discussed by Bahriawatiand and Carstensen [76]. To avoid clattering in Procedure 2, we introduce


These integrations over the element domains can be numerically evaluated by Gauss integration method.The following comments and Gauss quadrature rules over triangles [77] are considered for numerical examples of section 6

  • The order of the Gauss integration must be selected to accurately calculate the integrations. For instance, at least a six-point formula (three-point formula), which gives accurate results for up to fourth (second) degree of precision, is required to approximate a mass matrix with quadratic (linear) interpolation functions.

  • A reduction in the order of Gauss numerical integration from the order that is required to evaluate these matrices exactly, leads to inaccurate results with long lasting and high amplitude nonphysical oscillations.

  • To decrease the computational effort in case of a quadratic interpolation function, we use a sub-parametric element rather than an iso-parametric one, since all the elements in the numerical examples are constructed as straight edges, and the geometry can be interpreted to a lower degree than the nodal variable’s interpolation.

Moreover, in forming matrices of Procedure 2, we utilize the following properties of the mixed element


The right hand side equality is used in construction of matrix . The left hand side equality indicates that the velocity interpolation functions at edge , have the property of producing a unit flux through the edge and 0 flux through all other edges. Furthermore, according to our notation, in a typical triangle, the enumeration on the three vertex nodes is counterclockwise, while the edge enumeration is set such that the edge opposite to vertex node is the edge of the triangle.

1:Given definitions in (13)
2: ; ; ; ;
3:for  do
4:     Assemble , and from (23) into , and respectively
5:     Assemble into
6:     Assemble into
7:     Assemble into
8:end for
Procedure 2 Formation of coefficient matrices for two dimensional problems

6 Numerical Examples

To demonstrate the capabilities of our formulation to simulate the incompressible dynamical behavior of saturated porous media, we perform two benchmark numerical examples and verify them with analytical solutions (Example 1) and boundary element solutions (Example 2). As was noted in section 3.2, it has been found that when the hydraulic conductivity is very small, locking (in terms of failing the LBB condition) occurs [37, 56]. To explore this effect for our formulation, we consider a third set of examples to evaluate the stability of the method when locking is often a concern. Overall, we found that the presented numerical scheme is stable and sufficiently accurate for a wide range of material parameters including large and small values of hydraulic conductivity.

6.1 Example 1: half space – soil column

The first benchmark example is the application of a uniform load on the surface of a half space[65, 57, 5, 47, 3, 78]. Due to symmetry, only a representative column of small width, 0.1 m, is studied and the problem is of one dimension. The model is taken as an incompressible saturated soil with plane strain behavior. The material data are shown in Table I and are taken from the same benchmark example as in [47, 79]. The geometry and finite elements are illustrated in Figure 9, where the side walls and the bottom are impermeable and the displacements normal to their surfaces are constrained. The upper boundary is perfectly drained and subjected to a step load of 3 amplitude, which results in a total load of on top boundary. 400 finite elements with a crisscross mesh pattern are used to model the soil column. We also use a time step of (see section 6.1.2) that is adequately small for the purpose of wave propagation analysis. We have provided the results utilizing the element and a lumped mass matrix unless otherwise mentioned.

Example 1 14.516e+3 0.3 2000 1000 0.33 1e-2
Example 2 14.516e+3 0.3 2700 1000 0.42 1e-1 and 1e-4
Example 3 10 0.4 2667 1000 0.4 1e-7
Table I: Material properties for numercial examples
Figure 9: Example 1 —- a column of soil under step load: geometry and finite element mesh. We have isolated a 10 m length by 0.1 m width for analysis. The figure is not in scale. Three dots show the continuation.

With this example, we pursue many goals as listed below:

  1. Verification with analytical solutions

    • Short-term behavior: we compare the numerical results with the analytical solution of de Boer et al. [79] for wave propagation in an infinitely long column of incompressible saturated soil. The influence of the total depth of the column is negligible [47]; therefore, for our numerical analysis, we are able to isolate a 10-m length and still compare the short term behavior of our system with de Boer’s results for an infinite long column.

    • Long-term behavior: we compare long term behavior with Terzaghi’s equation for consolidation.

  2. Wave propagation analysis: according to the well known derivation of Biot [1], there are three types of waves that pass through a porous medium. Two dilatational waves (fast and slow) and one shear wave Accounting for apparent mass, there are two shear waves, one in fluid and one in solid skeleton[1].. It is proved analytically that as a consequence of the incompressibility condition, the fast dilatational wave propagates with an infinite speed and practically, only the slow dilatational wave exists [79, 80]. Also, note that for this one-dimensional problem, there is no shear waves.

  3. Pressure build-up and diffusion: we are interested in studying the initial and long term behavior of the pressure field.

  4. Effect of instantaneous loading: the impact and shock nature of the step load creates a challenging condition for the incompressible porous medium. The effect of this condition on the response history of the medium is studied, and some remedies are suggested to avoid the associated numerical artifacts.

  5. Effect of mass representation: it is suggested in the literature, that such numerical artifacts can be alleviated by appropriate representation of mass matrix. We explore different mass lumping methods and compare their results with those of the consistent mass matrix method.

  6. Energy balance analysis: in order to highlight the working of the presented numerical algorithm, we monitor the balance of the energy in the porous medium. In particular, by means of the energy balance equation, we can evaluate our time integration scheme and size of the time step.

6.1.1 Verification with analytical solutions

Shown in Figure 10 is the displacement response at the top boundary. Our numerical results match well with de Boer’s analytical solution in the transient region [79]. We point out that de Boer’s solution is for an infinitely long column where there is no reflection of wave from the bottom boundary; however, we isolate a finite length for computations, and therefore our numerical plots separates from de Boer’s solution after a while. Moreover, to verify the long term behavior, we use Terzaghi’s equation for steady state consolidation. According to that equation, the total settlement at the top boundary [81, equation 2-3] can be simplified for such one dimensional problems as:


where is the total load, is the length of the soil column, and is the the slope from a one dimensional test,


where a plain strain assumption is made.

Figure 10: Example 1 — displacement time history. Good agreement between numerical and the analytical computations is observed for both the transient behavior given by de Boer et al. [79] and the long term behavior given by Terzaghi 25

6.1.2 Wave propagation analysis

wave speed and choice of time step
The significance of the choice of an adequate time step is clear in a wave propagation problem. In order to represent the wave travel accurately, it is recommended that the time step is determined by the approximate Courant-Friedrichs-Lewy (CFL) condition (equation 9.102,[23])


This is the smallest time that is required for a dilatational wave to pass across any of the elements. Given the smallest element size, , we can compute knowing the wave speed . An estimate of the slow dilatational wave speed is given by de Boer et al.[79] as,


Using the parameters of Example 1 in the above relationship, we get . In the next section, we compute wave velocity using the plots of Figure 11 and demonstrate correlation with this estimate.

Setting the minimum mesh size, m, to be the short-edge length for triangles, equation (27) results in a time step of , so we apply a time step of .

Wave velocity analysis

Figure 11: Example 1— snapshots of wave propagation through the medium. The propagation speed reduces over the time as roughly calculate in (29).

We are interested in observing the passage of the slow dilatational wave which is clear in Figure 11 recorded at some snapshots of time. The time span was selected such that no reflections of the wave affect the results. In this plot, we see the propagation of the wave induced by the load at m and its propagation through the medium more and more as time progresses. Using this figure, the approximate wave velocity can be calculated as: . For instance, taking the first three plots we get:


Recall that an approximate value of 85 was calculated earlier using de Boer’s equation (28); that is close to the early time value of our numerical speed in Equation (29). However, our numerical results show a decreasing trend of wave speed as the wave travels through the media which is predicted by Schanz and Pryl [80] as well. One possible reason for this difference is the finite length of the column in our example compared to infinite length of de Boer analytical formulation.

Time histories of skeleton velocity at different locations , shown in Figure (b)b, also reveals the passage of the wave through the medium. We observe that moment when the wave reaches to a certain depth, the velocity profile of that depth separates from the rest of the profiles with zero velocity.

6.1.3 Pressure build-up and diffusion

Figure 14: Example 1— (a) Pressure time histories. Diffusion of pressure is apparent over the long time. (b) Zoom in of first moments. Pressure jumps upon loading at time zero. Also note the diffusion begins right after the passage of the stress wave.

Pore fluid pressure time histories computed at different depths is displayed in Figure (a)a for short and long term. Upon the instantaneous application of the load, a pressure build-up occurs at all points of the medium, indicating the propagation of the pressure at infinite speed [69](Figure (b)b). After the stress wave hits a depth, it induces a change in pressure which results in diffusion of the pressure back toward its initial state.

6.1.4 Effect of instantaneous loading

Figure 17: Example 1— (a) Solid skeleton velocity time histories for some representative depths. The trend, however, is the same at all locations with their velocities ultimately reach to zero. (b) Zoom in of first moments at some critical depths. The high frequency oscillations, byproducts of step loading, are more pronounced close to the top boundary, but they subside as time marches forward.

The time history of the skeleton velocity is shown in Figure 17. Over time, as energy dissipates and diffusion occurs, the velocity tends to zero (Figure (a)a). The trend in this plot is the same at all depths, so we only display the response at a few representative locations.

The step loading due to its abruptness creates a jump in pressure and inertia forces. This in turn excites some higher modes of vibration that survive only for a short time interval due to presence of physical damping in the form of fluid diffusion. These early high frequency oscillations are pronounced in the vicinity of the top boundary with rapid pressure change as shown in Figure (b)b. We recognize that these high frequency oscillations are a numerical artifact as their frequencies become unbounded upon mesh refinement. We note that decreasing the time step does not remove these spurious oscillations either. We also found that these nonphysical oscillations are not affected significantly (or removed) by finite element type, as presented in Figure 18. In these plots, we compare velocity response of the medium using and elements with the same number of elements. It is observed that the response using element has higher amplitude and frequency oscillations in comparison with element solution. However, as we move away form the loading boundary, the responses of both finite elements are indistinguishable.

Figure 18: Example 1 — a comparison of the solid skeleton velocity time histories. In the proximity of the loading boundary, the element exhibits higher amplitudes and frequency oscillations compared to the linear element. After a while, both elements converge to identical results.

To avoid these numerical oscillations, we explore different remedies – (1) Replacing the step load with a ramp load (even with a relatively short rise time) smooths out the initial oscillations (the results of this study are not shown here); (2) An alternative approach is to apply time integration methods containing algorithmic dissipation such as the Generalized method [82]. We note that the results shown here are computed by the constant average Newmark scheme where no algorithmic dissipation is used and consequently, the total energy of the system (20) is preserved. This fact is illustrated in Figure (b)b as well; (3) While such modifications of time integration algorithm can be easily implemented, we instead search for remedies at the element formulation level. In [83], it is suggested that such oscillations may be associated with representation of mass matrix. This is discussed in the next section 6.1.5. We also point out that these early numerical oscillations do not appear in displacement time histories since the integration process automatically filters out the high frequency components.

6.1.5 Effects of mass representation

In constructing the matrix defined in (13), we use the same shape functions for velocities and displacements, resulting in a consistent variational formulation as presented in (11). Accordingly, represents a consistent mass matrix, which is generally non-diagonal. An alternative to the consistent mass matrix is to employ a diagonal mass matrix based on direct lumping. In this method, the total mass of element is directly apportioned to nodal DOFs, ignoring any cross coupling. There are several ways of constructing a diagonally lumped mass matrix. Here, we test two of them, namely, mass lumping by nodal Lobatto quadrature [84] and the lumping method proposed by Hinton [85]. We note that in this example the lumping methods are applied on the mass matrices associated with the block related to the skeleton displacement, i.e, and The lumping process can be viewed as substituting a different interpolation function for velocity (as opposed to using the same interpolation function as used for displacement) to approximate the mass matrix as, . Accordingly, to build the coupled mass matrix we must use the same interpolation function in order to tune it with the lumped mass ., while we keep using the same mass matrix for the fluid velocity block, i.e., ; since there is no notion of fluid displacement in our formulation and instead the primary variable filed is fluid velocity. The resulting weighting factor of each DOF in the element mass matrix is illustrated in Figure (a)a. Note that over the linear triangle which is used in the element, both procedures give rise to the identical nodal weights.

Figure 22: Example 1 — (a) Weighting factor of each nodal DOF in the element lumped mass matrix. Note that over the linear triangle element, , both lumping procedures result in identical nodal weights. (b) Comparison of lumped and consistent mass matrix responses for the element. The consistent mass matrix invokes a greater amount of high frequency oscillations, especially upon the arrival of the wave front. (c) Unstable response of the quadratic triangle element, , when a Lobatto quadrature lumping method is used.

Shown in Figure (b)b, we provide the numerical results employing both the lumped and consistent mass matrices for the linear triangle element, . The lumped mass matrix results are displayed in solid line while the consistent mass matrix results are in dashed line. We see that the consistent mass matrix creates oscillations with higher frequencies compared to the lumped mass matrix response. In other words, it is more noisy compared to lumped mass response especially upon the arrival of the wave front. Nevertheless, in both cases, noise is damped out after a short while.

With the quadratic triangle element, , we see the same trend as for linear triangle element when the Hinton et al. method of lumping is applied. However, the use of the Lobatto quadrature lumping method gives rise to an unbounded and unstable solution as plotted in Figure (c)c. This numerical instability occurs regardless of the mesh size, the time step size, and the hydraulic conductivity. In fact, we observe that the Lobatto quadrature lumping method over the quadratic triangular element produces unstable numerical solution even in a typical dynamic elasticity problem. We also emphasize that when using a quadratic triangle element, a consistent load vector as defined in (13) should be exerted regardless of mass matrix begin lumped or consistent. Otherwise, very inaccurate results are obtained as reported in [23] as well.

6.1.6 Energy plots

We have developed the energy balance equation for the entire porous medium in section 4.2. As a means of highlighting the performance and accuracy of the presented time integration approach, the relative error in the energy balance of the solution is computed as shown in Figure (a)a. Left hand side (LHS) energy is defined as the summation of all the energy terms except the input energy as expressed in (22). The relative error is defined as the absolute value of the difference of the LHS energy and input energy divided by the input energy. The error is very small and the time step scheme is accurate. Also all the components of energy are illustrated in Figure (b)b.

Figure 25: Example 1 — (a) Energy balance error. (b) Components of energy as defined in (22)

6.1.7 Summary of findings

  • As shown in Figure 10, a good agreement between the numerical and the analytical results is observed for both the transient behavior given by de Boer et al. [79] for wave propagation analysis, and the long term behavior given by Terzaghi’s consolidation equation.

  • We observe a decreasing trend of wave speed as it travels through the medium.

  • Pressure diffusion begins right after the passage of the stress wave as shown in Figure (a)a.

  • Due to the instantaneous nature of the step loading, high frequency numerical modes are excited early in time. These are shown in the velocity time history of Figure (b)b. We suggest using a small ramp or a time integration method with algorithmic dissipation to alleviate this problem. Also a mass lumping method is shown to be effective in reduction of these types of numerical artifacts. We note that these high frequencies are not apparent in displacement time history.

  • The use of the Lobatto quadrature mass lumping method gives rise to an unbounded and unstable solution as shown in Figure (c)c.

  • The small relative error in the energy balance of numerical solution expresses the satisfactory performance and accuracy of the presented time integration method.

6.2 Example 2: soil block with a partially loaded surface

A key benchmark for verification of numerical scheme is the analysis of a saturated soil block with a partially loaded surface as shown in Figure 26. Due to symmetry only half of the problem is modeled. All the boundaries except the unloaded area of the upper boundary are impermeable. The displacements normal to the surfaces of the side walls and the bottom are constrained. We adopt the material properties that are used by Li et al. [