A continuous/discontinuous Galerkin framework for modeling coupled subsurface and surface water flow

We consider conjunctive surface-subsurface flow modeling, where surface water flow is described by the shallow water equations and ground water flow by Richards’ equation for the vadose zone. Coupling between the models is based on the continuity of flux and water pressure. Numerical approximation of the coupled model using the framework of discontinuous Galerkin (DG) methods is formulated. In the subsurface, the local discontinuous Galerkin (LDG) method is used to approximate ground water velocity and hydraulic head; a DG method is also used to approximate surface water velocity and elevation. This approach allows for a weak coupling of the models and the use of different approximating spaces and/or meshes within each regime. A simplified LDG method based on continuous approximations to water head is also described. Numerical results that investigate physical and numerical aspects of surface–subsurface flow modeling are presented. This work was supported by National Science Foundation grant DMS-0411413.     

双重介质流线模拟技术在天然裂缝性油藏开发中的应用

为了更好地把握天然裂缝性油藏的开发动态,提升大型复杂油藏模型的运算速度,建立了以流线模拟技术为基础的双重介质模型,对天然裂缝性油藏的井型、注水速率、注水周期等影响采收率的因素进行研究,确定了油藏的最佳井型、最佳注水速率、最佳注水周期等指导油藏开发的关键参数,证实了流线模拟技术可以更直观地显示三维驱替前缘。在流线模拟技术与有限差分方法运算时间的比较中,也证实流线模拟技术具有更快的运算速度。     

粘弹性聚合物驱普通稠油微观渗流数学模型

我国海上稠油资源比较丰富,但由于受到海上条件等因素限制,聚合物驱成为提高海上稠油采收率的主要方法.因此深化聚合物溶液驱稠油微观渗流机理对于进一步提高采收率具有十分重要的意义.目前关于粘弹性聚合物渗流机理的理论研究主要局限于弹性聚合物溶液的单相流体在微观孔道内流动特征研究,而针对粘弹性聚合物、油两相流体渗流机理的研究甚少,特别是针对稠油聚合物驱的相关研究未见报道.为此,借助于计算方法较为成熟的OpenFOAM开源平台开展了聚合物驱稠油两相流体渗流机理的研究;以收缩孔道为微观物理模型,建立了粘弹性聚合物溶液、普通稠油两相渗流连续性方程、运动方程及本构方程,并采用VOF（volume of fluid）界面追踪方法建立两相界面相方程;以OpenFOAM开源平台为基础,开发了粘弹性流体、幂律流体两相流体求解器;绘制了不同弹性聚合物溶液在微观孔道内驱油的饱和度分布、速度分布及应力分布特征.结果表明,相对于水驱,纯粘性聚合物溶液前缘突破时间慢,波及面积大,驱油效率高.相比于同等粘度的纯粘性聚合物溶液,粘弹性聚合物的弹性有助于挖潜凸角内的残余油,聚合物溶液的弹性越大,稠油驱油效率越高.随着聚合物溶液弹性的增强,第一法向应力增大,当聚合物溶液进入到孔道突变处时,其弹性发挥的作用最大,法向应力的值最大.研究结果可为矿场实施聚合物驱设计、筛选聚合物溶液提供重要的理论支持.      

A segregated flow scheme to control numerical dispersion for multi-component flow simulations

One of the challenges for reservoir simulation is numerical dispersion. For waterflooding applications the effect is controlled due to the self-sharpening nature of a Buckley–Leverett shock. However, for multi-component flow simulations, incorrect wavespeeds can develop leading to the excessive smearing of fronts because of the coupling of compositional dispersion with the fractional flow. Rather than implementing a higher-order discretization method, we propose a simple scheme based on segregation-in-flow within a gridblock to control numerical dispersion. We extend the method originally proposed for polymer flooding to augmented waterflooding simulations in general as well as simulations of miscible or near miscible gas injection. For compositional simulations of gas injection, this is done through a coupled limited-flash/upstream-exclusion assumption. To test the scheme, an in-house streamline simulator has been modified and validated for modeling low-salinity floods as well as ternary two-phase displacements. Simulation results presented with and without segregation demonstrate the potential of the approach as a heuristic method to control numerical dispersion in multi-component flow simulations.     

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.     

Computational issues of hybrid and multipoint mixed methods for groundwater flow in anisotropic media

In this work, lowest-order Raviart–Thomas and Brezzi–Douglas–Marini mixed methods are considered for groundwater flow simulations. Typically, mixed methods lead to a saddle-point problem, which is expensive to solve. Two approaches are numerically compared here to allow an explicit velocity elimination: (1) the well-known hybrid formulation leading to a symmetric positive definite system where the only unknowns are the Lagrange multipliers and (2) a more recent approach, inspired from the multipoint flux approximation method, reducing low-order mixed methods to cell-centered finite difference schemes. Selected groundwater flow scenarios are used for the comparison between hybrid and multipoint approaches. The simulations are performed in the bidimensional case with a general triangular discretization because of its practical interest for hydrogeologists.     

A numerical approach to study the properties of solutions of the diffusive wave approximation of the shallow water equations

In this paper, we study the properties of approximate solutions to a doubly nonlinear and degenerate diffusion equation, known in the literature as the diffusive wave approximation of the shallow water equations (DSW), using a numerical approach based on the Galerkin finite element method. This equation arises in shallow water flow models when special assumptions are used to simplify the shallow water equations and contains as particular cases the porous medium equation and the p-Laplacian. Diverse numerical schemes have been implemented to approximately solve the DSW equation and have been successfully applied as suitable models to simulate overland flow and water flow in vegetated areas such as wetlands; yet, no formal mathematical analysis has been carried out in order to study the properties of approximate solutions. In this study, we propose a numerical approach as a means to understand some properties of solutions to the DSW equation and, thus, to provide conditions for which the use of the DSW equation may be inappropriate from both the physical and the mathematical points of view, within the context of shallow water modeling. For analysis purposes, we propose a numerical method based on the Galerkin method and we obtain a priori error estimates between the approximate solutions and weak solutions to the DSW equation under physically consistent assumptions. We also present some numerical experiments that provide relevant information about the accuracy of the proposed numerical method to solve the DSW equation and the applicability of the DSW equation as a model to simulate observed quantities in an experimental setting.     

Modeling semi-steady state near-well flow performance for horizontal wells in anisotropic reservoirs

This paper presents a novel methodology to model semi-steady state horizontal well flow performance in an anisotropic reservoir taking into account flow in the near-well region for an arbitrary well trajectory. It is based on an analytical productivity model describing coupled axial reservoir flow and radial well inflow. In order to apply this model in an anisotropic reservoir, the permeability field relative to the radial direction perpendicular to the well trajectory and the axial direction along the well trajectory must first be determined. A classical space transformation is used in concert with rotational transforms to obtain a virtual isotropic model. The transformation preserves the volumes and pressures. It is not a novel concept, but different from previous approaches in the sense that it is only applied in the near-well domain to formulate an equally isotropic media. As a result, the use of this virtual isotropic model requires the Dietz shape factor for an ellipse, transformed from the original cylindrical near-well domain. The Dietz shape factors are determined numerically in this research. The semi-steady state well/near-well model is implemented in a numerical simulator incorporating formation anisotropy and wellbore hydraulics. The specific productivity index along the well trajectory is generated using the virtual configuration. Numerical results for different anisotropy ratios and also incorporating frictional losses in the well are presented. Furthermore, the well/near-well model is applied in coupling with streamline reservoir model for a water flooding case. This appears to be the first coupling of a well hydraulics model and a streamline simulator. It presents the application of the well/near-well model in integrated reservoir simulation in an efficient and accurate manner. The results demonstrate that the coupling approach with a streamline reservoir model and the well/near-well is of great potential for advanced well simulation efficiently.     

基坑潜水疏干的简化计算在昊华水泥厂中的应用

以流线、流面、汇点的概念为基础,对稳定流双井干扰和直线隔水边界附近涌水量理论公式进行对比分析,提出了二个虚拟界面,其中虚拟界面Ⅰ,运用流线、流面的性质,流线方程等给出证明;虚拟界面Ⅱ则通过半无限条形降落漏斗的分析,应用元流和总流的能量方程得到流量为零,流线为零的平面。在同样条件下,条形无限涌水量是半无限潜含水层涌水量的二倍。应用总流能量方程对三种情况水头损失的分析,解释了这种关系存在的合理性,得出虚拟界面Ⅱ,并以此得出该界面内的最大残余水头计算公式。将基坑降水运用虚拟界面简化为扇形,条形半无限含水层,从而实现单井预测,该方法应用到昊华水泥厂基坑降水中,预测效果理想。     

Accuracy and efficiency of time integration methods for 1D diffusive wave equation

The diffusive wave approximation of the Saint-Venant equations is commonly used in hydrological models to describe surface flow processes. Numerous numerical approaches can be used to solve this highly nonlinear equation. Nonlinear time integration schemes—also called methods of lines (MOL)—were proven very efficient to solve other nonlinear problems in geosciences but were never considered to deal with surface flow modeling with the diffusive wave equation. In this paper, we study the relative performance of different time and space integration schemes by comparing the results obtained with classical approaches and with nonlinear time integration approaches. The results show that (i) the integration method with a higher order in space shows high accuracy regarding an integrated indicator such as the global mass balance error but is less accurate regarding local indicators, and (ii) nonlinear time integration techniques perform better than classical ones. Overall, it seems that integration techniques combining nonlinear time integration and a low spatial order need to be considered when developing hydrological modeling tools owing to their simplicity of implementation and very good performance.     

Isogeometric analysis for unsaturated flow problems

Unsaturated flow problems in porous media often described by Richards’ equation are of great importance in many engineering applications. In this contribution, we propose a new numerical flow approach based on isogeometric analysis (IGA) for modeling the unsaturated flow problems. The non-uniform rational B-spline (NURBS) basis is utilized for spatial discretization whereas the stable implicit backward Euler method for time discretization. The nonlinear Richards’ equation is iteratively solved with the aid of the Newton–Raphson scheme. Owing to some desirable features of an efficient numerical flow approach, major advantages of the present formulation involve: (a) numerical oscillation at the wetting front can be avoided or facilitated, simply by using either an h-refinement or a lumped mass matrix technique; (b) higher-order exactness can be obtained due to the nature of the IGA features; (c) the approach is straightforward to implement and it does not need any transformation, e.g., Kirchhoff transformation or filter algorithm; and (d) in contrast to the Picard iteration scheme, which forms linear convergences, the proposed approach can however yield quadratic convergences by using the Newton–Raphson method for solving resultant nonlinear equations. Numerical model validation is analyzed by solving a three-dimensional unsaturated flow problem in soil, and its derived results are verified against analytical solutions. Numerical applications are then studied by considering three extensive examples with simple and complex configurations to further show the accuracy and applicability of the present IGA.     

Flood hazard assessment with multiparameter approach derived from coupled 1D and 2D hydrodynamic flow model

Hydrodynamic flow modeling is carried out using a coupled 1D and 2D hydrodynamic flow model in northern India where an industrial plant is proposed. Two flooding scenarios, one considering the flooding source at regional/catchment level and another considering all flooding sources at local level have been simulated. For simulating flooding scenario due to flooding of the upstream catchment, the probable maximum flood (PMF) in the main river is routed and its flooding impact at the plant site is studied, while at the local level flooding, in addition to PMF in the main river, the probable maximum precipitation at the plant site and breaches in the canals near the plant site have been considered. The flood extent, depth, level, duration and maximum flow velocity have been computed. Three parameters namely the flood depth, cross product of flood depth and velocity and flood duration have been used for assessing the flood hazard, and a flood hazard classification scheme has been proposed. Flood hazard assessment for flooding due to upstream catchment and study on local scale facilitates determination of plinth level for the plant site and helps in identifying the flood protection measures.     

A mixed finite element method for Hele-Shaw cell equations

A new numerical method to solve the system of equations describing two phase flow in a Hele-Shaw cell is presented. It combines a mixed finite element method, the method of subtraction of the singularity and a front tracking grid in a single computational strategy. This choice of discretization techniques is well motivated by the difficulties present in the system of equations and the physics of the problem. The new method was tested against analytical solutions and also by solving the Saffman–Taylor viscous fingering problem for finite and infinite mobility ratios. In both cases convergence under mesh refinement is achieved for the fingers developed from an initial sinusoidal interface. Finger splitting is observed for low values of the surface tension and high mobility ratio. Different explanations, based in our results, are provided for this phenomenon. This revised version was published online in July 2006 with corrections to the Cover Date.     

Fast computation of arrival times in heterogeneous media

We continue the work that was initiated in (K. H. Karlsen, K.-A. Lie, and N. H. Risebro. A fast marching method for reservoir simulation. Comp. Geo., 4(2) (2000)185–206) on a marching method for simulating two-phase incompressible immiscible flow of water and oil in a porous medium. We first present an alternative derivation of the marching method that reveals a strong connection to modern streamline methods. Then, through the study of three numerical test cases we present two deficiencies: (i) the original marching algorithm does not always compute the correct solution of the underlying difference equations, and (ii) the method gives largely inaccurate arrival times in the presence of large jumps within the upwind difference stencil. As a remedy of the first problem, we present a new advancing-front method, which is faster than the original marching method and guarantees a correct solution of the underlying discrete linear system. To cure the second problem, we present two adaptive strategies that avoid the use of finite-difference stencils containing large jumps in the arrival times. The original marching method was introduced as a fast tool for simulating two-phase flow scenarios in heterogeneous formations. The new advancing-front method has limited applicability in this respect, but may rather be used as a fast and relatively accurate method for computing arrival times and derived quantities in heterogeneous media.     

高堆石坝中期度汛挡水风险率估计

在考虑高堆石坝坝体中期临时挡水度汛不允许坝面过水条件下,综合考虑坝前洪水位与防洪度汛高程的随机性,构建高堆石坝工程中期度汛挡水风险数学模型.针对堆石坝施工系统的特点,考虑各月停工天数与日平均上升速度的随机性,建立了防洪度汛高程仿真计算模型.基于Monte-Carlo方法原理,耦合水文、水力和施工随机因素对风险模型进行求解,并研究了日均控制最低上升速度这一施工可控指标对风险率的影响.通过大渡河流域某大型工程实例分析表明,风险模型与计算方法可行、有效.     

A fully-coupled discontinuous Galerkin method for two-phase flow in porous media with discontinuous capillary pressure

In this paper, we formulate and test numerically a fully-coupled discontinuous Galerkin (DG) method for incompressible two-phase flow with discontinuous capillary pressure. The spatial discretization uses the symmetric interior penalty DG formulation with weighted averages and is based on a wetting-phase potential/capillary potential formulation of the two-phase flow system. After discretizing in time with diagonally implicit Runge-Kutta schemes, the resulting systems of nonlinear algebraic equations are solved with Newton’s method and the arising systems of linear equations are solved efficiently and in parallel with an algebraic multigrid method. The new scheme is investigated for various test problems from the literature and is also compared to a cell-centered finite volume scheme in terms of accuracy and time to solution. We find that the method is accurate, robust, and efficient. In particular, no postprocessing of the DG velocity field is necessary in contrast to results reported by several authors for decoupled schemes. Moreover, the solver scales well in parallel and three-dimensional problems with up to nearly 100 million degrees of freedom per time step have been computed on 1,000 processors.     

地下水流场三维流线可视化模拟与实现

流线对流场的描述是场描述的一种基本形式,在生成流线算法中质点追踪算法是简单而快速的。但是,质点追踪算法在流线行进过程中会产生歧义点。为了解决地下水流场三维流线可视化中质点追踪算法产生歧义点的问题,本文从流线的定义出发,以流线方向为K·J取最大值的方向为原则,提出了降维行进算法,并对一个理论模型定义的流场进行了流线可视化研究,实现了流场的三维流线模拟与可视化表述,其结果与理论研究相吻合。     

非结构网格上的三维浅水流动数值模型

针对当前复杂环境水流模拟的需求,建立了新型的基于特征型高分辨率数值算法的三维非结构网格浅水动力模型。模型采用有限体积法离散sigma坐标下的三维浅水方程,运用Roe黎曼近似解评估水平界面通量。模型网格拟合边界能力强,可根据需要局部加密;格式数值性能优良,具有守恒性、单调迎风性、高数值分辨率等特性。同时,应用干湿判别法处理动边界,以适应浅滩地形漫/露过程模拟的需要。封闭水池内部风生环流、干河床上溃坝过程和长江口实际潮流场的模拟从不同侧面展示了模型的特点,结果表明它能够准确地预测水流的三维流动结构,而且计算简单高效,具有良好的数值稳定性。     

A fast marching method for reservoir simulation

We present a fast marching level set method for reservoir simulation based on a fractional flow formulation of two-phase, incompressible, immiscible flow in two or three space dimensions. The method uses a fast marching approach and is therefore considerably faster than conventional finite difference methods. The fast marching approach compares favorably with a front tracking method as regards both efficiency and accuracy. In addition, it maintains the advantage of being able to handle changing topologies of the front structure. This revised version was published online in July 2006 with corrections to the Cover Date.     

Fast 3D Reservoir Simulation and Scale Up Using Streamtubes

This paper presents an implementation of a semianalytical method for oil recovery calculation in heterogeneous reservoirs that is both fast and accurate. The method defines streamline paths based on a conventional single-phase incompressible flow calculation. By calculating the time-of-flight for a particle along a streamline and assigning a volumetric flux to each streamline, the cumulative pore volume of a streamtube containing the streamline can be calculated. Subsequently, the streamtube geometries are kept constant and the effects of the time varying mobility distribution in two-phase flow are accounted for by varying the flow rate in each streamtube, based on fluid resistance changes along the streamtube. Oil recovery calculations are then done based on the 1D analytical Buckley–Leverett solution. This concept makes the method extremely fast and easy to implement, making it ideal to simulate large reservoirs generated by geostatiscal methods. The simulation results of a 3D heterogeneous reservoir are presented and compared with those of other simulators. The results shows that the new simulator is much faster than a traditional finite difference simulator, while having the same accuracy. The method also naturally handles the upscaling of absolute and relative permeability. We make use of these upscaling abilities to generate a coarse curvilinear grid that can be used in conventional simulators with a great advantage over conventional upscaled Cartesian grids. This paper also shows an upscaling example using this technique.