Robust numerical implementation of a 3D rate-dependent model for reservoir geomechanical simulations

A 3D elasto-plastic rate-dependent model for rock mechanics is formulated and implemented into a Finite Element (FE) numerical code. The model is based on the approach proposed by Vermeer and Neher (A soft soil model that accounts for creep. In: Proceedings of the International Symposium “Beyond 2000 in Computational Geotechnics,” pages 249-261, 1999). An original strain-driven algorithm with an Inexact Newton iterative scheme is used to compute the state variables for a given strain increment.The model is validated against laboratory measurements, checked on a simplified test case, and used to simulate land subsidence due to groundwater and hydrocarbon production. The numerical results prove computationally effective and robust, thus allowing for the use of the model on real complex geological settings.     

Upscaling of elastic properties for large scale geomechanical simulations

Large scale geomechanical simulations are being increasingly used to model the compaction of stress dependent reservoirs, predict the long term integrity of under‐ground radioactive waste disposals, and analyse the viability of hot‐dry rock geothermal sites. These large scale simulations require the definition of homogenous mechanical properties for each geomechanical cell whereas the rock properties are expected to vary at a smaller scale. Therefore, this paper proposes a new methodology that makes possible to define the equivalent mechanical properties of the geomechanical cells using the fine scale information given in the geological model. This methodology is implemented on a synthetic reservoir case and two upscaling procedures providing the effective elastic properties of the Hooke's law are tested. The first upscaling procedure is an analytical method for perfectly stratified rock mass, whereas the second procedure computes lower and upper bounds of the equivalent properties with no assumption on the small scale heterogeneity distribution. Both procedures are applied to one geomechanical cell extracted from the reservoir structure. The results show that the analytical and numerical upscaling procedures provide accurate estimations of the effective parameters. Furthermore, a large scale simulation using the homogenized properties of each geomechanical cell calculated with the analytical method demonstrates that the overall behaviour of the reservoir structure is well reproduced for two different loading cases. Copyright © 2004 John Wiley & Sons, Ltd.     

虚拟单元法求解渗流流量

提出了计算渗流场流量的虚拟单元法,推导了虚拟单元法有限元方程,并探讨了虚拟单元法的有限元实现方法,建立了虚拟单元法求解渗流流量的理论和方法体系。二维径向渗流和平底坝基算例比较分析表明,虚拟单元法具有计算精度高及有限元实现简单并容易编写通用程序等优点。虚拟单元法与中线法计算流量具有相同的精度,但在应用上虚拟单元法有一定优势。     

Boundary element analysis of geomechanical fracture

In this paper, a single-region BEM formulation for the three-dimensional analysis of fractures in geomechanics is developed. The technique allows the use of continuous elements in the discretization of the crack surfaces. Intially, an example from earthquake control theory is solved to demonstrated the validity of the technique for partially loaded crack surfaces. The method is then extended to deal with contact between crack surfaces by the incorporation of spring constraints. This allows two futher applications to be studied, one from hydraulic fracture and the other involving extraction of an ore seam. Copyright © 1999 John Wiley & Sons, Ltd.     

Ensemble smoothing of land subsidence measurements for reservoir geomechanical characterization

Geomechanical models are often used to predict the impact on land surface of fluid withdrawal from deep reservoirs, as well as investigating measures for mitigation. The ability to accurately simulate surface displacements, however, is often impaired by limited information on the geomechanical parameters characterizing the geological formations of interest. In this study, we employ an ensemble smoother, a data assimilation algorithm, to provide improved estimates of reservoir parameters through assimilation of measurements of both horizontal and vertical surface displacement into geomechanical model results. The method leverages the demonstrated potential of remote sensing techniques developed in the last decade to provide accurate displacement data for large areas of the land surface. For evaluation purposes, the methodology is applied to the case of a disk‐shaped reservoir embedded in a homogeneous, isotropic, and linearly elastic half space, subject to a uniform change in fluid pressure. Multiple sources of uncertainty are investigated, including the radius, R, the thickness, h, and the depth, c, of the reservoir; the pore pressure change, Δp; porous medium's vertical uniaxial compressibility, c_M, and Poisson's ratio, ν, and the ratio, s, between the compressibilities of the medium during loading and unloading cycles. Results from all simulations show that the ensemble smoother has the capability to effectively reduce the uncertainty associated with those parameters to which the variability and the spatial distribution of land surface displacements are most sensitive, namely, R, c, c_M, and s. These analyses demonstrate that the estimation of these parameters values depends on the number of measurements assimilated and the error assigned to the measurement values. Copyright © 2014 John Wiley & Sons, Ltd.     

Discrete element method simulations of bio-cemented sands

Microbially induced calcite precipitation (MICP) has emerged as a novel soil improvement method. In this paper, 3-D discrete element method (DEM) simulations are used to explore the behavior of MICP-cemented sands. Comparisons of the macro-scale response of numerical and physical specimens are made. Microstructure analyses indicate a shear band formed in the numerical specimens, consistent with physical experiments. The bond breakage pattern in numerical specimens is explored and compared to observed measurements from physical specimens. The relationship between dilatancy and stress-strain behavior is evaluated. The results indicate DEM is an effective technique to capture the mechanical behavior of MICP-cemented sand.     

Rock quality designation-fracture intensity index method for geomechanical classification

Rock quality designation (RQD) is a simple and effective way of rock mass classification from surface scanlines or preferably borehole measurements. A major drawback in its classical use is that only one representative RQD value is obtained from the field measurements per core run. However, it is shown in this paper that the field measurements along a scanline provide detailed information about the rock quality and the fracture intensity (FI) for a given number of joints. In order to be able to extract the complete information concealed within the field data, the RQD-fracture index diagram concept is proposed and applied to actual field scanline measurements from England. The basis of the suggested procedure is to obtain a series of all possible RQD and FI values from given intact length measurements. This procedure provides additional information as to the local zones of heterogeneities within the rock mass and a new way of rock classification on the basis of the combined effects of RQD and FI. It is also possible to calculate percentages of different rock qualities within the same rock mass.     

Application of the representer method for parameter estimation in numerical reservoir models

A data assimilation method was applied to estimate poorly known parameters (permeabilities) in a numerical reservoir model. Most variational methods for data assimilation are based on the assumption that the model is perfect except for the poorly known parameters. The representer method allows also for model errors, i.e. for uncertainties in the state variables (pressures and saturations). The method is based on minimizing a cost functional, assuming all the errors and parameters to be multivariate Gaussian random variables with given mean and covariances. The uncertain parameters and variables are expanded into a finite sum of basis functions called representers, and the gradients of the cost functional are obtained with an adjoint method. This approach gives an optimal parametrization in the sense that the final result is equal to the solution of the full inverse problem. The method was tested on a simple one-dimensional model to simulate two-phase (oil-water) flow through a heterogeneous reservoir. The results show that the method is able to provide an acceptable estimate of the permeability field. We used pressure measurements from a small number of observation wells in between the injection and production wells, but the representer method could be used equally well to assimilate data from other sources. The method appears to be a promising data assimilation tool for applications in reservoir engineering.     

Fast multipole displacement discontinuity method (FM‐DDM) for geomechanics reservoir simulations

The displacement discontinuity method (DDM) is frequently used in geothermal and petroleum applications for modeling the behavior of fractures in linear‐elastic rocks. The DDM requires O(N²) memory and O(N³) floating point operations (where N is the number of unknowns) to construct the coefficient matrix and solve the linear system of equations by direct methods. Therefore, the conventional implementation of the DDM is not computationally efficient for very large systems of cracks, often limiting its application to small‐scale problems. This work presents an approach for solving large‐scale fracture problems using the fast multipole method (FMM). The approach uses both the DDM and a kernel‐independent version of the FMM along with a preconditioned generalized minimal residual algorithm to accelerate the solution of linear systems of equations using desktop computers. Using the fundamental solutions for constant displacement discontinuity in a two‐dimensional elastic medium, several numerical examples involving fracture networks representing fractured reservoirs are treated. Numerical results show good agreement with analytical solutions and demonstrate the efficiency of the FMM implementation of the DDM for large‐scale simulations. Copyright © 2015 John Wiley & Sons, Ltd.     

Online conservative generalized multiscale finite element method for highly heterogeneous flow models

Computational Geosciences - In this work, we consider an online enrichment procedure in the context of the Generalized Multiscale Finite Element Method (GMsFEM) for the two-phase flow model in...     

Finite element methods for geothermal reservoir simulation

Two finite element algorithms suitable for long term simulation of geothermal reservoirs are presented. Both methods use a diagonal mass matrix and a Newton iteration scheme. The first scheme solves the 2N unsymmetric algebraic equations resulting from the finite element discretization of the equations governing the flow of heat and mass in porous media by using a banded equation solver. The second method, suitable for problems in which the transmissibility terms are small compared to the accumulation terms, reduces the set of N equations for the Newton corrections to a symmetric system. Comparison with finite difference schemes indicates that the proposed algorithms are competitive with existing methods.     

离散单元法在水库库岸滑坡稳定性分析中的应用

在现场勘察和室内试验的基础上,建立了离散单元法的物理力学模型,并且考虑到水库蓄水对地下水渗流场的影响,定量地研究了地下水动水压力对滑坡体稳定性的影响.把该模型应用于三峡库区秭归县下土地岭滑坡的稳定性评价中,计算结果表明,滑坡的失稳主要是由于水库蓄水的影响,滑坡体内的地下水渗流场发生了变化,在自重力、地下水动水压力的共同作用下,滑坡体的后缘和中部滑体首先出现剪切滑移,进而扩展到整个滑坡体,使滑坡产生整体失稳.     

Fabric symmetry of low anisotropic rocks inferred from ultrasonic sounding: Implications for the geomechanical models

The geomechanical models were established based on the absence or presence of certain rock fabric elements — texture (crystallographic preferred orientation), microstructure (shape preferred orientation) and microcracks (flat voids). The proposed models include both (i) the ideal material showing random texture and structure but no microcracks, i.e. the material which is hardly to be found in nature, and (ii) the materials possessing various combinations of fabric elements that show different spatial arrangements. The mutual relationship between those parameters and seismic and geomechanical properties are discussed.Selected models were experimentally verified during laboratory experiments. These consist of measurement of P-wave velocities in 132 independent directions under several confining pressures in the range 0.1–400 MPa. From measured data 3D P-wave patterns can be constructed and the influence of microcracks and of texture and structure on the rock seismic anisotropy can be determined. The seismic anisotropy established at different levels of confining pressure can be used for the interpretation of rock fabric symmetry of rocks showing low anisotropy in macroscale and for the selection of directions in which the geomechanical test can be performed. The measured P-wave velocities were then mathematically processed by using a fitting function which reflects contribution of P-wave velocity in the mineral skeleton of an ideal sample without microcracks extrapolated to the atmospheric pressure level from high confining pressure interval (ca. 200–400 MPa) (v₀), linear compressibility of the samples (k_v), and confining pressure during which most of the cracks are closed (P₀). These parameters improve the understanding of the response of various rock fabric elements on increasing confinement and corresponding changes in elasticity.The observed seismic and geomechanical anisotropies reflect intensity of the fabric of rock-forming minerals and microcracks. The magnitude of seismic anisotropy measured at atmospheric pressure corresponds to the anisotropy of static elastic modulus and is governed by the spatial arrangement of microcracks. The magnitude of strength anisotropy (uniaxial compressive strength) correlates more likely to the seismic anisotropy determined at high confining pressure and is connected to the preferred orientations (either CPO or SPO or both) of rock-forming minerals.     

Variational principles and finite element simulations for thermo-elastic consolidation

In this paper, a general variational principle for the initial boundary value problem of quasi-static thermoelastic consolidation is developed by assuming infinitesimal deformation and an incompressible fluid flowing through a linearly elastic solid. By manipulating the coupling operators, an extended form of the variational pronciple is derved. The associated finite element formulation based on this principle is presented and numerical applications for plane strain thermo-elastic consolidation are revealed.     

FV-MHMM method for reservoir modeling

The present paper proposes a new family of multiscale finite volume methods. These methods usually deal with a dual mesh resolution, where the pressure field is solved on a coarse mesh, while the saturation fields, which may have discontinuities, are solved on a finer reservoir grid, on which petrophysical heterogeneities are defined. Unfortunately, the efficiency of dual mesh methods is strongly related to the definition of up-gridding and down-gridding steps, allowing defining accurately pressure and saturation fields on both fine and coarse meshes and the ability of the approach to be parallelized. In the new dual mesh formulation we developed, the pressure is solved on a coarse grid using a new hybrid formulation of the parabolic problem. This type of multiscale method for pressure equation called multiscale hybrid-mixed method (MHMM) has been recently proposed for finite elements and mixed-finite element approach (Harder et al. 2013). We extend here the MH-mixed method to a finite volume discretization, in order to deal with large multiphase reservoir models. The pressure solution is obtained by solving a hybrid form of the pressure problem on the coarse mesh, for which unknowns are fluxes defined on the coarse mesh faces. Basis flux functions are defined through the resolution of a local finite volume problem, which accounts for local heterogeneity, whereas pressure continuity between cells is weakly imposed through flux basis functions, regarded as Lagrange multipliers. Such an approach is conservative both on the coarse and local scales and can be easily parallelized, which is an advantage compared to other existing finite volume multiscale approaches. It has also a high flexibility to refine the coarse discretization just by refinement of the lagrange multiplier space defined on the coarse faces without changing nor the coarse nor the fine meshes. This refinement can also be done adaptively w.r.t. a posteriori error estimators. The method is applied to single phase (well-testing) and multiphase flow in heterogeneous porous media.     

Empirical assessment of the critical time increment in explicit particulate discrete element method simulations

This contribution considers the critical time increment (Δt_crit) to achieve stable simulations using particulate discrete element method (DEM) codes that adopt a Verlet-type time integration scheme. The Δt_crit is determined by considering the maximum vibration frequency of the system. Based on a series of parametric studies, Δt_crit is shown to depend on the particle mass (m), the maximum contact stiffness (K_max), and the maximum particle coordination number (C_N,max). Empirical expressions relating Δt_crit to m, K_max, and C_N,max are presented; while strictly only valid within the range of simulation scenarios considered here, these can inform DEM analysts selecting appropriate Δt_crit values.     

基于强度折减法的库岸滑坡三维有限元分析

通过现场踏勘调研,根据滑坡的环境地质条件,详细探讨了该库岸滑坡的成因机制和演化过程,并重点分析了滑坡失稳的影响因素。建立滑坡体的三维有限元模型,运用有限元强度折减理论对其稳定性进行了计算分析。结果表明,滑坡体在天然状况下是稳定的,库水位上升至正常蓄水位时将超过滑坡前缘高程,滑坡安全系数降低,不满足规范规定的最小安全系数的要求,在考虑降雨入渗、库水位骤降时,滑坡进一步趋于危险;在考虑蓄水+地震工况时,库岸滑坡有失稳的可能。据此提出了削坡减载治理措施,治理后滑坡在各工况下的安全系数均满足规范的要求,研究成果对滑坡的治理具有较大的参考价值。     

A stabilized extended finite element framework for hydraulic fracturing simulations

We present a stabilized extended finite element formulation to simulate the hydraulic fracturing process in an elasto‐plastic medium. The fracture propagation process is governed by a cohesive fracture model, where a trilinear traction‐separation law is used to describe normal contact, cohesion and strength softening on the fracture face. Fluid flow inside the fracture channel is governed by the lubrication equation, and the flow rate is related to the fluid pressure gradient by the ‘cubic’ law. Fluid leak off happens only in the normal direction and is assumed to be governed by the Carter's leak‐off model. We propose a ‘local’ U‐P (displacement‐pressure) formulation to discretize the fluid‐solid coupled system, where volume shape functions are used to interpolate the fluid pressure field on the fracture face. The ‘local’ U‐P approach is compatible with the extended finite element framework, and a separate mesh is not required to describe the fluid flow. The coupled system of equations is solved iteratively by the standard Newton‐Raphson method. We identify instability issues associated with the fluid flow inside the fracture channel, and use the polynomial pressure projection method to reduce the pressure oscillations resulting from the instability. Numerical examples demonstrate that the proposed framework is effective in modeling 3D hydraulic fracture propagation. Copyright © 2016 John Wiley & Sons, Ltd.