Higher-order compositional modeling with Fickian diffusion in unstructured and anisotropic media

We present advances in compositional modeling of two-phase multi-component flow through highly complex porous media. Higher-order methods are used to approximate both mass transport and the velocity and pressure fields. We employ the Mixed Hybrid Finite Element (MHFE) method to simultaneously solve, to the same order, the pressure equation and Darcy's law for the velocity. The species balance equation is approximated by the discontinuous Galerkin (DG) approach, combined with a slope limiter. In this work we present an improved DG scheme where phase splitting is analyzed at all element vertices in the two-phase regions, rather than only as element averages. This approximation is higher-order than the commonly employed finite volume method and earlier DG approximations. The method reduces numerical dispersion, allowing for an accurate capture of shock fronts and lower dependence on mesh quality and orientation. Further new features are the extension to unstructured grids and support for arbitrary permeability tensors (allowing for both scalar heterogeneity, and shear anisotropy). The most important advancement in this work is the self-consistent modeling of two-phase multi-component Fickian diffusion. We present several numerical examples to illustrate the powerful features of our combined MHFE–dg method with respect to lower-order calculations, ranging from simple two component fluids to more challenging real problems regarding CO₂ injection into a vertical domain saturated with a multi-component petroleum fluid.     

Analysis of hydrodynamic conditions in adjacent free and heterogeneous porous flow domains

The existence of a free‐flow domain (e.g. a liquid layer) adjacent to a porous medium is a common occurrence in many environmental and petroleum engineering problems. The porous media may often contain various forms of heterogeneity, e.g. layers, fractures, micro‐scale lenses, etc. These heterogeneities affect the pressure distribution within the porous domain. This may influence the hydrodynamic conditions at the free–porous domain interface and, hence, the combined flow behaviour. Under steady‐state conditions, the heterogeneities are known to have negligible effects on the coupled flow behaviour. However, the significance of the heterogeneity effects on coupled free and porous flow under transient conditions is not certain. In this study, numerical simulations have been carried out to investigate the effects of heterogeneous (layered) porous media on the hydrodynamics conditions in determining the behaviour of combined free and porous regimes. Heterogeneity in the porous media is introduced by defining a domain composed of two layers of porous media with different values of intrinsic permeability. The coupling of the governing equations of motion in free and porous domains has been achieved through the well‐known Beavers and Joseph interfacial condition. Of special interest in this work are porous domains with flow‐through ends. They represent the general class of problems where large physical domains are truncated to smaller sections for ease of mathematical analysis. However, this causes a practical difficulty in modelling such systems. This is because the information on flow behaviour, i.e. boundary conditions at the truncated sections, is usually not available. Use of artificial boundary conditions to solve these problems effectively implies the imposition of conditions that do not necessarily match with the solutions required for the interior of the domain. This difficulty is resolved in this study by employing ‘stress‐free boundary conditions’ at the open ends of the domains, which have been shown to provide accurate results by a number of previous workers. Copyright © 2005 John Wiley & Sons, Ltd.     

Implementation of a second‐order paraxial boundary condition for a water‐saturated layered half‐space in plane strain

A half‐space finite element and a transmitting boundary are developed for a water‐saturated layered half‐space using a paraxial boundary condition. The exact dynamic stiffness of a half‐space in plane strain is derived and a second‐order paraxial approximation of the stiffness is obtained. A half‐space finite element and a transmitting boundary are then formulated. The development is verified by comparison of the dynamic stiffness of impermeable and permeable rigid strip foundations with other published results. The advantage of using the paraxial boundary condition in comparison with the rigid boundary condition is examined. It is shown that the paraxial boundary condition offers significant gain and the resulting half‐space finite element and transmitting boundary can represent the effects of a water‐saturated layered half‐space with good accuracy and efficiency. In addition, the numerical method described herein maintains the strengths and advantages of the finite element method and can be easily applied to demanding problems of soil–structure interaction in a water‐saturated layered half‐space. Copyright © 2010 John Wiley & Sons, Ltd.     

水平分层孔隙介质中点源激发的震电波场数值模拟及分析

为认识地震波诱导的电磁场的特性,本文研究地震波在孔隙介质中由于动电效应引起的电磁场.基于Pride弹性-电磁耦合方程组推导了双力偶震源对应的位移-应力-电磁场间断向量的表达式,模拟了双力偶源激发的震电波场.作为比较,还模拟了爆炸点源激发的震电波场.结果表明：存在伴随纵波的电场,其各分量的波形与固相位移对应分量的波形相似,但相位相反;存在伴随横波的磁场,其波形与固相位移波形相似;双力偶震源还激发出了独立传播的辐射电磁场,其速度比纵波至少高一个数量级,几乎是“瞬间同时”到达了每个接收器,但是其强度比伴随电磁场小得多,且随着源距增大而迅速减小.本文研究还表明：伴随纵波的电场强度不仅与地震纵波幅度和动电耦合系数有关,还与由介质孔隙结构决定的流-固两相动力协调性有关,存在一种动力协调介质,纵波在这种介质中不引起电场.     

A non-reflecting boundary condition for the finite element modeling of infinite reservoir with layered sediment

The design of seismic resistant concrete gravity dam necessitates accurate determination of hydrodynamic pressure developed in the adjacent reservoir. The hydrodynamic pressure developed on structure is dependent on the physical characteristics of the boundaries surrounding the reservoir including reservoir bottom. The sedimentary material in the reservoir bottom absorbs energy at the bottom, which will affect the hydrodynamic pressure at the upstream face of the dam. The fundamental parameter characterizing the effect of absorption of hydrodynamic pressure waves at the reservoir bottom due to sediment is the reflection coefficient. The wave reflection coefficient is determined from parameters based on sediment layer thickness, its material properties and excitation frequencies. An analytical or a closed-form solution cannot account for the arbitrary geometry of the dam or reservoir bed profile. This problem can be efficiently tackled with finite element technique. The need for an accurate truncation boundary is felt to reduce the computational domain of the unbounded reservoir system. An efficient truncation boundary condition (TBC) which accounts for the reservoir bottom effect is proposed for the finite element analysis of infinite reservoir. The results show the efficiency of the proposed truncation boundary condition.     

饱和多孔介质三维时域黏弹性人工边界与动力反应分析的显式有限元法

本文基于Biot的饱和多孔介质本构方程,考察具有辐射阻尼的外行球面波,推导了饱和多孔介质三维黏弹性人工边界的法向和切向边界方程;在已有的饱和多孔介质二维显式有限元数值计算方法基础上,提出该理论的三维方法,并开发了实现该三维方法的有限元程序.算例表明饱和多孔介质三维时域黏弹性人工边界与动力反应分析的显式有限元法具有较好的精度和稳定性.     

含孔隙介质的分层半空间表面瑞利波的衰减特性

以分层半空间内部含有一层孔隙介质为物理模型进行数值计算,研究半空间表面瑞利波的传播和衰减特性.为更加接近实际,结合瑞利波的激发特性,确定了瑞利波的主衰减曲线,并主要以此进行规律分析.针对速度递增和含低速层这两种典型的地质模型,讨论了瑞利波的传播衰减在不同地质模型下的特性,并分析了各自的规律.结果表明,在这两种模型下瑞利波的主衰减曲线都受孔隙介质所处空间位置影响产生比较明显的变化,但衰减系数极大值对应的波长与模型的表层厚度存在较明显的线性对应关系,利用这一关系,可以在实际勘探中快速得到表层介质厚度.另外,通过对比分析还可以看到,瑞利波主衰减曲线随孔隙介质的孔隙度和渗透率的变化都强于主频散曲线的变化,表明衰减曲线对孔隙度和渗透率的变化更加敏感,理论上更加适合进行介质参数反演工作.综合对比结果,我们认为瑞利波主衰减曲线中包含了更丰富的介质参数信息,如果能够有效利用,将可以提高瑞利波勘探的准确性和应用范围.

     

Wave-induced dynamic response of saturated multi-layer porous media: Analytical solutions and validity regions of various formulations in non-dimensional parametric space

In this paper, dynamic response of saturated-layered porous media under harmonic waves is evaluated through a semi-analytical solution. The coupled differential equations governing the dynamics of saturated or nearly saturated porous media such as soils containing all the inertial terms of solid and fluid phases are presented for a multi-layer system. Possible simplifications of the equations which are called formulations are introduced based upon the presence of inertial terms associated with the phases. The semi-analytical solutions to the response of multiple layers for all the formulations are presented in terms of pore water pressure and stress variations considering a set of non-dimensional parameters and their respective ratios. Validity of the formulations is presented in a non-dimensional parametric space. The maximum discrepancies in the pore pressure response of the formulations leading to validity regions are illustrated for typical dynamic problems. Subsequently, the effects of layering and drainage conditions on these regions are also presented. The proposed semi-analytical solution may be served as a benchmark one for validating the coupled numerical solutions, which can be used to deal with real scientific and geo-engineering problems in the emerging field of computational geomechanics.     

Lebedev网格改进差分系数TTI介质正演模拟方法研究

本文采用一种新的交错网格-Lebedev网格（LG）进行TTI介质的正演模拟研究,避免了Virieux标准交错网格（SSG）算法在处理TTI、单斜等各向异性介质时波场插值引入的数值误差,提高了模拟精度.在方法实现过程中,本文针对有限差分正演模拟面临的网格频散与边界反射两个关键性问题分别做了优化,并通过模型试算验证了它们的有效性与可行性：（1）结合最小二乘思想推导出新的频散改进差分系数（DIC）,该系数比Taylor系数更能有效地压制粗网格引起的数值频散,可以节约内存,提高计算效率;（2）将分裂的多轴完全匹配层（M-PML）吸收边界条件引入到LG算法中,解决了传统PML边界条件在某些各向异性介质中的不稳定现象并且具有较好的边界吸收效果.     

Accuracy of kinematic wave and diffusion wave approximations for space-independent flows with lateral inflow neglected in the momentum equation

Error equations for the kinematic wave and diffusion wave approximations with lateral inflow neglected in the momentum equation are derived under simplified conditions for space-independent flows. These equations specify error as a function of time in the flow hydrograph. The kinematic wave, diffusion wave and dynamic wave solutions are parameterized through a dimensionless parameter γ which is dependent on the initial conditions. This parameter reflects the effect of initial flow depth, channel-bed slope, lateral inflow and channel roughness when the initial condition is non-vanishing; and it reflects the effect of bed slope, channel roughness and acceleration due to gravity when the initial condition is vanishing. The error equations are found to be the Riccati equation. The structure of the error equations in the case when the momentum equation neglects lateral inflow is different from that when the lateral inflow is included.     

Stochastic modeling of solute transport in 3-D heterogeneous porous media with random source condition

During probabilistic analysis of flow and transport in porous media, the uncertainty due to spatial heterogeneity of governing parameters are often taken into account. The randomness in the source conditions also play a major role on the stochastic behavior in distribution of the dependent variable. The present paper is focused on studying the effect of both uncertainty in the governing system parameters as well as the input source conditions. Under such circumstances, a method is proposed which combines with stochastic finite element method (SFEM) and is illustrated for probabilistic analysis of concentration distribution in a 3-D heterogeneous porous media under the influence of random source condition. In the first step SFEM used for probabilistic solution due to spatial heterogeneity of governing parameters for a unit source pulse. Further, the results from the unit source pulse case have been used for the analysis of multiple pulse case using the numerical convolution when the source condition is a random process. The source condition is modeled as a discrete release of random amount of masses at fixed intervals of time. The mean and standard deviation of concentration is compared for the deterministic and the stochastic system scenarios as well as for different values of system parameters. The effect of uncertainty of source condition is also demonstrated in terms of mean and standard deviation of concentration at various locations in the domain.     

含方形凹陷半无限非均匀介质波动问题FDM模拟

采用规则网格有限差分方法对二维平面弹性波动方程进行差分离散,得到相应的弹性波动方程的有限差分方程,再将弹性波动方程的差分格式与吸收边界、自由边界的离散形式结合形成弹性波动方程有限差分方程解决问题的主体,将其应用于含方形凹陷半无限非均匀介质的模型中进行数值模拟,得到此离散化模型中不同时刻不同节点的位移值。针对具体算例,运用上述方法结合科学计算软件MATLAB和结果后处理软件DIFEM ISOLINE PLOTER得到不同时刻的水平方向位移等值线图与接收器测量点处的合成位移记录,讨论非均匀介质、吸收边界、方形凹陷等对波动特性的影响。     

Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures

Contrast in capillary pressure of heterogeneous permeable media can have a significant effect on the flow path in two-phase immiscible flow. Very little work has appeared on the subject of capillary heterogeneity despite the fact that in certain cases it may be as important as permeability heterogeneity. The discontinuity in saturation as a result of capillary continuity, and in some cases capillary discontinuity may arise from contrast in capillary pressure functions in heterogeneous permeable media leading to complications in numerical modeling. There are also other challenges for accurate numerical modeling due to distorted unstructured grids because of the grid orientation and numerical dispersion effects. Limited attempts have been made in the literature to assess the accuracy of fluid flow modeling in heterogeneous permeable media with capillarity heterogeneity. The basic mixed finite element (MFE) framework is a superior method for accurate flux calculation in heterogeneous media in comparison to the conventional finite difference and finite volume approaches. However, a deficiency in the MFE from the direct use of fractional flow formulation has been recognized lately in application to flow in permeable media with capillary heterogeneity. In this work, we propose a new consistent formulation in 3D in which the total velocity is expressed in terms of the wetting-phase potential gradient and the capillary potential gradient. In our formulation, the coefficient of the wetting potential gradient is in terms of the total mobility which is smoother than the wetting mobility. We combine the MFE and discontinuous Galerkin (DG) methods to solve the pressure equation and the saturation equation, respectively. Our numerical model is verified with 1D analytical solutions in homogeneous and heterogeneous media. We also present 2D examples to demonstrate the significance of capillary heterogeneity in flow, and a 3D example to demonstrate the negligible effect of distorted meshes on the numerical solution in our proposed algorithm.     

Geostatistical modeling of clay spatial distribution in siliciclastic rock samples using the plurigaussian simulation method

In order to implement secondary and enhanced oil recovery processes in complex terrigenous formations as is usual in turbidite deposits, a precise knowledge of the spatial distribution of shale grains is a crucial element for the fluid flow prediction. The reason of this is that the interaction of water with shale grains can significantly modify their size and/or shape, which in turn would cause porous space sealing with the subsequent impact in the flow. In this work, a methodology for stochastic simulations of spatial grains distributions obtained from scanning electron microscopy images of siliciclastic rock samples is proposed. The aim of the methodology is to obtain stochastic models would let us investigate the shale grain behavior under various physico-chemical interactions and flux regimes, which in turn, will help us get effective petrophysical properties (porosity and permeability) at core scale. For stochastic spatial grains simulations a plurigaussian method is applied, which is based on the truncation of several standard Gaussian random functions. This approach is very flexible, since it allows to simultaneously manage the proportions of each grain category in a very general manner and to rigorously handle their spatial dependency relationships in the case of two or more grain categories. The obtained results show that the stochastically simulated porous media using the plurigaussian method adequately reproduces the proportions, basic statistics and sizes of the pore structures present in the studied reference images.     

Analytical and numerical analysis of pressure drawdown in a poroelastic reservoir with complete overburden effect considered

Starting from an analytical reservoir model that incorporates full interaction with an elastic overburden, a new hybrid mathematical approach is developed by combining two numerical discretization methods. A tabular reservoir (petroleum reservoir or an aquifer) in an infinite or semi-infinite domain is viewed as a macroscopic displacement discontinuity, allowing use of the efficient displacement discontinuity mathematical method to calculate stresses and displacements that arise because of pressure changes. A 3-D finite element method using a poroelastic formulation is used to discretize the reservoir itself. By coupling the displacement discontinuity and finite element methods, a 3-D large-scale poroelastic reservoir can be simulated within an infinite or semi-infinite domain. The numerical model has been verified through comparison to known solutions, and some time-dependent pressure drawdown problems are analyzed. Results indicate that including the complete overburden (reservoir surroundings) response has a significant effect on pressure drawdown in a poroelastic reservoir during pumping, and should be incorporated in appropriate applications such as well test equations and subsidence analyses.     

基于第三类边界条件的大地电磁三维矢量有限元任意各向异性正演

地球内部介质的各向异性对地球物理场解译有很大影响, 研究各向异性介质中大地电磁响应具有重要的意义.边界条件是影响电磁场正演精度的一个关键因素, 其中第一类边界条件需要将底面边界设置在离异常体足够远的地方, 面临着计算规模大、求解速度慢的问题.相比第一类边界条件, 第三类边界条件具有计算规模更小、计算精度更高的优点, 在三维各向同性正演中被广泛使用.然而, 各向异性使得第三类边界条件理论变得更为复杂, 目前尚未见到采用第三类边界条件的大地电磁三维各向异性正演.本文推导出各向异性介质中三维矢量有限元底界面的第三类边界条件, 通过模型计算对算法的可靠性和精度进行了验证.结果表明, 当精度相同时, 第三类边界条件的计算规模可以更小; 而当底界面离异常体较近时, 第三类边界条件有着更高的精度.

     

Simulation of normal distributed smooth fields by Karhunen-Loéve expansion in combination with kriging

Simulation of multigaussian stochastic fields can be made after a Karhunen-Loéve expansion of a given covariance function. This method is also called simulation by Empirical Orthogonal Functions. The simulations are made by drawing stochastic coefficients from a random generator. These numbers are multiplied with eigenfunctions and eigenvalues derived from the predefined covariance model. The number of eigenfunctions necessary to reproduce the stochastic process within a predefined variance error, turns out to be a cardinal question. Some ordinary analytical covariance functions are used to evaluate how quickly the series of eigenfunctions can be truncated. This analysis demonstrates extremely quick convergence to 99.5% of total variance for the 2nd order exponential (‘gaussian’) covariance function, while the opposite is true for the 1st order exponential covariance function. Due to these convergence characteristics, the Karhunen-Loéve method is most suitable for simulating smooth fields with ‘gaussian’ shaped covariance functions. Practical applications of Karhunen-Loéve simulations can be improved by spatial interpolation of the eigenfunctions. In this paper, we suggest interpolation by kriging and limits for reproduction of the predefined covariance functions are evaluated.     

Simulation of normal distributed smooth fields by Karhunen-Loéve expansion in combination with kriging

Simulation of multigaussian stochastic fields can be made after a Karhunen-Loéve expansion of a given covariance function. This method is also called simulation by Empirical Orthogonal Functions. The simulations are made by drawing stochastic coefficients from a random generator. These numbers are multiplied with eigenfunctions and eigenvalues derived from the predefined covariance model. The number of eigenfunctions necessary to reproduce the stochastic process within a predefined variance error, turns out to be a cardinal question. Some ordinary analytical covariance functions are used to evaluate how quickly the series of eigenfunctions can be truncated. This analysis demonstrates extremely quick convergence to 99.5% of total variance for the 2nd order exponential (‘gaussian’) covariance function, while the opposite is true for the 1st order exponential covariance function. Due to these convergence characteristics, the Karhunen-Loéve method is most suitable for simulating smooth fields with ‘gaussian’ shaped covariance functions. Practical applications of Karhunen-Loéve simulations can be improved by spatial interpolation of the eigenfunctions. In this paper, we suggest interpolation by kriging and limits for reproduction of the predefined covariance functions are evaluated.     

最优系数有限单元法色散介质GPR频率域正演模拟

工程实际勘探对象如土壤、岩石等多为色散介质,雷达波在其中传播时易发生衰减与畸变,应用常规有限单元法（Finite Element Method,FEM）方法进行数值模拟时,存在数值频散现象.为此,作者以色散介质为研究对象,开展最优系数有限单元法探地雷达（Ground Penetrating Radar,GPR）频率域正演.首先,分析了有限元质量、刚度矩阵的约束条件对有限元求解精度的影响,基于归一化相速度与1的误差最小策略,利用最小二乘法,仅需三个优化参数求取最优的有限元刚度矩阵与质量矩阵.四种不同方法的频散曲线分析及精度对比实验结果表明,优化矩阵在单位波长仅需4.8个网格点下便可达到误差小于0.2%的精度;而一致、集中和折衷矩阵不仅需要更多的网格点,且误差较大.然后,将精确完全匹配层（Exact Perfectly Matched Layer,EPML）吸收边界条件引入最优系数频域有限单元（Finite Element Frequency Domain,FEFD）算法中,简化了吸收参数优化过程,取5层即可达到常规完全匹配层（Perfectly Matched Layer,PML）的10层的吸收效果,能够有效提升正演效率.并将基于EPML的最优系数有限单元法算法引入到城市道路病害模型正演中,实验表明：本文算法能有效压制频散并实现实际色散介质高精度模拟,模拟结果更接近波在地下介质中的实际传播特性.

     

Influence of boundary conditions on the out‐of‐plane response of brick masonry walls in buildings with RC slabs

In modern unreinforced masonry buildings with stiff RC slabs, walls of the top floor are most susceptible to out‐of‐plane failure. The out‐of‐plane response depends not only on the acceleration demand and wall geometry but also on the static and kinematic boundary conditions of the walls. This paper discusses the influence of these boundary conditions on the out‐of‐plane response through evaluation of shake table test results and numerical modelling. As a novum, it shows that the in‐plane response of flanking elements, which are orthogonal to the wall whose out‐of‐plane response is studied, has a significant influence on the vertical restraint at the top of the walls. The most critical configuration exists if the flanking elements are unreinforced masonry walls that rock. In this case, the floor slabs can uplift, and the out‐of‐plane load‐bearing walls loose the vertical restraint at the top. Numerical modelling confirms this experimentally observed behaviour and shows that slab uplift and the difference in base and top excitation have a strong influence on the out‐of‐plane response of the walls analysed. Copyright © 2016 John Wiley & Sons, Ltd.