One‐dimensional dynamics of saturated incompressible porous media: analytical solutions and influence of inertia terms

The complexity of formulations for the hydromechanical coupled mechanics of porous media is typically minimised by simplifying assumptions such as neglecting the effect of inertia terms. For example, three formulations commonly employed to model practical problems are classified as fully dynamic, simplified dynamic and quasi‐static. Thus, depending on the porous media conditions, each formulation will have advantages and limitations. This paper presents a comprehensive analysis of these limitations when solving one‐dimensional fully saturated porous media problems in addition to a new solution that considers a more general loading situation. A phase diagram is developed to assist on the selection of which formulation is more appropriate and convenient regarding particular cases of porosity and hydraulic conductivity values. Non‐dimensional formulations are proposed to achieve this goal. Results using the analytical solutions are compared against numerical values obtained with the finite element method, and the effect of porosity is investigated. Copyright © 2016 John Wiley & Sons, Ltd.     

Effect of local nonlinearity in cohesionless soil on optimal radius minimizing fixed-head pile bending by inertial and kinematic interactions

This study presents the effects of a local nonlinearity in cohesionless soil upon the optimal radius minimizing the bending strains of a vertical, cylindrical fixed-head pile embedded in a layered soil stratum in a soil–pile–structure system where the kinematic interaction dominates. The seismic deformation method (SDM) with discretized numerical models is applied since the SDM is a static numerical method that can easily consider realistic conditions of layered soil strata and the nonlinearity of the soil. In the numerical models, the local nonlinearity of the soil in the vicinity of the pile is represented by subgrade springs having bi-linear skeleton curves with a simple hysteretic loop. Various amplitudes of the lateral displacements of the soil and the lateral forces at the head of the pile are considered as numerical parameters. The results of parametric analyses reveal the presence of an optimal pile radius that locally minimizes the bending strains of the piles under strong nonlinearity of the soil, and the optimal pile radius tends to increase as the degree of nonlinearity increases. Criteria are presented for predicting the increment in the optimal radius of soil–pile–structure systems under strong nonlinearity in the soil.     

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.     

Higher resolution total velocity <Emphasis Type="Italic">Vt</Emphasis> and <Emphasis Type="Italic">Va</Emphasis> finite-volume formulations on cell-centred structured and unstructured grids

Novel cell-centred finite-volume formulations are presented for incompressible and immiscible two-phase flow with both gravity and capillary pressure effects on structured and unstructured grids. The Darcy-flux is approximated by a control-volume distributed multipoint flux approximation (CVD-MPFA) coupled with a higher resolution approximation for convective transport. The CVD-MPFA method is used for Darcy-flux approximation involving pressure, gravity, and capillary pressure flux operators. Two IMPES formulations for coupling the pressure equation with fluid transport are presented. The first is based on the classical total velocity Vt fractional flow (Buckley Leverett) formulation, and the second is based on a more recent Va formulation. The CVD-MPFA method is employed for both Vt and Va formulations. The advantages of both coupled formulations are contrasted. The methods are tested on a range of structured and unstructured quadrilateral and triangular grids. The tests show that the resulting methods are found to be comparable for a number of classical cases, including channel flow problems. However, when gravity is present, flow regimes are identified where the Va formulation becomes locally unstable, in contrast to the total velocity formulation. The test cases also show the advantages of the higher resolution method compared to standard first-order single-point upstream weighting.     

Application of MODFLOW and geographic information system to groundwater flow simulation in North China Plain, China

MODFLOW is a groundwater modeling program. It can be compiled and remedied according to the practical applications. Because of its structure and fixed data format, MODFLOW can be integrated with Geographic Information Systems (GIS) technology for water resource management. The North China Plain (NCP), which is the politic, economic and cultural center of China, is facing with water resources shortage and water pollution. Groundwater is the main water resource for industrial, agricultural and domestic usage. It is necessary to evaluate the groundwater resources of the NCP as an entire aquifer system. With the development of computer and internet information technology it is also necessary to integrate the groundwater model with the GIS technology. Because the geological and hydrogeological data in the NCP was mainly in MAPGIS format, the powerful function of GIS of disposing of and analyzing spatial data and computer languages such as Visual C and Visual Basic were used to define the relationship between the original data and model data. After analyzing the geological and hydrogeological conditions of the NCP, the groundwater flow numerical simulation modeling was constructed with MODFLOW. On the basis of GIS, a dynamic evaluation system for groundwater resources under the internet circumstance was completed. During the process of constructing the groundwater model, a water budget was analyzed, which showed a negative budget in the NCP. The simulation period was from 1 January 2002 to 31 December 2003. During this period, the total recharge of the groundwater system was 49,374 × 10⁶ m³ and the total discharge was 56,530 × 10⁶ m³ the budget deficit was −7,156 × 10⁶ m³. In this integrated system, the original data including graphs and attribution data could be stored in the database. When the process of evaluating and predicting groundwater flow was started, these data were transformed into files that the core program of MODFLOW could read. The calculated water level and drawdown could be displayed and reviewed online.     

A new locally conservative numerical method for two-phase flow in heterogeneous poroelastic media

We construct a new class of locally conservative numerical methods for two-phase immiscible flow in heterogeneous poroelastic media. Within the framework of the so-called iteratively coupled methods and fixed-stress split algorithm we develop mixed finite element methods for the flow and geomechanics subsystems which furnish locally conservative Darcy velocity and transient porosity input fields for the transport problem for the water saturation. Such hyperbolic equation is decomposed within an operator splitting technique based on a predictor–corrector scheme with the predictor step discretized by a higher-order non-oscillatory finite volume central scheme. The proposed scheme adopts an inhomogeneous dual mesh with variable cell size ruled by the local wave speed of propagation to compute numerical fluxes at cell edges. In the limit of small time steps the central scheme gives rise to a semidiscrete formulation for the water saturation capable of incorporating heterogeneous porosity fields and generalized flux functions including the water transport due to the solid phase velocity. Numerical simulations of a water-flooding problem in secondary oil recovery are presented for different realizations of the input random fields (permeability, Young modulus and initial porosity). Comparison between the accuracies of the proposed approach and the traditional one-way coupled hydro-geomechanical formulation are presented. The effects of the cross-correlation between the input random fields and compaction drive mechanism upon finger growth and breakthrough curves are also analyzed. A notable feature of the formulation proposed herein is the accurate prediction of the influence of geomechanical effects upon the unstable movement of the water front, whose evolution is dictated by rock heterogeneity and unfavorable viscosity ratio, without deteriorating the local conservative character of the numerical schemes.     

Peridynamics simulations of geomaterial fragmentation by impulse loads

In this work, we have developed nonlinear peridynamics models of drained and saturated geomaterials, and applied them to simulations of dynamic fragmentation and ejecta formation due to impulse loads. First, we have re‐phrased and re‐interpreted the non‐local state‐based peridynamics formulation to connect the non‐local integral operator with the local differential operator. Second, we have implemented the Drucker–Prager plasticity model in state‐based peridynamics at finite strain. A peridynamics version of the Hughes and Winget algorithm is derived for the constitutive update. Third, we have developed a peridynamics U‐p formulation for saturated geomaterials. Fourth, numerical simulations have been carried out to verify the proposed peridynamics formulations in the simulation of geomaterial fragmentation induced by impulse loads. Comparisons of analytical and numerical results indicate that the peridynamics model has the ability to both match traditional continuum compression examples, as well as simulate complex geomaterial fragmentation processes resulting from impulse loads. Copyright © 2015 John Wiley & Sons, Ltd.     

Compositional flow modeling using a multi-point flux mixed finite element method

We present a general compositional formulation using multi-point flux mixed finite element (MFMFE) method on general hexahedral grids. The mixed finite element framework allows for local mass conservation, accurate flux approximation, and a more general treatment of boundary conditions. The multi-point flux inherent in MFMFE scheme allows the usage of a full permeability tensor. The proposed formulation is an extension of single and two-phase flow formulations presented by Wheeler and Yotov, SIAM J. Numer. Anal. 44(5), 2082–2106 (35) with similar convergence properties. Furthermore, the formulation allows for black oil, single-phase and multi-phase incompressible, slightly and fully compressible flow models utilizing the same design for different fluid systems. An accurate treatment of diffusive/dispersive fluxes owing to additional velocity degrees of freedom is also presented. The applications areas of interest include gas flooding, CO ₂ sequestration, contaminant removal, and groundwater remediation.     

Fully coupled generalized hybrid-mixed finite element approximation of two-phase two-component flow in porous media. Part I: formulation and properties of the mathematical model

This paper is a prequel to that of Marchand et al. (Comput Geosci 16:691–708, 2012), where an efficient and accurate hybrid-mixed finite element approximation for a system of time-dependent nonlinear conservation equations has been formulated, implemented, and tested, which are general enough to represent most of the existing formulations for two-component liquid–gas flow in porous medium with phase exchange, also allowing for any (dis)appearance of one of the phases. Temperature variation is neglected, but capillary effects are included by extended Darcy’s law, and Fickian diffusion is taken into account. The efficiency and stability of the numerical method of Lake (1989) relies on an equivalent reformulation of the otherwise commonly used model in terms of new principal variables and subsequent static (flash) equations allowing more generally for any (dis)appearance of one of the phases without the need of variable switching or unphysical quantities. In particular, the formulation in terms of complementarity conditions allows for an efficient and stable solution by the semismooth Newton’s method.     

Hybrid discretization of multi-phase flow in porous media in the presence of viscous,gravitational, and capillary forces

Multi-phase flow in porous media in the presence of viscous, gravitational, and capillary forces is described by advection diffusion equations with nonlinear parameters of relative permeability and capillary pressures. The conventional numerical method employs a fully implicit finite volume formulation. The phase-potential-based upwind direction is commonly used in computing the transport terms between two adjacent cells. The numerical method, however, often experiences non-convergence in a nonlinear iterative solution due to the discontinuity of transmissibilities, especially in transition between co-current and counter-current flows. Recently, Lee et al. (Adv. Wat. Res. 82, 27–38, 2015) proposed a hybrid upwinding method for the two-phase transport equation that comprises viscous and gravitational fluxes. The viscous part is a co-current flow with a one-point upwinding based on the total velocity and the buoyancy part is modeled by a counter-current flow with zero total velocity. The hybrid scheme yields C¹-continuous discretization for the transport equation and improves numerical convergence in the Newton nonlinear solver. Lee and Efendiev (Adv. Wat. Res. 96, 209–224, 2016) extended the hybrid upwind method to three-phase flow in the presence of gravity. In this paper, we present the hybrid-upwind formula in a generalized form that describes two- and three-phase flows with viscous, gravity, and capillary forces. In the derivation of the hybrid scheme for capillarity, we note that there is a strong similarity in mathematical formulation between gravity and capillarity. We thus greatly utilize the previous derivation of the hybrid upwind scheme for gravitational force in deriving that for capillary force. Furthermore, we also discuss some mathematical issues related to heterogeneous capillary domains and propose a simple discretization model by adapting multi-valued capillary pressures at the end points of capillary pressure curves. We demonstrate this new model always admits a consistent solution that is within the discretization error. This new generalized hybrid scheme yields a discretization method that improves numerical stability in reservoir simulation.     

Three‐dimensional analysis of the seismic behaviour of micropiles used in the reinforcement of saturated soil

This paper includes a numerical study of the behaviour of micropiles used for the reinforcement of saturated soil. Analysis is carried out using the (u–p) formulation (displacement for the solid phase and pore‐pressure for the fluid phase) implemented in a three‐dimensional finite element program. The soil behaviour is described by means of a cyclic elastoplastic constitutive relation which was developed within the framework of the bounding surface concept. The paper is composed of three parts. The first one is concerned with a presentation of the numerical model; the second includes analysis of the seismic behaviour of a single micropile; the last part deals with the group effect under seismic loading. Copyright © 2001 John Wiley & Sons, Ltd.     

A fully implicit and consistent finite element framework for modeling reservoir compaction with large deformation and nonlinear flow model. Part I: theory and formulation

Reservoir depletion results in rock failure, wellbore instability, hydrocarbon production loss, oil sand production, and ground surface subsidence. Specifically, the compaction of carbonate reservoirs with soft rocks often induces large plastic deformation due to rock pore collapse. On the other hand, following the compaction of reservoirs and failure of rock formations, the porosity and permeability of formations will, in general, decrease. These bring a challenge for reservoir simulations because of high nonlinearity of coupled geomechanics and fluid flow fields. In this work, we present a fully implicit, fully coupled, and fully consistent finite element formulation for coupled geomechanics and fluid flow problems with finite deformation and nonlinear flow models. The Pelessone smooth cap plasticity model, an important material model to capture rock compaction behavior and a challenging material model for implicit numerical formulations, is incorporated in the proposed formulation. Furthermore, a stress-dependent permeability model is taken into account in the formulation. A co-rotational framework is adopted for finite deformation, and an implicit material integrator for cap plasticity models is consistently derived. Furthermore, the coupled field equations are consistently linearized including nonlinear flow models. The physical theories, nonlinear material and flow models, and numerical formulations are the focus of part I of this work. In part II, we verify the proposed numerical framework and demonstrate the performance of our numerical formulation using several numerical examples including a field reservoir with soft rocks undergoing serious compaction.     

Limit Analysis of Unsaturated Soil Slope Stability Considering Intermediate Principal Stress and Strength Nonlinearity

Natural and man-made slopes are mostly in unsaturated soils, and the saturated or dry condition is just a particular case. This paper investigates the stability of unsaturated soil slopes. The influence of intermediate principal stress is considered by introducing the unified shear strength formulation, and the strength nonlinearity of unsaturated soils is addressed by means of two different methods. Based on the upper bound theorem of limit analysis, analytical expressions of the critical height of slope are deduced. Relevant numerical solutions are obtained through optimization and parametric analysis is then conducted. It is found that: the critical height of slope increases with the increase of UST parameter b or matric suction, while it decreases with the increase of slope angles. Besides, as regards the strength nonlinearity, it’s suggested to adopt Method II when parameter \(m\) can be reasonably determined.     

Stabilized single-point 4-node quadrilateral element for dynamic analysis of fluid saturated porous media

An accurate and efficient low-order quadrilateral mixed u?Cp element suitable for dynamic analysis of fluid saturated porous media is presented. The element uses physical hourglass stabilization to facilitate single-point integration for the solid phase, and non-residual stabilization of the fluid phase to circumvent instability in the incompressible-impermeable limit due to the use of equal-order interpolation for the displacement and pressure fields. Element behavior is verified and demonstrated through several numerical examples.     

Smeared crack approach: back to the original track

This paper briefly reviews the formulations used over the last 40 years for the solution of problems involving tensile cracking, with both the discrete and the smeared crack approaches. The paper focuses on the smeared approach, identifying as its main drawbacks the observed mesh‐size and mesh‐bias spurious dependence when the method is applied ‘straightly’. A simple isotropic local damage constitutive model is considered, and the (exponential) softening modulus is regularized according to the material fracture energy and the element size. The continuum and discrete mechanical problems corresponding to both the weak discontinuity (smeared cracks) and the strong discontinuity (discrete cracks) approaches are analysed and the question of propagation of the strain localization band (crack) is identified as the main difficulty to be overcome in the numerical procedure. A tracking technique is used to ensure stability of the solution, attaining the necessary convergence properties of the corresponding discrete finite element formulation. Numerical examples show that the formulation derived is stable and remarkably robust. As a consequence, the results obtained do not suffer from spurious mesh‐size or mesh‐bias dependence, comparing very favourably with those obtained with other fracture and continuum mechanics approaches. Copyright © 2006 John Wiley & Sons, Ltd.     

Explicit meshfree solution for large deformation dynamic problems in saturated porous media

In this paper a new methodology to simulate saturated soils subjected to dynamic loadings under large deformation regime (locally up to 40% in equivalent plastic strain) is presented. The coupling between solid and fluid phases is solved through the complete formulation of the Biot’s equations. The additional novelty lies in the employment of an explicit time integration scheme of the \(u-w\) (solid displacement–relative fluid displacement) formulation which enables us to take advantage of such explicit schemes. Shape functions based on the principle of maximum entropy implemented in the framework of Optimal Transportation Meshfree schemes are utilized to solve both elastic and plastic problems.     

Numerical modeling of three‐phase dissolution of underground cavities using a diffuse interface model

Natural evaporite dissolution in the subsurface can lead to cavities having critical dimensions in the sense of mechanical stability. Geomechanical effects may be significant for people and infrastructures because the underground dissolution may lead to subsidence or collapse (sinkholes). The knowledge of the cavity evolution in space and time is thus crucial in many cases. In this paper, we describe the use of a local nonequilibrium diffuse interface model for solving dissolution problems involving multimoving interfaces within three phases, that is, solid–liquid–gas as found in superficial aquifers and karsts. This paper generalizes developments achieved in the fluid–solid case, that is, the saturated case [1]. On one hand, a local nonequilibrium dissolution porous medium theory allows to describe the solid–liquid interface as a diffuse layer characterized by the evolution of a phase indicator (e.g., porosity). On the other hand, the liquid–gas interface evolution is computed using a classical porous medium two‐phase flow model involving a phase saturation, that is, generalized Darcy's laws. Such a diffuse interface model formulation is suitable for the implementation of a finite element or finite volume numerical model on a fixed grid without an explicit treatment of the interface movement. A numerical model has been implemented using a finite volume formulation with adaptive meshing (e.g., adaptive mesh refinement), which improves significantly the computational efficiency and accuracy because fine gridding may be attached to the dissolution front. Finally, some examples of three‐phase dissolution problems including density effects are also provided to illustrate the interest of the proposed theoretical and numerical framework. Copyright © 2014 John Wiley & Sons, Ltd.     

全长黏结式锚杆沿程应力分布模拟方法

全长黏结式锚杆是地下洞室常用的一种围岩支护手段。研究全长黏结式锚杆沿程应力分布特性对于锚杆-围岩组合体的稳定性评价具有重要的作用。利用单元嵌入计算格式,建立锚杆-围岩组合体中锚杆的等效模拟方法。通过轴对称四边形单元模拟锚杆内部变形,据此推导该单元内力、切线刚度及锚杆模型固有内力、计算刚度,最终通过整体及局部坐标的变换建立锚杆模型附加内力及切线刚度的计算方法。锚杆的非线性特性主要由接触面失效、杆体屈服引起,该方法能够较方便地应用隐式有限元计算格式,对围岩及锚杆的非线性特性进行统一分析。该方法在反映锚杆对围岩支护效果的同时,能够较好地反映锚杆的沿程应力分布特性。     

On persistent primary variables for numerical modeling of gas migration in a nuclear waste repository

Numerical simulation of gas migration driven by compressible two-phase partially miscible flow in porous media is of major importance for safety assessment of deep geological repositories for long-lived high-level nuclear waste. We present modeling of compositional liquid and gas flow for numerical simulations of hydrogen migration in deep geological radioactive waste repository based on persistent primary variables. Two-phase flow is considered, with incompressible liquid and compressible gas, which includes capillary effects, gas dissolution, and diffusivity. After discussing briefly the existing approaches to deal with phase appearance and disappearance problem, including a persistent set of variables already considered in a previous paper (Bourgeat et al., Comput Geosci 13(1):29–42, 2009), we focus on a new variant of the primary variables: dissolved hydrogen mass concentration and liquid pressure. This choice leads to a unique and consistent formulation in liquid saturated and unsaturated regions, which is well adapted to heterogeneous media. We use this new set of variable for numerical simulations and show computational evidences of its adequacy to simulate gas phase appearance and disappearance in different but typical situations for gas migration in an underground radioactive waste repository.     

A phase-field method for the direct simulation of two-phase flows in pore-scale media using a non-equilibrium wetting boundary condition

Advances in pore-scale imaging (e.g., μ-CT scanning), increasing availability of computational resources, and recent developments in numerical algorithms have started rendering direct pore-scale numerical simulations of multi-phase flow on pore structures feasible. Quasi-static methods, where the viscous and the capillary limit are iterated sequentially, fall short in rigorously capturing crucial flow phenomena at the pore scale. Direct simulation techniques are needed that account for the full coupling between capillary and viscous flow phenomena. Consequently, there is a strong demand for robust and effective numerical methods that can deliver high-accuracy, high-resolution solutions of pore-scale flow in a computationally efficient manner. Direct simulations of pore-scale flow on imaged volumes can yield important insights about physical phenomena taking place during multi-phase, multi-component displacements. Such simulations can be utilized for optimizing various enhanced oil recovery (EOR) schemes and permit the computation of effective properties for Darcy-scale multi-phase flows.We implement a phase-field model for the direct pore-scale simulation of incompressible flow of two immiscible fluids. The model naturally lends itself to the transport of fluids with large density and viscosity ratios. In the phase-field approach, the fluid-phase interfaces are expressed in terms of thin transition regions, the so-called diffuse interfaces, for increased computational efficiency. The conservation law of mass for binary mixtures leads to the advective Cahn–Hilliard equation and the condition that the velocity field is divergence free. Momentum balance, on the other hand, leads to the Navier–Stokes equations for Newtonian fluids modified for two-phase flow and coupled to the advective Cahn–Hilliard equation. Unlike the volume of fluid (VoF) and level-set methods, which rely on regularization techniques to describe the phase interfaces, the phase-field method facilitates a thermodynamic treatment of the phase interfaces, rendering it more physically consistent for the direct simulations of two-phase pore-scale flow. A novel geometric wetting (wall) boundary condition is implemented as part of the phase-field method for the simulation of two-fluid flows with moving contact lines. The geometric boundary condition accurately replicates the prescribed equilibrium contact angle and is extended to account for dynamic (non-equilibrium) effects. The coupled advective Cahn–Hilliard and modified Navier–Stokes (phase-field) system is solved by using a robust and accurate semi-implicit finite volume method. An extension of the momentum balance equations is also implemented for Herschel–Bulkley (non-Newtonian) fluids. Non-equilibrium-induced two-phase flow problems and dynamic two-phase flows in simple two-dimensional (2-D) and three-dimensional (3-D) geometries are investigated to validate the model and its numerical implementation. Quantitative comparisons are made for cases with analytical solutions. Two-phase flow in an idealized 2-D pore-scale conduit is simulated to demonstrate the viability of the proposed direct numerical simulation approach.