A Godunov method for localization in elastoplastic granular flow

The equations governing the elastic-plastic deformation of granular materials are typically hyperbolic, or contain small-magnitude damping or rate effects. A finite element algorithm is the standard method for the numerical integration of these systems. In particular, finite elements allow great flexibility in the design of grid geometry. However, modern finite difference methods for hyperbolic systems have been successful in aerodynamics computations, resolving wave structures more sharply than finite element schemes. In this paper we develop a finite difference scheme for granular flow problems. We report on a second-order Godunov-type scheme for the integration of hyperbolic equations for the elastoplastic deformation of a simple model of granular flow. The Godunov method includes a characteristic tracing step in the integration, providing minimal wave dispersion, and a slope limiting step, preventing unphysical oscillations. The granular flow model we consider is hyperbolic, but hyperbolicity is lost at a large value of accumulated plastic strain. This loss of hyperbolicity is a tell-tale signal for the formation of a shear band within the sample. Typically, when systems lose hyperbolicity a regularization mechanism is added to the model equations in order to maintain the well posedness of the system. These regularizations include viscosity, viscoplasticity, higher-order gradient effects or stress coupling. Here we appeal to a very different kind of regularization. When the system loses hyperbolicity and a shear band forms, we treat the band as an internal boundary, and impose jump conditions at this boundary. Away from the band, the system remains hyperbolic and the integration step proceeds as usual.     

非均匀层状介质一维波动方程精确解的有限差分算法

平面波的传播问题通常可以归结为一维波动方程的定解问题。在非均匀介质中,即使简单的一维波动方程也需要借助于数值方法获得近似解。3层5点古典差分格式是计算偏微分方程一种常用算法,作为一种显式迭代格式,需要满足稳定性条件 ,其中 为波速, 为空间采样间隔, 为时间采样间隔。当 时, ,古典差分格式达到临界稳定状态。在这种情况下,平面波在 时间内的传播距离恰好等于空间采样间隔,差分格式真实地反映了平面波的传播原理,因而可以得到一维波动方程的精确解。但是,由于在非均匀介质中存在不连续的波阻抗界面,此方法不适于计算非均匀介质的波场。为了将临界稳定情况下的古典差分格式推广应用至非均匀层状介质,提出了一种能够处理波阻抗界面的有限差分格式,并应用傅里叶分析法得到其稳定性条件。模型算例验证了此算法的正确性。     

Implementation of a Locally Conservative Numerical Subgrid Upscaling Scheme for Two-Phase Darcy Flow

We present a locally mass conservative scheme for the approximation of two-phase flow in a porous medium that allows us to obtain detailed fine scale solutions on relatively coarse meshes. The permeability is assumed to be resolvable on a fine numerical grid, but limits on computational power require that computations be performed on a coarse grid. We define a two-scale mixed finite element space and resulting method, and describe in detail the solution algorithm. It involves a coarse scale operator coupled to a subgrid scale operator localized in space to each coarse grid element. An influence function (numerical Greens function) technique allows us to solve these subgrid scale problems independently of the coarse grid approximation. The coarse grid problem is modified to take into account the subgrid scale solution and solved as a large linear system of equations posed over a coarse grid. Finally, the coarse scale solution is corrected on the subgrid scale, providing a fine grid representation of the solution. Numerical examples are presented, which show that near-well behavior and even extremely heterogeneous permeability barriers and streaks are upscaled well by the technique.     

Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media

Fully implicit time-space discretizations applied to the two-phase Darcy flow problem leads to the systems of nonlinear equations, which are traditionally solved by some variant of Newton’s method. The efficiency of the resulting algorithms heavily depends on the choice of the primary unknowns since Newton’s method is not invariant with respect to a nonlinear change of variable. In this regard, the role of capillary pressure/saturation relation is paramount because the choice of primary unknowns is restricted by its shape. We propose an elegant mathematical framework for two-phase flow in heterogeneous porous media resulting in a family of formulations, which apply to general monotone capillary pressure/saturation relations and handle the saturation jumps at rocktype interfaces. The presented approach is applied to the hybrid dimensional model of two-phase water-gas Darcy flow in fractured porous media for which the fractures are modelled as interfaces of co-dimension one. The problem is discretized using an extension of vertex approximate gradient scheme. As for the phase pressure formulation, the discrete model requires only two unknowns by degree of freedom.     

A technique to test finite difference schemes to model some geophysical processes in a geological structure

A preliminary problem to solve in the set-up of a mathematical model simulating a geophysical process is the choice of a suitable discrete scheme to approximate the governing differential equations. This paper presents a simple technique to test finite difference schemes used in the modeling of geophysical processes occurring in a geological structure. This technique consists in generating analytical solutions similar to the ones characterizing a geophysical process, given general information on some relevant parameters. Useful information for the choice of the discrete scheme to employ in the mathematical model simulating the original geophysical process can be obtained from the comparison between the analytical solution and the approximated numerical solutions generated by means of different discrete schemes. Two classes of numerical examples approximating the differential equation that governs the steady state earth's heat flow have been treated using three different finite differences schemes. The first class of examples deals with media whose phenomenological parameters vary as continuous space functions; the second class, instead, deals with media whose phenomenological parameters vary as discontinuous space functions. The finite difference schemes that have been utilized are: Centered Finite Difference Scheme (CDS), Arithmetic Mean Scheme (AMS), and Harmonic Mean Scheme (HMS).The numerical simulations showed that: the CDS may yield physically inconsistent solutions if the lattice internodal distance is too large, but in case of phenomenological parameters varying as a continuous function, this pitfall can be avoided increasing the lattice node refinement. In case of phenomenological parameters varying as a discontinuous function, instead, the CDS may yield physically inconsistent solutions for any lattice-node refinement. The HMS produced good results for both classes of examples showing to be a scheme suitable to model situations like these.     

Multiscale simulations of porous media flows in flow-based coordinate system

In this paper, we propose a multiscale technique for the simulation of porous media flows in a flow-based coordinate system. A flow-based coordinate system allows us to simplify the scale interaction and derive the upscaled equations for purely hyperbolic transport equations. We discuss the applications of the method to two-phase flows in heterogeneous porous media. For two-phase flow simulations, the use of a flow-based coordinate system requires limited global information, such as the solution of single-phase flow. Numerical results show that one can achieve accurate upscaling results using a flow-based coordinate system.     

A phase-field method for the direct simulation of two-phase flows in pore-scale media using a non-equilibrium wetting boundary condition

Advances in pore-scale imaging (e.g., μ-CT scanning), increasing availability of computational resources, and recent developments in numerical algorithms have started rendering direct pore-scale numerical simulations of multi-phase flow on pore structures feasible. Quasi-static methods, where the viscous and the capillary limit are iterated sequentially, fall short in rigorously capturing crucial flow phenomena at the pore scale. Direct simulation techniques are needed that account for the full coupling between capillary and viscous flow phenomena. Consequently, there is a strong demand for robust and effective numerical methods that can deliver high-accuracy, high-resolution solutions of pore-scale flow in a computationally efficient manner. Direct simulations of pore-scale flow on imaged volumes can yield important insights about physical phenomena taking place during multi-phase, multi-component displacements. Such simulations can be utilized for optimizing various enhanced oil recovery (EOR) schemes and permit the computation of effective properties for Darcy-scale multi-phase flows.We implement a phase-field model for the direct pore-scale simulation of incompressible flow of two immiscible fluids. The model naturally lends itself to the transport of fluids with large density and viscosity ratios. In the phase-field approach, the fluid-phase interfaces are expressed in terms of thin transition regions, the so-called diffuse interfaces, for increased computational efficiency. The conservation law of mass for binary mixtures leads to the advective Cahn–Hilliard equation and the condition that the velocity field is divergence free. Momentum balance, on the other hand, leads to the Navier–Stokes equations for Newtonian fluids modified for two-phase flow and coupled to the advective Cahn–Hilliard equation. Unlike the volume of fluid (VoF) and level-set methods, which rely on regularization techniques to describe the phase interfaces, the phase-field method facilitates a thermodynamic treatment of the phase interfaces, rendering it more physically consistent for the direct simulations of two-phase pore-scale flow. A novel geometric wetting (wall) boundary condition is implemented as part of the phase-field method for the simulation of two-fluid flows with moving contact lines. The geometric boundary condition accurately replicates the prescribed equilibrium contact angle and is extended to account for dynamic (non-equilibrium) effects. The coupled advective Cahn–Hilliard and modified Navier–Stokes (phase-field) system is solved by using a robust and accurate semi-implicit finite volume method. An extension of the momentum balance equations is also implemented for Herschel–Bulkley (non-Newtonian) fluids. Non-equilibrium-induced two-phase flow problems and dynamic two-phase flows in simple two-dimensional (2-D) and three-dimensional (3-D) geometries are investigated to validate the model and its numerical implementation. Quantitative comparisons are made for cases with analytical solutions. Two-phase flow in an idealized 2-D pore-scale conduit is simulated to demonstrate the viability of the proposed direct numerical simulation approach.     

Development of a Darcy-Brinkman model to simulate water flow and tracer transport in a heterogeneous karstic aquifer (Val d’Orléans, France)

Darcy’s law is the equation of reference widely used to model aquifer flows. However, its use to model karstic aquifers functioning with large pores is problematic. The physics occurring within the karstic conduits requires the use of a more representative macroscopic equation. A hydrodynamic model is presented which is adapted to the karstic aquifer of the Val d’Orléans (France) using two flow equations: (1) Darcy’s law, used to describe water flow within the massive limestone, and (2) the Brinkman equation, used to model water flow within the conduits. The flow equations coupled with the transport equation allow the prediction of the karst transfer properties. The model was tested by using six dye tracer tests and compared to a model that uses Darcy’s law to describe the flow in karstic conduits. The simulations show that the conduit permeability ranges from 5?×?10^?6 to 5.5?×?10^?5?m² and the limestone permeability ranges from 8?×?10^?11 to 6?×?10^?10?m². The dispersivity coefficient ranges from 23 to 53 m in the conduits and from 1 to 5 m in the limestone. The results of the simulations carried out using Darcy’s law in the conduits show that the dispersion towards the fractures is underestimated.     

Numerical simulation of two-phase flow in conceptualized fractures

Two-phase flow in fractured rock is an important phenomenon related to a range of practical problems, including non-aqueous phase liquid contamination of groundwater. Although fractured rocks consist of fracture networks, the study of two-phase flow in a single fracture is a pre-requisite. This paper presents a conceptual and numerical model of two-phase flow in a variable fracture. The void space of the fracture is conceptualized as a system of independent channels with position-dependent apertures. Fundamental equations, governing two-phase displacement in each channel, are derived to represent the interface positions and fractional flows in the fracture. For lognormal aperture distributions, simple approximations to fractional flows are obtained in analytical form by assuming void occupancy based on a local capillary allowability criterion. The model is verified by analytical solutions including two-phase flow in a parallel-plate fracture, and used to study the impacts of aperture variation, mobility ratio and fracture orientation on properties of two-phase flow. Illustrative examples indicate that aperture variation may control the distribution of wetting and non-wetting fluids within the fracture plane and hence the ability of the fracture to transmit these fluids. The presence of wetting fluid does little to hinder non-wetting fluid flow in fractures with large aperture variations, whereas a small volume of non-wetting fluid present in the fracture can significantly reduce wetting fluid flow. Large mobility ratios and high fracture slope angles facilitates the migration of non-wetting fluid through fractures.     

Discretization of flow contain nitrate in porous media by finite volume technique

One of the most important pollutants of groudwaters is nitrate. Different human activities including the application of chemical fertilizers in agriculture, causes the emission of nitrate into groudwaters. In this paper, the dynamic effect of soil moisture on carbon and nitrogen cycles has been analyzed by presenting a connection between soil moisture sample and nonlinear differential equations. At present, wide researches are carried out on modeling soil moisture control in solution flows contain nitrate. In order to do so, separation of energy conservation law equations is carried out by a particular method. The mathematical model governing the nitrate containing current in non-isotropic environment has been presented in the form of combined equations. Equation for distribution in multiple environments and Darcy rule has been considered in this model. Then, using finite volume method, separation of flows contain nitrate in porous media is carried out. The current flux is obtained from central difference approximations or upwind approximation. Mashad plain has been considered for case study at this research. Carrying out calibration operation, the measured results have been contrasted with numerical results of finite volume method. After testing the model, it is possible to foresee the way of nitrate changes in other nodes of calculation network. Using these forecasts, the quality of drinking water for several next years is determined. Carrying out numerical modeling by finite volume method, it is found out that the quality of drinking water of Mashad plain would be suitable for the next ten years.     

Uncoupling of coupled flows in soil—a finite element method

Coupled flow of water, chemicals, heat and electrical potential in soil are of significance in a variety of circumstances. The problem is characterized by the coupling between different flows, i.e. a flow of one type driven by gradients of other types, and by the dual nature of certain flows, i.e. combined convection and conduction. Effective numerical solutions to the problem are challenged due to the coupling and the dual nature. In this paper, we first present a general expression that can be used to represent various types of coupled flows in soil. A finite element method is then proposed to solve the generalized coupled flows of convection-conduction pattern. The unknown vector is first decomposed into two parts, a convective part forming a hyperbolic system and a conductive part forming a parabolic system. At each time step, the hyperbolic system is solved analytically to give an initial solution. To solve the multi-dimensional hyperbolic system, we assume that a common eigenspace exists for the coefficient matrices, so that the system can be uncoupled by transforming the unknown vector to the common eigenspace. The uncoupled system is solved by the method of characteristics. Using the solution of the hyperbolic system as the initial condition, we then solve the parabolic system by a Galerkin finite element method for space discretization and a finite difference scheme for time stepping. The proposed technique can be used for solving multi-dimensional, transient, coupled or simultaneous flows of convection-conduction type. Application to a flow example shows that the technique indeed exhibits optimality in convergence and in stability.     

An interface-condition substitution strategy for theoretical study of dissolution-timescale reactive infiltration instability in fluid-saturated porous rocks

A theoretical study of reactive infiltration instability is conducted on the dissolution timescale. In the present theoretical study, the transient behavior of a dissolution-timescale reactive infiltration system needs to be considered, so that the upstream region of the chemical dissolution front should be finite. In addition, the chemical dissolution front of finite thickness should be considered on the dissolution timescale. Owing to these different considerations, it is very difficult, even in some special cases, to derive the first-order perturbation solutions of the reactive infiltration system on the dissolution timescale. To overcome this difficulty, an interface-condition substitution strategy is proposed in this paper. The basic idea behind the proposed strategy is that although the first-order perturbation equations in the downstream region cannot be directly solved in a purely mathematical manner, they should hold at the planar reference front, which is the interface between the upstream region and the downstream region. This can lead to two new equations at the interface. The main advantage of using the proposed interface-condition substitution strategy is that through using the original interface conditions as a bridge, the perturbation solutions for the dimensionless acid concentration, dimensionless Darcy velocity, and their derivatives involved in the two new equations at the interface can be evaluated just by using the obtained analytical solutions in the upstream region. The proposed strategy has been successfully used to derive the dimensionless growth rate, which is the key issue associated with the theoretical study of dissolution-timescale reactive infiltration instability in fluid-saturated porous rocks.     

双相介质中纵波方程的高阶有限差分解法

从双相介质中的纵波方程出发,导出了求解双相各向同性介质中纵波方程的高阶差分格式,给出了吸收边界条件和稳定性条件,在此基础上实现了双相各向同性介质中纵波方程的高阶有限差分法正演模拟,数值模拟结果表明,这种算法能在少量增加计算量的前提下大大提高精度,算法可同时应用于叠前和叠后的数值模拟。     

Using a shifted Grünwald-Letnikov scheme for the Caputo derivative to study anomalous solute transport in porous medium

We propose an extension of the shifted Grünwald-Letnikov method to solve fractional partial differential equations in the Caputo sense with arbitrary fractional order derivative α and with an advective term. The method uses the relation between Caputo and Riemann-Liouville definitions, the shifted Grünwald-Letnikov, and the traditional backward and forward finite difference method. The stability of the method is investigated for the implicit and explicit scheme with homogeneous boundary conditions, and a stability criterion is found for the advective-dispersive equation. An application of the method is used to solve contaminant diffusion and advective-dispersive problems. The numerical solution for the fractional diffusion and fractional advection-dispersion is compared with their respective analytical solutions for different time and space grid refinements. The diffusion simulation exhibited a good fit between the analytical and numerical solutions, with the explicit scheme going from stable to unstable as the time and space refinement changes. The fractional advection-dispersion application produced small deviations from the analytical solution. These deviations, however, are analogous to the numerical dispersions encountered in conventional finite difference solutions of the advection-dispersion equation. The new method is also compared with the traditional L2 method. Notably, an example that involves asymmetrical fractional conditions, a fractional diffusivity that depends on time, and a source term show how the methods compare. Overall, this study assesses the quality and easiness of use of the numerical method.     

Computations with finite element methods for the Brinkman problem

Various finite element families for the Brinkman flow (or Stokes–Darcy flow) are tested numerically. Particularly, the effect of small permeability is studied. The tested finite elements are the MINI element, the Taylor–Hood element, and the stabilized equal order methods. The numerical tests include both a priori analysis and adaptive methods.     

Sequential implicit vertex approximate gradient discretization of incompressible two-phase Darcy flows with discontinuous capillary pressure

This work deals with sequential implicit schemes for incompressible and immiscible two-phase Darcy flows which are commonly used and well understood in the case of spatially homogeneous capillary pressure functions. To our knowledge, the stability of this type of splitting schemes solving sequentially a pressure equation followed by the saturation equation has not been investigated so far in the case of discontinuous capillary pressure curves at different rock type interfaces. It will be shown here to raise severe stability issues for which stabilization strategies are investigated in this work. To fix ideas, the spatial discretization is based on the Vertex Approximate Gradient (VAG) scheme accounting for unstructured polyhedral meshes combined with an Hybrid Upwinding (HU) of the transport term and an upwind positive approximation of the capillary and gravity fluxes. The sequential implicit schemes are built from the total velocity formulation of the two-phase flow model and only differ in the way the conservative VAG total velocity fluxes are approximated. The stability, accuracy and computational cost of the sequential implicit schemes studied in this work are tested on oil migration test cases in 1D, 2D and 3D basins with a large range of capillary pressure parameters for the drain and barrier rock types. It will be shown that usual splitting strategies fail to capture the right solutions for highly contrasted rock types and that it can be fixed by maintaining locally the pressure saturation coupling at different rock type interfaces in the definition of the conservative total velocity fluxes. The numerical investigation of the sequential schemes is also extended to the widely used finite volume Two-Point Flux Approximation spatial discretization.

     

Numerical well model for non-Darcy flow through isotropic porous media[*]This work was supported in part by the EPA under grant # R 825207-01-1, by the State of Texas under ARP/ATP grant # 010366-168 and by a gift grant from Mobil Oil Corp.

Numerical simulation of fluid flow in a hydrocarbon reservoir has to account for the presence of wells. The pressure of a grid cell containing a well is different from the average pressure in that cell and different from the bottom-hole pressure for the well [17]. This paper presents a study of grid pressures obtained from the simulation of single phase flow through an isotropic porous medium using different numerical methods. Well equations are proposed for Darcy flow with Galerkin finite elements and mixed finite elements. Furthermore, high velocity (non-Darcy) flow well equations are developed for cell-centered finite difference, Galerkin finite element and mixed finite element techniques.     

Three dimensional simulation of fluid flow in X‐ray CT images of porous media

A numerical scheme is developed in order to simulate fluid flow in three dimensional (3‐D) microstructures. The governing equations for steady incompressible flow are solved using the semi‐implicit method for pressure‐linked equations (SIMPLE) finite difference scheme within a non‐staggered grid system that represents the 3‐D microstructure. This system allows solving the governing equations using only one computational cell. The numerical scheme is verified through simulating fluid flow in idealized 3‐D microstructures with known closed form solutions for permeability. The numerical factors affecting the solution in terms of convergence and accuracy are also discussed. These factors include the resolution of the analysed microstructure and the truncation criterion. Fluid flow in 2‐D X‐ray computed tomography (CT) images of real porous media microstructure is also simulated using this numerical model. These real microstructures include field cores of asphalt mixes, laboratory linear kneading compactor (LKC) specimens, and laboratory Superpave gyratory compactor (SGC) specimens. The numerical results for the permeability of the real microstructures are compared with the results from closed form solutions. Copyright © 2004 John Wiley & Sons, Ltd.     

Adjoint analysis of Buckley-Leverett and two-phase flow equations

This paper analyzes the adjoint equations and boundary conditions for porous media flow models, specifically the Buckley-Leverett equation, and the compressible two-phase flow equations in mass conservation form. An adjoint analysis of a general scalar hyperbolic conservation law whose primal solutions include a shock jump is initially presented, and the results are later specialized to the Buckley-Leverett equation. The non-convexity of the Buckley-Leverett flux function results in adjoint characteristics that are parallel to the shock front upstream of the shock and emerge from the shock front downstream of the shock. Thus, in contrast to the behavior of Burgers’ equation where the adjoint is continuous at a shock, the Buckley-Leverett adjoint, in general, contains a discontinuous jump across the shock. Discrete adjoint solutions from space-time discontinuous Galerkin finite element approximations of the Buckley-Leverett equation are shown to be consistent with the derived closed-form analytical solutions. Furthermore, a general result relating the adjoint equations for different (though equivalent) primal equations is used to relate the two-phase flow adjoints to the Buckley-Leverett adjoint. Adjoint solutions from space-time discontinuous Galerkin finite element approximations of the two-phase flow equations are observed to obey this relationship.     

TP or not TP,that is the question

We give here a comparative study on the mathematical analysis of two (classes of) discretization schemes for the computation of approximate solutions to incompressible two-phase flow problems in homogeneous porous media. The first scheme is the well-known finite volume scheme with a two-point flux approximation, classically used in industry. The second class contains the so-called approximate gradient schemes, which include finite elements with mass lumping, mixed finite elements, and mimetic finite differences. Both (classes of) schemes are nonconforming and can be expressed using discrete function and gradient reconstructions within a variational formulation. Each class has its specific advantages and drawbacks: monotony properties are natural with the two-point finite volume scheme, but meshes are restricted due to consistency issues; on the contrary, gradient schemes can be used on general meshes, but monotony properties are difficult to obtain.