A new mixed finite element method for poro‐elasticity

Development of robust numerical solutions for poro‐elasticity is an important and timely issue in modern computational geomechanics. Recently, research in this area has seen a surge in activity, not only because of increased interest in coupled problems relevant to the petroleum industry, but also due to emerging applications of poro‐elasticity for modelling problems in biomedical engineering and materials science. In this paper, an original mixed least‐squares method for solving Biot consolidation problems is developed. The solution is obtained via minimization of a least‐squares functional, based upon the equations of equilibrium, the equations of continuity and weak forms of the constitutive relationships for elasticity and Darcy flow. The formulation involves four separate categories of unknowns: displacements, stresses, fluid pressures and velocities. Each of these unknowns is approximated by linear continuous functions. The mathematical formulation is implemented in an original computer program, written from scratch and using object‐oriented logic. The performance of the method is tested on one‐ and two‐dimensional classical problems in poro‐elasticity. The numerical experiments suggest the same rates of convergence for all four types of variables, when the same interpolation spaces are used. The continuous linear triangles show the same rates of convergence for both compressible and entirely incompressible elastic solids. This mixed formulation results in non‐oscillating fluid pressures over entire domain for different moments of time. The method appears to be naturally stable, without any need of additional stabilization terms with mesh‐dependent parameters. Copyright © 2007 John Wiley & Sons, Ltd.     

A FIC‐based stabilized mixed finite element method with equal order interpolation for solid–pore fluid interaction problems

A new mixed displacement‐pressure element for solving solid–pore fluid interaction problems is presented. In the resulting coupled system of equations, the balance of momentum equation remains unaltered, while the mass balance equation for the pore fluid is stabilized with the inclusion of higher‐order terms multiplied by arbitrary dimensions in space, following the finite calculus (FIC) procedure. The stabilized FIC‐FEM formulation can be applied to any kind of interpolation for the displacements and the pressure, but in this work, we have used linear elements of equal order interpolation for both set of unknowns. Examples in 2D and 3D are presented to illustrate the accuracy of the stabilized formulation for solid–pore fluid interaction problems. Copyright © 2016 John Wiley & Sons, Ltd.     

Analytical solution for one‐dimensional consolidation of unsaturated soils using eigenfunction expansion method

This paper introduces an exact analytical solution for governing flow equations for one‐dimensional consolidation in unsaturated soil stratum using the techniques of eigenfunction expansion and Laplace transformation. The homogeneous boundary conditions adopted in this study are as follows: (i) a one‐way drainage system of homogenous soils, in which the top surface is considered as permeable to air and water, whereas the base is an impervious bedrock; and (ii) a two‐way drainage system where both soil ends allow free dissipation of pore‐air and pore‐water pressures. In addition, the analytical development adopts initial conditions capturing both uniform and linear distributions of the initial excess pore pressures within the soil stratum. Eigenfunctions and eigenvalues are parts of the general solution and can be obtained based on the proposed boundary conditions. Besides, the Laplace transform method is adopted to solve the first‐order differential equations. Once equations with transformed domain are all obtained, the final solutions, which are proposed to be functions of time and depth, can be achieved by taking an inverse Laplace transform. To verify the proposed solution, two worked examples are provided to present the consolidation characteristics of unsaturated soils based on the proposed method. The validation of the recent results against other existing analytical solutions is graphically demonstrated. Copyright © 2013 John Wiley & Sons, Ltd.     

Coupled formulation and algorithms for the simulation of non‐planar three‐dimensional hydraulic fractures using the generalized finite element method

This paper presents a coupled hydro‐mechanical formulation for the simulation of non‐planar three‐dimensional hydraulic fractures. Deformation in the rock is modeled using linear elasticity, and the lubrication theory is adopted for the fluid flow in the fracture. The governing equations of the fluid flow and elasticity and the subsequent discretization are fully coupled. A Generalized/eXtended Finite Element Method (G/XFEM) is adopted for the discretization of the coupled system of equations. A Newton–Raphson method is used to solve the resulting system of nonlinear equations. A discretization strategy for the fluid flow problem on non‐planar three‐dimensional surfaces and a computationally efficient strategy for handling time integration combined with mesh adaptivity are also presented. Several three‐dimensional numerical verification examples are solved. The examples illustrate the generality and accuracy of the proposed coupled formulation and discretization strategies. Copyright © 2015 John Wiley & Sons, Ltd.     

Iteratively coupled solution strategies for a four‐field mixed finite element method for poroelasticity

In this paper, we consider numerical algorithms for modeling of the time‐dependent coupling between the fluid flow and deformation in elastic porous media. Here, we employ a four‐field formulation which uses the total stress, displacement, flux, and pressure as its primary variables and satisfies Darcy's law and linear elasticity in mixed weak form. We present four different iteratively coupled methods, known as drained, undrained, fixed‐strain, and fixed‐stress splits, in which the diffusion operator is separated from the elasticity operator and the two subproblems are solved in a staggered way while ensuring convergence of the solution at each time step. A‐priori convergence results for each iterative coupling which differs from those found when using a traditional two‐field or three‐field formulation are presented. We also present some numerical results to support the convergence estimates and to show the accuracy and efficiency of the algorithms. Copyright © 2016 John Wiley & Sons, Ltd.     

Three‐dimensional transient thermo‐hydro‐mechanical fundamental solutions of unsaturated soils

The governing differential equations of unsaturated soils considering the thermo‐poro‐mechanical behaviour consist of equilibrium, moisture air and heat transfer equations. In this paper at first, following some necessary simplifications, the thermal three‐dimensional fundamental solution for an unsaturated deformable porous medium with linear elastic behaviour in Laplace transform domain is presented. Subsequently, the closed‐form time domain fundamental solutions are derived by analytical inversion of the Laplace transform domain solutions. Then a set of numerical results are presented, which demonstrate the accuracies and some salient features of the derived analytical transient fundamental solutions. Finally, the closed‐form time domain fundamental solution will be verified mathematically by comparison with the previously introduced corresponding fundamental solution. Copyright © 2009 John Wiley & Sons, Ltd.     

On finite volume method implementation of poro‐elasto‐plasticity soil model

Accurate prediction of the interactions between the nonlinear soil skeleton and the pore fluid under loading plays a vital role in many geotechnical applications. It is therefore important to develop a numerical method that can effectively capture this nonlinear soil‐pore fluid coupling effect. This paper presents the implementation of a new finite volume method code of poro‐elasto‐plasticity soil model. The model is formulated on the basis of Biot's consolidation theory and combined with a perfect plasticity Mohr‐Coulomb constitutive relation. The governing equation system is discretized in a segregated manner, namely, those conventional linear and uncoupled terms are treated implicitly, while those nonlinear and coupled terms are treated explicitly by using any available values from previous time or iteration step. The implicit–explicit discretization leads to a linearized and decoupled algebraic system, which is solved using the fixed‐point iteration method. Upon the convergence of the iterative method, fully nonlinear coupled solutions are obtained. Also explored in this paper is the special way of treating traction boundary in finite volume method compared with FEM. Finally, three numerical test cases are simulated to verify the implementation procedure. It is shown in the simulation results that the implemented solver is capable of and efficient at predicting reasonable soil responses with pore pressure coupling under different loading situations. Copyright © 2015 John Wiley & Sons, Ltd.     

Analytical study of mine closure behaviour in a poro‐visco‐elastic medium

This paper is interested in the hydro‐mechanical behaviour of an underground cavity abandoned at the end of its service life. It is an extension of a previous study that accounted for a poro‐elastic behaviour of the rock mass (Int. J. Comput. Geomech. 2007; DOI: 10.1016/j.compgeo.2007.11.003 ). Deterioration of the lining support with time leads to the transfer of the loading from the exterior massif to the interior backfill. The in situ material has a poro‐visco‐elastic constitutive behaviour while the backfill is poro‐elastic, both saturated with water. This loading transfer is accompanied by an inward cavity convergence, thereby compressing the backfill, and induces an outward water flow. This leads to a complex space–time evolution of pore pressures, displacements and stresses, which is not always intuitive. In its general setting, a semi‐explicit solution to this problem is developed, using Laplace transform, the inversion being performed numerically. Analytical inversion leading to a quasi‐explicit solution in the time domain is possible by identifying the characteristic creep and relaxation times of volumetric strains with those of the deviatoric strains, on the basis of a parametric study. A few numerical examples are given to illustrate the hydro‐mechanical behaviour of the cavity and highlight the influence of key parameters (e.g. stiffness of backfill, lining deterioration rate, etc.). Further studies accounting for more general material behaviours for the backfill and external ground are ongoing. Copyright © 2008 John Wiley & Sons, Ltd.     

A closed form analytical solution for two‐dimensional plane strain consolidation of unsaturated soil stratum

This paper discusses the excess pore‐air and pore‐water pressure dissipations and the average degree of consolidation in the 2D plane strain consolidation of an unsaturated soil stratum using eigenfunction expansion and Laplace transformation techniques. In this study, the application of a constant external loading on a soil surface is assumed to immediately generate uniformly or linearly distributed initial excess pore pressures. The general solutions consisting of eigenfunctions and eigenvalues are first proposed. The Laplace transform is then applied to convert the time variable t in partial differential equations into the Laplace complex argument s. Once the domain is obtained, a simplified set of equations with variable s can be achieved. The final analytical solutions can be computed by taking a Laplace inverse. The proposed equations predict the two‐dimensional consolidation behaviour of an unsaturated soil stratum capturing the uniformly and linearly distributed initial excess pore pressures. This study investigates the effects of isotropic and anisotropic permeability conditions on variations of excess pore pressures and the average degree of consolidation. Additionally, isochrones of excess pore pressures along vertical and horizontal directions are presented. It is found that the initial distribution of pore pressures, varying with depth, results in considerable effects on the pore‐water pressure dissipation rate whilst it has insignificant effects on the excess pore‐air pressure dissipation rate. Copyright © 2015 John Wiley & Sons, Ltd.     

Analysis of undrained cavity expansion in elasto‐plastic soils with non‐linear elasticity

A large strain analysis of undrained expansion of a spherical/cylindrical cavity in a soil modelled as non‐linear elastic modified Cam clay material is presented. The stress–strain response of the soil is assumed to obey non‐linear elasticity until yielding. A power‐law characteristic or a hyperbolic stress–strain curve is used to describe the gradual reduction of soil stiffness with shear strain. It is assumed that, after yielding, the elasto‐plastic behaviour of the soil can be described by the modified Cam clay model. Based on a closed‐form stress–strain response in undrained condition, a numerical solution is obtained with the aid of simple numerical integration technique. The results show that the stresses and the pore pressure in the soil around an expanded cavity are significantly affected by the non‐linear elasticity, especially if the soil is overconsolidated. The difference between large strain and small strain solutions in the elastic zone is not significant. The stresses and the pore pressure at the cavity wall can be expressed as an approximate closed‐form solution. Copyright © 2001 John Wiley & Sons, Ltd.     

Semi‐analytical solution to one‐dimensional consolidation for unsaturated soils with semi‐permeable drainage boundary under time‐dependent loading

This paper presents semi‐analytical solutions to Fredlund and Hasan's one‐dimensional consolidation of unsaturated soils with semi‐permeable drainage boundary under time‐dependent loadings. Two variables are introduced to transform two coupled governing equations of pore‐water and pore‐air pressures into an equivalent set of partial differential equations, which are easily solved by the Laplace transform. The pore‐water pressure, pore‐air pressure and settlement are obtained in the Laplace domain. Crump's method is adopted to perform the inverse Laplace transform in order to obtain semi‐analytical solutions in time domain. It is shown that the present solutions are more general and have a good agreement with the existing solutions from literatures. Furthermore, the current solutions can also be degenerated into conventional solutions to one‐dimensional consolidation of unsaturated soils with homogeneous boundaries. Finally, several numerical examples are provided to illustrate consolidation behavior of unsaturated soils under four types of time‐dependent loadings, including instantaneous loading, ramp loading, exponential loading and sinusoidal loading. Parametric studies are illustrated by variations of pore‐air pressure, pore‐water pressure and settlement at different values of the ratio of air–water permeability coefficient, depth and loading parameters. Copyright © 2017 John Wiley & Sons, Ltd.     

Some numerical issues using element‐free Galerkin mesh‐less method for coupled hydro‐mechanical problems

A new formulation of the element‐free Galerkin (EFG) method is developed for solving coupled hydro‐mechanical problems. The numerical approach is based on solving the two governing partial differential equations of equilibrium and continuity of pore water simultaneously. Spatial variables in the weak form, i.e. displacement increment and pore water pressure increment, are discretized using the same EFG shape functions. An incremental constrained Galerkin weak form is used to create the discrete system equations and a fully implicit scheme is used for discretization in the time domain. Implementation of essential boundary conditions is based on a penalty method. Numerical stability of the developed formulation is examined in order to achieve appropriate accuracy of the EFG solution for coupled hydro‐mechanical problems. Examples are studied and compared with closed‐form or finite element method solutions to demonstrate the validity of the developed model and its capabilities. The results indicate that the EFG method is capable of handling coupled problems in saturated porous media and can predict well both the soil deformation and variation of pore water pressure over time. Some guidelines are proposed to guarantee the accuracy of the EFG solution for coupled hydro‐mechanical problems. Copyright © 2008 John Wiley & Sons, Ltd.     

Pore‐scale modeling of fluid‐particles interaction and emerging poromechanical effects

A micro‐hydromechanical model for granular materials is presented. It combines the discrete element method for the modeling of the solid phase and a pore‐scale finite volume formulation for the flow of an incompressible pore fluid. The coupling equations are derived and contrasted against the equations of conventional poroelasticity. An analogy is found between the discrete element method pore‐scale finite volume coupling and Biot's theory in the limit case of incompressible phases. The simulation of an oedometer test validates the coupling scheme and demonstrates the ability of the model to capture strong poromechanical effects. A detailed analysis of microscale strain and stress confirms the analogy with poroelasticity. An immersed deposition problem is finally simulated and shows the potential of the method to handle phase transitions. Copyright © 2013 John Wiley & Sons, Ltd.     

Poro‐elasto‐plastic response of an unconsolidated formation confined with stiff seal rocks under radial injection

In this article, we evaluate geomechanics of fluid injection from a fully penetrating vertical well into an unconsolidated formation confined with stiff seal rocks. The coupled behavior of an isotropic, homogeneous sand layer is studied under injection pressures that are high enough to induce plasticity yet not fracturing. Propagation of the significant influence zone surrounding the injection borehole, quantified by the extent of the plastic domain in the elasto‐plastic model, is examined for the first time. First, a new fully coupled axisymmetric numerical model is developed. A comprehensive assessment is performed on pore pressures, stresses/strains, and failure planes during the entire transient period of an injection cycle. Results anticipate existence of five distinctive zones in terms of plasticity state: liquefaction at wellbore; two inner plastic domains surrounding the wellbore, where failure occurs along two planes and major principal stress is in vertical direction; remaining of the plastic domain, where formation fails along one plane and major principal stress is in radial direction; and a non‐plastic region. Failure mechanism at the wellbore is found to be shear followed by liquefaction. Next, a novel methodology is proposed based on which new weakly coupled poro‐elasto‐plastic analytical solutions are derived for all three stress/strain components. Unlike previous studies, extension of the plastic zone is obtained as a function of injection pressure, incorporating plasticity effects on the subsequent elastic domain. Solutions, proven to be a good approximation of numerical simulations, offer a huge advantage as the run time of coupled numerical simulations is considerably long. Copyright © 2016 John Wiley & Sons, Ltd.     

Fast iterative solution of large undrained soil‐structure interaction problems

In view of rapid developments in iterative solvers, it is timely to re‐examine the merits of using mixed formulation for incompressible problems. This paper presents extensive numerical studies to compare the accuracy of undrained solutions resulting from the standard displacement formulation with a penalty term and the two‐field mixed formulation. The standard displacement and two‐field mixed formulations are solved using both direct and iterative approaches to assess if it is cost‐effective to achieve more accurate solutions. Numerical studies of a simple footing problem show that the mixed formulation is able to solve the incompressible problem ‘exactly’, does not create pressure and stress instabilities, and obviate the need for an ad hoc penalty number. In addition, for large‐scale problems where it is not possible to perform direct solutions entirely within available random access memory, it turns out that the larger system of equations from mixed formulation also can be solved much more efficiently than the smaller system of equations arising from standard formulation by using the symmetric quasi‐minimal residual (SQMR) method with the generalized Jacobi (GJ) preconditioner. Iterative solution by SQMR with GJ preconditioning also is more elegant, faster, and more accurate than the popular Uzawa method. Copyright © 2003 John Wiley & Sons, Ltd.     

Fundamental solutions for a fluid‐saturated,transversely isotropic,poroelastic solid

The fundamental solutions were obtained for step‐like point forces acting in three orthogonal directions and an instantaneous fluid point source in a fluid‐saturated, porous, infinite solid of transversely isotropic elasticity and permeability. After expressing the governing equations in the form of matrix in the Laplace space, we employed Kupradze's method together with the triple Fourier transforms. This method reduces the simultaneous partial differential equations with respect to three displacement components and a pore fluid pressure to a differential equation in terms of only one potential scalar function, which can be operationally solved in the transformed space. After the Laplace inversion of the potential, the residue theorem was applied to its Fourier inverse transform with respect to one of the transformation variables. The Fourier transforms with respect to two other variables were rewritten into the Hankel transforms. Copyright © 2002 John Wiley & Sons, Ltd.     

Two‐dimensional finite element analysis of gravitational and lateral‐driven deformation in sedimentary basins

The paper deals with the modeling of some aspects, such as the formulation of constitutive equations for sediment material or finite element approach for basin analysis, related to mechanical compaction in sedimentary basins. In addition to compaction due to gravity forces and pore‐pressure dissipation, particular emphasis is given to the study of deformation induced by tectonic sequences. The numerical model relies upon the implementation of a comprehensive constitutive model for the sediment material formulated within the framework of finite poroplasticity. The theoretical model accounts for both hydromechanical and elasticity–plasticity coupling due to the effects of irreversible large strains. From the numerical viewpoint, a finite element procedure specifically devised for dealing with sedimentary basins as open systems allows to simulate within a two‐dimensional setting the process of sediment accretion or erosion. Several basin simulations are presented. The main objective is to analyze the behavior of a sedimentary basin during the different phases of its life cycle: accretion phase, pore‐pressure dissipation phase and compressive/extensional tectonic motions. Copyright © 2013 John Wiley & Sons, Ltd.     

A continuum‐discrete model using Darcy's law: formulation and verification

This paper presents a numerical scheme for fluid‐particle coupled discrete element method (DEM), which is based on poro‐elasticity. The motion of the particles is resolved by means of DEM. While within the proposition of Darcian regime, the fluid is assumed as a continuum phase on a Eulerian mesh, and the continuity equation on the fluid mesh for a compressible fluid is solved using the FEM. Analytical solutions of traditional soil mechanics examples, such as the isotropic compression and one‐dimensional upward seepage flow, were used to validate the proposed algorithm quantitatively. The numerical results showed very good agreement with the analytical solutions, which show the correctness of this algorithm. Sensitivity studies on the effect of some influential factors of the coupling scheme such as pore fluid bulk modulus, volumetric strain calculation, and fluid mesh size were performed to display the accuracy, efficiency, and robustness of the numerical algorithm. It is revealed that the pore fluid bulk modulus is a critical parameter that can affect the accuracy of the results. Because of the iterative coupling scheme of these algorithms, high value of fluid bulk modulus can result in instability and consequently reduction in the maximum possible time‐step. Furthermore, the increase of the fluid mesh size reduces the accuracy of the calculated pore pressure. This study enhances our current understanding of the capacity of fluid‐particle coupled DEM to simulate the mechanical behavior of saturated granular materials. Copyright © 2014 John Wiley & Sons, Ltd.     

Dynamics of unsaturated soils using various finite element formulations

Unsaturated soils are three‐phase porous media consisting of a solid skeleton, pore liquid, and pore gas. The coupled mathematical equations representing the dynamics of unsaturated soils can be derived based on the theory of mixtures. Solution of these fully coupled governing equations for unsaturated soils requires tremendous computational resources because three individual phases and interactions between them have to be taken into account. The fully coupled equations governing the dynamics of unsaturated soils are first presented and then two finite element formulations of the governing equations are presented and implemented within a finite element framework. The finite element implementation of all the terms in the governing equations results in the complete formulation and is solved for the first time in this paper. A computationally efficient reduced formulation is obtained by neglecting the relative accelerations and velocities of liquid and gas in the governing equations to investigate the effects of fluid flow in the overall behavior. These two formulations are used to simulate the behavior of an unsaturated silty soil embankment subjected to base shaking and compared with the results from another commonly used partially reduced formulation that neglects the relative accelerations, but takes into account the relative velocities. The stress–strain response of the solid skeleton is modeled as both elastic and elastoplastic in all three analyses. In the elastic analyses no permanent deformations are predicted and the displacements of the partially reduced formulation are in between those of the reduced and complete formulations. The frequency of vibration of the complete formulation in the elastic analysis is closer to the predominant frequency of the base motion and smaller than the frequencies of vibration of the other two analyses. Proper consideration of damping due to fluid flows in the complete formulation is the likely reason for this difference. Permanent deformations are predicted by all three formulations for the elastoplastic analyses. The complete formulation, however, predicts reductions in pore fluid pressures following strong shaking resulting in somewhat smaller displacements than the reduced formulation. The results from complete and reduced formulations are otherwise comparable for elastoplastic analyses. For the elastoplastic analysis, the partially reduced formulation leads to stiffer response than the other two formulations. The likely reason for this stiffer response in the elastoplastic analysis is the interpolation scheme (linear displacement and linear pore fluid pressures) used in the finite element implementation of the partially reduced formulation. Copyright © 2008 John Wiley & Sons, Ltd.     

Effect of wave non‐linearity on the standing‐wave‐induced seabed response

A standing wave in front of a seawall may reach a height more than twice of its incident component. When excess pore pressure occurs, it may even induce seabed instability, hence endangering the structure. This issue was studied previously using only linear wave theory. In this paper, standing‐wave theory to a second‐order approximation is applied, in order to demonstrate the differences between these two solutions. The spatial and temporal variations in the instantaneous pore pressure are first calculated, in addition to their vertical distributions. The effects of wave height, water depth and the degree of soil saturation on pore pressure distributions are then discussed, followed by the net pore pressure averaged over one wave cycle. The results suggest the existence of a residual pore pressure in the seabed and its net pore pressure can be used to estimate the wave‐induced liquefaction potential in a soil column. It also indicates that, in deep water, the second‐order solution predicts that a negative pore pressure at an antinode which may be greater than a positive pressure. Overall, the second‐order solution is found to agree better with the experimental results of the pore pressures available, compared to the linear solution. Copyright © 2000 John Wiley & Sons, Ltd.