我国北方型裂隙岩溶水流及污染物运移数值模拟研究进展

简要回顾了近几十年来我国北方型裂隙岩溶水流及污染物运移数值模拟方面的研究成果,评述了裂隙岩溶水资源评价和污染物运移模拟的计算方法和模拟技术,其中对一些与数值模拟相关的重要问题、近年来研究的热点问题,如模型参数识别及其尺度效应问题、数值弥散和数值振荡问题、Non-Darcy流和Non-Fickian运移问题以及随机方法研究等进行了重点评述,对已有的解决办法进行了比较和对照,并探讨了今后的研究趋势.     

A mixed solution for two-dimensional unsteady flow in fractured porous media

A mixed finite element–boundary element solution for the analysis of two-dimensional flow in porous media composed of rock blocks and discrete fractures is described. The rock blocks are modelled implicitly by using boundary elements whereas finite elements are adopted to model the discrete fractures. The computational procedure has been implemented in a hybrid code which has been validated first by comparing the numerical results with the closed-form solution for flow in a porous aquifer intercepted by a vertical fracture only. Then, a more complex problem has been solved where a pervious, homogeneous and isotropic matrix containing a net of fractures is considered. The results obtained are shown to describe satisfactorily the main features of the flow problem under study. © 1997 by John Wiley & Sons, Ltd.     

A dual-porosity model considering the displacement effect for incompressible two-phase flow

The dual-porosity model is usually employed to simulate the flow in fractured reservoirs. However, its original form for the multiphase flow does not consider the displacement effect under macropressure gradient. Especially for the incompressible multiphase flow, it predicts zero transfer term between fracture and matrix, which is unreasonable. To improve this, a modified double-porosity model is proposed for incompressible two-phase flow, in which the displacement effect is considered and the corresponding shape factor is derived. For the anisotropic case, the shape factor of displacement depends upon the velocity direction. The accuracy and the efficiency of the proposed dual-porosity model are indicated through numerical tests.     

模拟裂隙多孔介质中变饱和渗流的广义等效连续体方法

描述了一种计算裂隙多孔介质中变饱和渗流的广义等效连续体方法。这种方法忽略裂隙的毛细作用,设定一个与某孔隙饱和度相对应的综合饱和度极限值,并假定：（1）如果裂隙多孔介质的综合饱和度小于该极限值,水只在孔隙中存在并流动,而裂隙中则没有水的流动;（2）如果综合饱和度等于或大于该极限值,水将进入裂隙,并在裂隙内运动。分析比较了等效连续体模型的不同计算方法,并给出了一个模拟裂隙岩体中变饱和渗流与传热耦合问题的应用算例。结果表明,所述方法具有一般性,可以有效地模拟裂隙多孔介质中变饱和渗流的基本特征。     

A TIME-STEPPING FINITE ELEMENT METHOD FOR ANALYSIS OF CONTAMINANT TRANSPORT IN FRACTURED POROUS MEDIA

This paper describes the development of a finite element method for analysing contaminant transport in double-porosity geomaterials using a time-stepping approach. In many cases, double-porosity models may be used to represent fractured rock formations and fissured soils. A distinctive feature of utilizing this kind of model is that it is not necessary to have an intimate knowledge of the nature, distribution and properties of individual fractures and fracture arrangement since the fracture geometry and details are considered only in an averaged or equivalent continuum sense. The flux exchange that occurs between the fluid in the fractures and in the solid matrix is represented by a linear heriditary process. This has the consequence that in order to carry the solution forward from time t to t+Δt, it is necessary to know and to store the complete contaminant history up to time t. This paper shows that all the hereditary information necessary to carry the solution forward is contained in the values of certain hereditary variables at time t so that it is not necessary to store the complete time history and consequently a more efficient numerical process can be developed.     

One‐dimensional contaminant transport through a deforming porous medium: theory and a solution for a quasi‐steady‐state problem

Contaminant migration through soil is usually modelled mathematically using the dispersion–advection equation. This type of model finds application when planning the remediation of contaminated land, predicting the movement of polluted groundwater and designing engineered landfills. Usually the analysis assumes that the porous media through which the contaminant migrates is stationary. However, the construction of landfills on clay soils means that the soil beneath the landfill will undergo time‐dependent deformation as the soil consolidates. To date, there are no published data on the effect a deforming porous media may have on contaminant transport beneath a landfill; indeed, there appears to be no theory of contaminant migration through a deforming soil. In this paper, a one‐dimensional theory of contaminant migration through a saturated deforming porous media is developed based on a small and large strain analysis of a consolidating soil and conservation of contaminant mass. By selection of suitable parameters, the new transport equation reduces to the familiar one‐dimensional dispersion–advection equation for a saturated soil with linear, reversible, equilibrium controlled sorption of the contaminant onto the soil skeleton. Analytic solutions to a quasi‐steady‐state contaminant transport problem for a deforming media are presented, and a preliminary assessment made of the potential importance of soil deformation on the results of a contaminant migration analysis. Copyright © 2000 John Wiley & Sons, Ltd.     

Transverse and longitudinal fluid flow modelling in fractured porous media with non‐matching meshes

A new discrete fracture model is introduced to simulate the steady‐state fluid flow in discontinuous porous media. The formulation uses a multi‐layered approach to capture the effect of both longitudinal and transverse permeability of the discontinuities in the pressure distribution. The formulation allows the independent discretisation of mesh and discontinuities, which do not need to conform. Given that the formulation is developed at the element level, no additional degrees of freedom or special integration procedures are required for coupling the non‐conforming meshes. The proposed model is shown to be reliable regardless of the permeability of the discontinuity being higher or lower than the surrounding domain. Four numerical examples of increasing complexity are solved to demonstrate the efficiency and accuracy of the new technique when compared with results available in the literature. Results show that the proposed method can simulate the fluid pressure distribution in fractured porous media. Furthermore, a sensitivity analysis demonstrated the stability regarding the condition number for wide range values of the coupling parameter.     

A semi‐analytical solution of surface flow and subsurface flow on a slope land under a uniform rainfall excess

A new approach is proposed to analyze the surface flow and subsurface flow passing over a pervious ground under a uniform rainfall excess. The flow field is divided into two regions that are called water layer and soil layer. To figure out the hydraulic behavior of overland flow on an inclined plane under a rainfall event, the simplified Navier–Stokes equations are employed for the surface water flow, and the flow inside the soil layer is porous media flow, which is governed by Biot's (1956, 1962) theory of poroelasticity. The velocity distribution of overland flow is nonzero at the ground surface. The relation between water depth and slope length was developed first. The profile of surface water flow was then found backwards from the downstream end of the flow section by the Runge–Kutta method. After that, the flow velocity and flow discharge of each layer could also be obtained via the water depth. Finally, the variation of fluid shear stress inside the soil layer is also discussed. Copyright © 2011 John Wiley & Sons, Ltd.     

Theoretical and numerical study of the steady‐state flow through finite fractured porous media

This paper investigates the two‐dimensional flow problem through an anisotropic porous medium containing several intersecting curved fractures. First, the governing equations of steady‐state fluid flow in a fractured porous body are summarized. The flow follows Darcy's law in matrix and Poiseuille's law in fractures. An infinite transversal permeability is considered for the fractures. A multi‐region boundary element method is used to derive a general pressure solution as a function of discharge through the fractures and the pressure and the normal flux on the domain boundary. The obtained solution fully accounts for the interaction and the intersection between fractures. A numerical procedure based on collocation method is presented to compute the unknowns on the boundaries and on the fractures. The numerical solution is validated by comparing with finite element solution or the results obtained for an infinite matrix. Pressure fields in the matrix are illustrated for domains containing several interconnected fractures, and mass balance at the intersection points is also checked. Copyright © 2013 John Wiley & Sons, Ltd.     

采用复合单元法模拟裂隙多孔介质变饱和流动

采用复合单元法建立了模拟裂隙多孔介质变饱和流动的数值模型。该模型具有以下特点:裂隙不需要离散成特定单元,而是根据几何位置插入到孔隙基质单元中形成复合单元;在复合单元中,分别建立裂隙流和孔隙基质流的计算方程,二者通过裂隙-基质界面产生联系并整合成复合单元方程;复合单元方程具有和常规有限单元方程相同的格式,因此,可以使用常规有限单元方程的求解技术。采用欠松弛迭代、集中质量矩阵以及自适应时步调节等技术,开发了裂隙多孔介质变饱和流动计算程序。通过模拟一维干土入渗和复杂裂隙含水层内的流动问题,验证了该模型的合理性和适用性。模拟结果为进一步认识非饱和裂隙含水层地下水流动特性提供了理论依据。     

A fracture mapping and extended finite element scheme for coupled deformation and fluid flow in fractured porous media

This paper presents a fracture mapping (FM) approach combined with the extended finite element method (XFEM) to simulate coupled deformation and fluid flow in fractured porous media. Specifically, the method accurately represents the impact of discrete fractures on flow and deformation, although the individual fractures are not part of the finite element mesh. A key feature of FM‐XFEM is its ability to model discontinuities in the domain independently of the computational mesh. The proposed FM approach is a continuum‐based approach that is used to model the flow interaction between the porous matrix and existing fractures via a transfer function. Fracture geometry is defined using the level set method. Therefore, in contrast to the discrete fracture flow model, the fracture representation is not meshed along with the computational domain. Consequently, the method is able to determine the influence of fractures on fluid flow within a fractured domain without the complexity of meshing the fractures within the domain. The XFEM component of the scheme addresses the discontinuous displacement field within elements that are intersected by existing fractures. In XFEM, enrichment functions are added to the standard finite element approximation to adequately resolve discontinuous fields within the simulation domain. Numerical tests illustrate the ability of the method to adequately describe the displacement and fluid pressure fields within a fractured domain at significantly less computational expense than explicitly resolving the fracture within the finite element mesh. Copyright © 2013 John Wiley & Sons, Ltd.     

多孔介质油水两相k~s~p关系数学模型的实验研究

数学模型是研究相对渗透率与饱和度关系曲线的重要方法。采用自行开发设计的人工平面多孔介质模型,测定了相对渗透率与饱和度的关系曲线。多孔介质选择粒径为0.5~1mm、1~2mm的标准砂,纯净的水为湿润相,用3号苏丹红染色的93#汽油为非湿润相,组成多孔介质油水两相流动系统。采用Van Genuchten and Mualeum(VGM)和Brooks-Corey-Burdine(BCB)两种数学模型计算相对渗透率与饱和度的关系曲线,通过比较两种数学模型计算结果之间和模型计算结果与实测结果的差异以及模型的应用、多相渗流系统自身特征,得出VGM、BCB两种数学模型计算结果符合实际情况,VGM模型应用过程更为简便,但VGM模型具有一定适用条件;在砂性多孔介质中,BCB模型计算相对渗透率与饱和度关系曲线更准确。     

JHomogenizer: a computational tool for upscaling permeability for flow in heterogeneous porous media

This paper presents an object-oriented programming approach for the design of numerical homogenization programs, called JHomogenizer. It currently includes five functional modules to compute effective permeability and simple codes for computing solutions for flow in porous media. Examples with graphical output are shown to illustrate some functionalities of the program. A series of numerical examples demonstrates the effectiveness of the methodology for two-phase flow in heterogeneous reservoirs. The software is freely available, and the open architecture of the program facilitates further development and can adapt to suit specific needs easily and quickly.     

Modelling of steady‐state fluid flow in 3D fractured isotropic porous media: application to effective permeability calculation

In this paper, 3D steady‐state fluid flow in a porous medium with a large number of intersecting fractures is derived numerically by using collocation method. Fluid flow in the matrix and fractures is described by Darcy's law and Poiseuille's law, respectively. The recent theoretical development presented a general potential solution to model the steady‐state flow in fractured porous media under a far‐field condition. This solution is a hypersingular integral equation with pressure field in the fracture surfaces as the main unknown. The numerical procedure can resolve the problem for any form of fractures and also takes into account the interactions and the intersection between fractures. Once the pressure field and then the flux field in fractures have been determined, the pressure field in the porous matrix is computed completely. The basic problem of a single fracture is investigated, and a semi‐analytical solution is presented. Using the solution obtained for a single fracture, Mori‐Tanaka and self‐consistent schemes are employed for upscaling the effective permeability of 3D fractured porous media. Copyright © 2012 John Wiley & Sons, Ltd.     

On zero-thickness interface elements for diffusion problems

The present study focuses on the hydraulic behaviour of joints, and, specially, on its numerical implementation in terms of the FEM analysis using a discrete fracture flow approach. Fluid flow through discontinuities has traditionally been modelled using special elements of zero-thickness, which we can classify into single, double and triple-nodded. Single node elements are the simplest and consist of ‘line’ or ‘pipe’ elements which are superimposed onto the standard continuum mesh and that can only model the longitudinal conductivity through the discontinuity. On the other hand, some authors have included the influence of a transversal conductivity, and the subsequent localized potential drop, by using triple node interface elements. In those, the two nodes of the adjacent continuum elements represent the potentials in the pore system on each side of the interface, and a third node in the middle represents the average potential of the fluid in the channel represented by the discontinuity. Finally, double node interface elements have also been proposed, which have the advantage of making it possible to use the same FE mesh for both mechanical and flow analysis. In some cases the influence of a transversal conductivity is not considered and, therefore, although geometrically double-nodded, these elements belong to the single node type and when time comes to solve the system the two nodes must have the same potential, which can only be obtained by the ‘trick’ of prescribing the equivalence of these two d.o.f. before solving the global system of equations. This limitation may, however, be avoided by assuming that the potential in the channel is the average of the two sides of the interface. Based in this simple assumption, an alternative flow interface model has been recently developed and implemented, which preserves both longitudinal and transversal conductivities. An application example is developed and solved with the three types of interfaces described. The results offer useful information regarding the range of applicability and limitations of the new double-nodded interface element proposed. Copyright © 2004 John Wiley & Sons, Ltd.     

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.     

On the velocity covariance for steady flows in heterogeneous porous formations and its application to contaminants transport

We consider groundwater steady flow in a heterogeneous porous formation of random and stationary log-conductivity Y = ln K, characterized by the mean 〈Y〉, and the two point correlation function C _Y which in turn has finite, and different horizontal and vertical integral scales I and I _v , respectively. The fluid velocity V, driven by a given head drop applied at the boundary, has constant mean value U ≡ (U, 0, 0). Approximate explicit analytical expressions for transverse velocity covariances are derived. The adopted methodology follows the approach developed by Dagan and Cvetkovic (Spatial moments of kinetically sorbing plume in a heterogeneous aquifers, Water Resour. Res. 29 (1993) 4053) to obtain a similar result for the longitudinal velocity covariance. Indeed, the approximate covariances of transverse velocities are determined by requiring that they have the exact first order variances as well as zero integral scale (G. Dagan, Flow and Transport in Porous Formations (Springer, 1989)) , and provide the exact asymptotic limits of the displacement covariance of the fluid particles obtained by Russo (On the velocity covariance and transport modeling in heterogeneous anisotropic porous formations 1. Saturated flow, Water Resour. Res., 31 (1995) 129). Comparisons with numerical results show that the proposed expressions compare quite well in the early time regime, and for Ut/I >100. Since most of the applications, like assessing the effective mobility of contaminants or quantifying the potential hazards of nuclear repositories, require predictions over higher times the proposed approximate expressions provide acceptable results. The main advantage related to such expressions is that they allow obtaining closed analytical forms of spatial moments pertaining to kinetically sorbing contaminant plumes avoiding the very heavy computational effort which is generally demanded. For illustration purposes, we consider the movement of one contaminant species, and show how our approximate spatial moments compare with the numerical simulations.     

分布热源作用下裂隙岩体渗流-传热的拉氏变换-格林函数半解析计算方法

以裂隙岩体高放射性核废物地下处置库性能评估为目标,提出了分布热源作用下单裂隙岩体渗流-传热的简化概念模型、控制微分方程和拉氏变换-格林函数半解析法,为进一步采用半解析法计算分布热源作用下多裂隙岩体的渗流-传热问题奠定了基础。针对单裂隙岩体的渗流-传热问题,建立考虑岩石内热源和二维热传导的控制微分方程,利用拉氏变换域微分方程的基本解建立格林函数积分方程,采用解析法处理其中的奇点,通过数值积分和拉氏数值逆变换求解,计算任意时刻裂隙水和岩石的温度分布。通过算例,与基于岩石一维热传导假定的解析解进行了对比,并计算分析了分布热源作用下单裂隙岩体的渗流-传热特征及其对裂隙开度、岩石热传导系数和热流集度的敏感度。算例表明,（1）就裂隙水温度而言,由于考虑了岩石的二维热传导,拉氏变换-格林函数半解析解小于基于岩石一维热传导假定的解析解;（2）裂隙水温度和岩石温度对裂隙开度和热流集度的敏感度较大,对岩石热传导系数的敏感度较小。     

A pore-scale numerical model for flow through porous media

A pore-scale numerical model based on Smoothed Particle Hydrodynamics (SPH) is described for modelling fluid flow phenomena in porous media. Originally developed for astrophysics applications, SPH is extended to model incompressible flows of low Reynolds number as encountered in groundwater flow systems. In this paper, an overview of SPH is provided and the required modifications for modelling flow through porous media are described, including treatment of viscosity, equation of state, and no-slip boundary conditions. The performance of the model is demonstrated for two-dimensional flow through idealized porous media composed of spatially periodic square and hexagonal arrays of cylinders. The results are in close agreement with solutions obtained using the finite element method and published solutions in the literature. Copyright © 1999 John Wiley & Sons, Ltd.     

A finite element method for modeling coupled flow and deformation in porous fractured media

Modeling the flow in highly fractured porous media by finite element method (FEM) has met two difficulties: mesh generation for fractured domains and a rigorous formulation of the flow problem accounting for fracture/matrix, fracture/fracture, and fracture/boundary fluid mass exchanges. Based on the recent theoretical progress for mass balance conditions in multifractured porous bodies, the governing equations for coupled flow and deformation in these bodies are first established in this paper. A weak formulation for this problem is then established allowing to build a FEM. Taking benefit from recent development of mesh‐generating tools for fractured media, this weak formulation has been implemented in a numerical code and applied to some typical problems of hydromechanical coupling in fractured porous media. It is shown that in this way, the FEM that has proved its efficiency to model hydromechanical phenomena in porous media is extended with all its performances (calculation time, couplings, and nonlinearities) to fractured porous media. Copyright © 2015 John Wiley & Sons, Ltd.