Numerical Homogenization of the Rigidity Tensor in Hooke''s Law Using the Node-Based Finite Element Method

Combining a geological model with a geomechanical model, it generally turns out that the geomechanical model is built from units that are at least a 100 times larger in volume than the units of the geological model. To counter this mismatch in scales, the geological data model's heterogeneous fine-scale Young's moduli and Poisson's ratios have to be upscaled to one equivalent homogeneous coarse-scale rigidity. This coarse-scale rigidity relates the volume-averaged displacement, strain, stress, and energy to each other, in such a way that the equilibrium equation, Hooke's law, and the energy equation preserve their fine-scale form on the coarse scale. Under the simplifying assumption of spatial periodicity of the heterogeneous fine-scale rigidity, homogenization theory can be applied. However, even then the spatial variability is generally so complex that exact solutions cannot be found. Therefore, numerical approximation methods have to be applied. Here the node-based finite element method for the displacement as primary variable has been used. Three numerical examples showing the upper bound character of this finite element method are presented.     

Nodal and Mixed Finite Elements for the Numerical Homogenization of 3D Permeability

The aim of upscaling is to determine equivalent homogeneous parameters at a coarse-scale from a spatially oscillating fine-scale parameter distribution. To be able to use a limited number of relatively large grid-blocks in numerical oil reservoir simulators or groundwater models, upscaling of the permeability is frequently applied. The spatial fine-scale permeability distribution is generally obtained from geological and geostatistical models. After upscaling, the coarse-scale permeabilities are incorporated in the relatively large grid-blocks of the numerical model. If the porous rock may be approximated as a periodic medium, upscaling can be performed by the method of homogenization. In this paper the homogenization is performed numerically, which gives rise to an approximation error. The complementarity between two different numerical methods – the conformal-nodal finite element method and the mixed-hybrid finite element method – has been used to quantify this error. These two methods yield respectively upper and lower bounds for the eigenvalues of the coarse-scale permeability tensor. Results of 3D numerical experiments are shown, both for the far field and around wells.     

An Approach for Modeling Rock Discontinuous Mechanical Behavior Under Multiphase Fluid Flow Conditions

In this paper, the two computer codes TOUGH2 and RDCA (for “rock discontinuous cellular automaton”) are integrated for coupled hydromechanical analysis of multiphase fluid flow and discontinuous mechanical behavior in heterogeneous rock. TOUGH2 is a well-established code for geohydrological analysis involving multiphase, multicomponent fluid flow and heat transport; RDCA is a numerical model developed for simulating the nonlinear and discontinuous geomechanical behavior of rock. The RDCA incorporates the discontinuity of a fracture independently of the mesh, such that the fracture can be arbitrarily located within an element, while the fluid pressure calculated by TOUGH2 can be conveniently applied to fracture surfaces. We verify and demonstrate the coupled TOUGH–RDCA simulator by modeling a number of simulation examples related to coupled multiphase flow and geomechanical processes associated with the deep geological storage of carbon dioxide—including modeling of ground surface uplift, stress-dependent permeability, and the coupled multiphase flow and geomechanical behavior of fractures intersecting the caprock.     

Upscaling Hydraulic Conductivities in Cross-Bedded Formations

Modern geostatistical techniques allow the generation of high-resolution heterogeneous models of hydraulic conductivity containing millions to billions of cells. Selective upscaling is a numerical approach for the change of scale of fine-scale hydraulic conductivity models into coarser scale models that are suitable for numerical simulations of groundwater flow and mass transport. Selective upscaling uses an elastic gridding technique to selectively determine the geometry of the coarse grid by an iterative procedure. The geometry of the coarse grid is built so that the variances of flow velocities within the coarse blocks are minimum. Selective upscaling is able to handle complex geological formations and flow patterns, and provides full hydraulic conductivity tensor for each block. Selective upscaling is applied to a cross-bedded formation in which the fine-scale hydraulic conductivities are full tensors with principal directions not parallel to the statistical anisotropy of their spatial distribution. Mass transport results from three coarse-scale models constructed by different upscaling techniques are compared to the fine-scale results for different flow conditions. Selective upscaling provides coarse grids in which mass transport simulation is in good agreement with the fine-scale simulations, and consistently superior to simulations on traditional regular (equal-sized) grids or elastic grids built without accounting for flow velocities.     

Multiscale finite-volume method for density-driven flow in porous media

The multiscale finite-volume (MSFV) method has been developed to solve multiphase flow problems on large and highly heterogeneous domains efficiently. It employs an auxiliary coarse grid, together with its dual, to define and solve a coarse-scale pressure problem. A set of basis functions, which are local solutions on dual cells, is used to interpolate the coarse-grid pressure and obtain an approximate fine-scale pressure distribution. However, if flow takes place in presence of gravity (or capillarity), the basis functions are not good interpolators. To treat this case correctly, a correction function is added to the basis function interpolated pressure. This function, which is similar to a supplementary basis function independent of the coarse-scale pressure, allows for a very accurate fine-scale approximation. In the coarse-scale pressure equation, it appears as an additional source term and can be regarded as a local correction to the coarse-scale operator: It modifies the fluxes across the coarse-cell interfaces defined by the basis functions. Given the closure assumption that localizes the pressure problem in a dual cell, the derivation of the local problem that defines the correction function is exact, and no additional hypothesis is needed. Therefore, as in the original MSFV method, the only closure approximation is the localization assumption. The numerical experiments performed for density-driven flow problems (counter-current flow and lock exchange) demonstrate excellent agreement between the MSFV solutions and the corresponding fine-scale reference solutions.     

Statistical assignment of upscaled flow functions for an ensemble of geological models

Upscaled flow functions are often needed to account for the effects of fine-scale permeability heterogeneity in coarse-scale simulation models. We present procedures in which the required coarse-scale flow functions are statistically assigned to an ensemble of upscaled geological models. This can be viewed as an extension and further development of a recently developed ensemble level upscaling (EnLU) approach. The method aims to efficiently generate coarse-scale flow models capable of reproducing the ensemble statistics (e.g., cumulative distribution function) of fine-scale flow predictions for multiple reservoir models. The most expensive part of standard coarsening procedures is typically the generation of upscaled two-phase flow functions (e.g., relative permeabilities). EnLU provides a means for efficiently generating these upscaled functions using stochastic simulation. This involves the use of coarse-block attributes that are both fast to compute and correlate closely with the upscaled two-phase functions. In this paper, improved attributes for use in EnLU, namely the coefficient of variation of the fine-scale single-phase velocity field (computed during computation of upscaled absolute permeability) and the integral range of the fine-scale permeability variogram, are identified. Geostatistical simulation methods, which account for spatial correlations of the statistically generated upscaled functions, are also applied. The overall methodology thus enables the efficient generation of coarse-scale flow models. The procedure is tested on 3D well-driven flow problems with different permeability distributions and variable fluid mobility ratios. EnLU is shown to capture the ensemble statistics of fine-scale flow results (water and oil flow rates as a function of time) with similar accuracy to full flow-based upscaling methods but with computational speedups of more than an order of magnitude.     

Successful application of multiscale methods in a real reservoir simulator environment

For the past 10 years or so, a number of so-called multiscale methods have been developed as an alternative approach to upscaling and to accelerate reservoir simulation. The key idea of all these methods is to construct a set of prolongation operators that map between unknowns associated with cells in a fine grid holding the petrophysical properties of the geological reservoir model and unknowns on a coarser grid used for dynamic simulation. The prolongation operators are computed numerically by solving localized flow problems, much in the same way as for flow-based upscaling methods, and can be used to construct a reduced coarse-scale system of flow equations that describe the macro-scale displacement driven by global forces. Unlike effective parameters, the multiscale basis functions have subscale resolution, which ensures that fine-scale heterogeneity is correctly accounted for in a systematic manner. Among all multiscale formulations discussed in the literature, the multiscale restriction-smoothed basis (MsRSB) method has proved to be particularly promising. This method has been implemented in a commercially available simulator and has three main advantages. First, the input grid and its coarse partition can have general polyhedral geometry and unstructured topology. Secondly, MsRSB is accurate and robust when used as an approximate solver and converges relatively fast when used as an iterative fine-scale solver. Finally, the method is formulated on top of a cell-centered, conservative, finite-volume method and is applicable to any flow model for which one can isolate a pressure equation. We discuss numerical challenges posed by contemporary geomodels and report a number of validation cases showing that the MsRSB method is an efficient, robust, and versatile method for simulating complex models of real reservoirs.     

Multiscale finite-volume formulation for multiphase flow in porous media: black oil formulation of compressible,three-phase flow with gravity

Most practical reservoir simulation studies are performed using the so-called black oil model, in which the phase behavior is represented using solubilities and formation volume factors. We extend the multiscale finite-volume (MSFV) method to deal with nonlinear immiscible three-phase compressible flow in the presence of gravity and capillary forces (i.e., black oil model). Consistent with the MSFV framework, flow and transport are treated separately and differently using a sequential implicit algorithm. A multiscale operator splitting strategy is used to solve the overall mass balance (i.e., the pressure equation). The black-oil pressure equation, which is nonlinear and parabolic, is decomposed into three parts. The first is a homo geneous elliptic equation, for which the original MSFV method is used to compute the dual basis functions and the coarse-scale transmissibilities. The second equation accounts for gravity and capillary effects; the third equation accounts for mass accumulation and sources/ sinks (wells). With the basis functions of the elliptic part, the coarse-scale operator can be assembled. The gravity/capillary pressure part is made up of an elliptic part and a correction term, which is computed using solutions of gravity-driven local problems. A particular solution represents accumulation and wells. The reconstructed fine-scale pressure is used to compute the fine-scale phase fluxes, which are then used to solve the nonlinear saturation equations. For this purpose, a Schwarz iterative scheme is used on the primal coarse grid. The framework is demonstrated using challenging black-oil examples of nonlinear compressible multiphase flow in strongly heterogeneous formations.     

Cell to node projections: An assessment of error

Reservoir simulators typically use cell‐centered finite volume schemes and do not model directly the coupling of the flow processes with the geomechanics. Coupling of geomechanics with fluid flow can be important in many cases, but introducing fully coupled geomechanical effects in those simulators is not a trivial issue, because the geomechanics is better done by using the Galerkin vertex‐centered finite element methods by which the solid displacements are computed at the vertices of the cells. This creates difficulties in interfacing cell variables with nodal variables. Uncoupled or loosely coupled models are used by many researchers/practitioners by which a reservoir model is coupled to a geomechanical model by staggering in‐time flow and deformation via a sophisticated interface that repeatedly calls first flow and then mechanics. The method therefore requires projection of the reservoir cell variables onto the nodes of the geomechanics Galerkin finite element mesh. In this note, we attempt to quantify the errors associated with cell to node projection operations. For that purpose, we use a simple model of the pressure equation for a heterogeneous medium in one dimension. We are able to derive the exact analytical solution for this problem for both nodal and cell pressures. This allows us to compute the errors due to projection analytically, function of meshing refinement and permeability field variations. We compute upper and lower bounds for the errors, and analyze their magnitude for a variety of cases. We conclude that, in general, cell to node projection operations lead to substantial errors. Copyright © 2010 John Wiley & Sons, Ltd.     

Using Multiple-Point Geostatistics for Tracer Test Modeling in a Clay-Drape Environment with Spatially Variable Conductivity and Sorption Coefficient

This study investigates the effect of fine-scale clay drapes on tracer transport. A tracer test was performed in a sandbar deposit consisting of cross-bedded sandy units intercalated with many fine-scale clay drapes. The heterogeneous spatial distribution of the clay drapes causes a spatially variable hydraulic conductivity and sorption coefficient. A fluorescent tracer (sodium naphthionate) was injected in two injection wells and ground water was sampled and analyzed from five pumping wells. To determine (1) whether the fine-scale clay drapes have a significant effect on the measured concentrations and (2) whether application of multiple-point geostatistics can improve interpretation of tracer tests in media with complex geological heterogeneity, this tracer test is analyzed with a local three-dimensional ground-water flow and transport model in which fine-scale sedimentary heterogeneity is modeled using multiple-point geostatistics. To reduce memory needs and calculation time for the multiple-point geostatistical simulation step, this study uses the technique of direct multiple-point geostatistical simulation of edge properties. Instead of simulating pixel values, model cell edge properties indicating the presence of irregularly shaped surfaces are simulated using multiple-point geostatistical simulations. Results of a sensitivity analysis show under which conditions clay drapes have a significant effect on the concentration distribution. Calibration of the model against measured concentrations from the tracer tests reduces the uncertainty on the clay-drape parameters. The calibrated model shows which features of the breakthrough curves can be attributed to the geological heterogeneity of the aquifer and which features are caused by other processes.     

模拟三维裂纹问题的多尺度扩展有限元法

扩展有限元法模拟裂纹时独立于网格,因此该方法是目前求解裂纹问题最有效的数值方法。为了在计算代价不大的情况,实现大型结构分析中考虑小裂纹或提高裂纹附近精度,在裂纹附近一般采用小尺度单元,其他区域采用大尺度单元。提出了分析三维裂纹问题的多尺度扩展有限元法,在需要的地方采用小尺度单元。基于点插值构造了六面体任意节点单元。所有尺度单元都采用8节点六面体单元,这样六面体任意节点单元可方便有效地连接不同尺度单元。采用互作用积分法计算三维应力强度因子。边裂纹和中心圆裂纹算例分析结果表明,该方法是正确和有效的。     

Kalman filter finite element method applied to dynamic ground motion

The purpose of this paper is to investigate the estimation of dynamic elastic behavior of the ground using the Kalman filter finite element method. In the present paper, as the state equation, the balance of stress equation, the strain–displacement equation and the stress–strain equation are used. For temporal discretization, the Newmark ¼ method is employed, and for the spatial discretization the Galerkin method is applied. The Kalman filter finite element method is a combination of the Kalman filter and the finite element method. The present method is adaptable to estimations not only in time but also in space, as we have confirmed by its application to the Futatsuishi quarry site. The input data are the measured velocity, acceleration, etc., which may include mechanical noise. It has been shown in numerical studies that the estimated velocity, acceleration, etc., at any other spatial and temporal point can be obtained by removing the noise included in the observation. Copyright © 2008 John Wiley & Sons, Ltd.     

A Combined Three-Dimensional Geological-Geostatistical-Numerical Model of Underground Excavations in Rock

Summary. This paper exploits geological and borehole geotechnical data obtained in the exploratory phase of a tunneling project to investigate in a first place if the kriging interpolation scheme may effectively reproduce the spatial variability of rock mass quality (Rock Mass Rating, RMR) in the vicinity of tunnels. For this purpose a quick solver in Fortran has been developed that performs variography analysis of 3D spatial data, fast kriging estimations of RMR between borehole sampling locations at the centroids of the elements of the numerical model, and model validation. For the purpose of an integrated underground excavation design, a step further is made by incorporating into the 3D mechanical numerical model of the rock mass, the three-dimensional (3D) solid geological model, thus coupling the geology with the ground (geotechnical) model (i.e. each element of the numerical model is assigned a geological material). The mechanical properties of each finite difference cell (or Representative Elementary Volume) of the ground model were then prescribed according to its geological type, the spatial heterogeneity of the rock mass expressed quantitatively with the kriging model, and the upscaling calculations of the mechanical properties of the intact rocks determined in the laboratory, based on the size-effect (strength dependence on size) and Damage Theory. Furthermore, a preliminary numerical simulation of the advance of unsupported tunnels in the model of the heterogeneous rock mass was performed for illustration purposes.     

Discontinuous analysis of softening solids under impact loading

Softening solids are analysed under impact loading using a new numerical method which allows displacement discontinuities to propagate arbitrarily through a finite element mesh. The Dirac‐delta distributions that arise in the strain field of classical continuum theory in the presence of strain softening are interpreted as discontinuities in the displacement field. A new finite element procedure with Heaviside jumps added to the underlying displacement interpolation basis is able to capture displacement jumps independent of the spatial discretisation. The amplitudes of displacement jumps are represented by extra degrees of freedom at existing nodes. Numerical results for mode‐I and mode‐II failure due to impact loading are presented. The numerical results highlight the objectivity of the approach with respect to spatial discretisation under dynamic loading conditions. Copyright © 2001 John Wiley & Sons, Ltd.     

Numerical modelling of regional faults in land subsidence prediction above gas/oil reservoirs

The stress variation induced by gas/oil production may activate pre‐existing regional faults. This may enhance the expected land subsidence due to the generation of mechanically weak points close to the producing field. A class of elasto‐plastic interface elements (IE), specifically designed to address the mechanical behaviour of faults over a regional scale, is integrated into a finite element (FE) geomechanical model and used to investigate the role exerted by active faults in anthropogenic land subsidence. The importance of regional faults depends on a variety of factors including depth of the depleted reservoir, fault number, orientation and size, geomechanical properties of porous medium, pore pressure drawdown induced by fluid production, etc. With the aid of some representative examples, a useful indication is provided as to where and how fault activation may influence both magnitude and extent of the land subsidence bowl above producing gas/oil reservoirs, pointing to a generally limited impact on the ground surface. The simulation of a real faulted gas reservoir in a complex 3‐D setting shows that the proposed IE can be simply and efficiently incorporated into a FE geomechanical model, thus improving the quality of the stress and displacement prediction. Copyright © 2007 John Wiley & Sons, Ltd.     

非均匀多孔介质等效渗透率的普适表达式

多孔介质渗流是普遍的物理过程,涉及地下工程、地热开采、环境工程等各行各业,尤其是工程建设,常面临防渗问题。由于地质条件的复杂性,工程区域地层受到成岩、压实、风化、生物作用等各种影响,故渗流性质复杂,常需要对建设区域的渗流状况进行数值模拟,从而为工程的设计施工提供决策依据。数值仿真结果依赖于对地层介质关键参数的选取,但目前工程多将其视为均匀介质处理,对于介质的非均匀特性考虑较少。文章旨在研究非均质多孔介质渗透率空间分布与等效渗透率的关系。基于连续介质假定、达西定律以及非均匀多孔介质渗透率空间分布函数,建立一维到三维的达西渗流问题模型,通过求解偏微分方程和理论推导,得到基于渗透率空间分布函数的等效渗透率理论表达式,并与有限元计算的数值解进行对比分析,结果表明理论值和数值解误差很小,证明等效渗透率的表达式的合理性。利用该成果可通过多点局部渗透率的测定构建渗透率空间分布函数,从而对整体渗流区域的渗透性质进行快速计算和评估,从而简化异常复杂的工程地质模型以减少计算量需求,对于工程仿真的快速计算和结果评估有重要意义。     

Quantitative Geoscience and Geological Big Data Development: A Review

After long-term development, mathematical geology has today become an independent discipline. Big Data science, which has become a new scientific paradigm in the 21 st century, gives rise to the geological Big Data, i.e. mathematical geology and quantitative geoscience. Thanks to a robust macro strategy for big data, China’s quantitative geoscience and geological big data’s rapid development meets present requirements and has kept up with international levels. This paper presents China’s decade-long achievements in quantitative prediction and assessment of mineral resources, geoscience information and software systems, geological information platform development, etc., with an emphasis on application of geological big data in informatics, quantitative mineral prediction, geological environment and disaster management, digital land survey, digital city, etc. Looking ahead, mathematical geology is moving towards "Digital Geology", "Digital Land" and "Geological Cloud", eventually realizing China’s grand "Digital China" blueprint, and these valuable results will be showcased on the international academic arena.     

Main active faults and seismic activity along the Yangtze River Economic Belt: Based on remote sensing geological survey

The Yangtze River Economic Belt (YREB) spans three terrain steps in China and features diverse topography that is characterized by significant differences in geological structure and present-day crustal deformation. Active faults and seismic activity are important geological factors for the planning and development of the YREB. In this paper, the spatial distribution and activity of 165 active faults that exist along the YREB have been compiled from previous findings, using both remote-sensing data and geological survey results. The crustal stability of seven particularly noteworthy typical active fault zones and their potential effects on the crustal stability of the urban agglomerations are analyzed. The main active fault zones in the western YREB, together with the neighboring regional active faults, make up an arc fault block region comprising primarily of Sichuan-Yunnan and a “Sichuan-Yunnan arc rotational-shear active tectonic system” strong deformation region that features rotation, shear and extensional deformation. The active faults in the central-eastern YREB, with seven NE-NNE and seven NW-NWW active faults (the “7-longitudinal, 7-horizontal” pattern), macroscopically make up a “chessboard tectonic system” medium-weak deformation region in the geomechanical tectonic system. They are also the main geological constraints for the crustal stability of the YREB.