Finite difference techniques for modelling geothermal reservoirs

Two finite difference algorithms suitable for long-time simulation of the exploitation of a two-phase geothermal reservoir are presented. One is based on the hopscotch method proposed by Saul'yev¹ and analysed further by Gordon² and Gourlay.³ The other is based on the well-known ADI method. Both methods use a Newton–Raphson iterative technique in order to obtain accurate solutions of the non-linear difference equations involved. Rapid convergence of the iterative schemes occurs both for single-phase and two-phase reservoir problems. One- and two-dimensional model problems are presented.     

水气二相流方程的一种离散数值解法

提出了水气二相流方程的一种数值解法.在利用有限元方法求解水气二相流方程时,引入了离散Newton迭代方法,用于非线性有限元方程组的线性化处理,将这一步计算的收敛阶由原有研究的线性收敛提高到平方收敛,并避免了直接应用Newton迭代方法给编程带来的不便.同时在求解两相的有限元方程组时,采用两相方程组并行迭代的方法,与联立计算相比节省了大量的内存空间.     

Mathematical and numerical filtration–advection–dispersion model of miscible grout propagation in saturated porous media

The development of a predictive model of behaviour of porous media during injection of miscible grout, taking into account convection, dilution and filtration of grout solution with interstitial water, as well as consolidation aspects, is presented. Model assumptions are reviewed and discussed first. During the establishment of the model, we insist on surface terms and their physical relevance in expressing adsorption effects. Constitutive laws such as Fick's law for diffusive mass transport, hydrodynamic dispersion tensor dealing with miscibility, are modified by taking into account filtration effects. A new surface term appears in mass balance equations as a consequence of filtration. According to the filtration laws used, an initial filtration rate is estimated on the basis of a one‐dimensional experimental campaign. The field equations are discretized by using Galerkin finite element and θ‐scheme standard method. For transport equation, Streamline Upwind Petrov Galerkin method is employed to prevent numerical oscillations. Lastly, confrontation of numerical results with laboratory experiments constitutes a first step to validate the model on a realistic basis. Copyright © 2001 John Wiley & Sons, Ltd.     

Numerical study on finite element implementation of hypoplastic models

A comprehensive numerical study on finite element implementation of hypoplastic models is presented. Two crucial aspects, local integration of the constitutive equations (the local problem) and forming tangent operators for Newton–Raphson iteration (the global problem), are investigated. For solving the local problem, different integration algorithms, including explicit and implicit methods, are examined using tri-axial compression tests and incremental stress response envelopes, as well as typical boundary value problems. For solving global problems, three different ways of generating the tangent operator are compared. The numerical evidences indicate that, in terms of accuracy, efficiency and robustness, explicit methods with substepping and error control are the best choices for constitutive integration of hypoplastic models while the so-called continuum tangent operators have certain advantages over two other types of numerically-generated consistent tangent operators.     

Geothermal heat extraction by water circulation through a large crack in dry hot rock mass

One proposed geothermal heat extraction scheme relies on water circulation in a large vertical crack created by hydraulic fracturing in a hot dry impermeable rock mass. Water flow, heat convection and crack opening widths are analysed by finite elements. Governing field equations of the problem are first set up rigorously and then various small terms are identified and neglected, retaining the effects of pressure gradient, buoyancy, velocity head (kinetic energy) and head loss due to viscous friction in the water flow equation, and the effects of heat convection in water and heat conduction in rock in the heat transfer equation. The finite element scheme for water flow is based on a variational principle that is typical for diffusion problems, and for heat transfer it is based on the method of least-square residuals. The system of differential equations is highly non-linear. The non-linear terms and coefficients are treated in the fiaite element analysis as constant; the finite element analysisof, the steady-state pressures, fluxes and temperatures is then iterated, evaluating all non-linear terms and coefficients on the basis of the solution obtained in the previous iteration. Numerically calculated fields at various times after the start ofcooling are presented. They indicate some features favourable for the geothermal scheme, such as formation of eddy currents, and downward flux of water toward hotter rock. However, other important questions would have to be solved to gain full understanding, of this proposed geothermal scheme.     

Coupled HM analysis using zero‐thickness interface elements with double nodes. Part I: Theoretical model

In recent years, the authors have proposed a new double‐node zero‐thickness interface element for diffusion analysis via the finite element method (FEM) (Int. J. Numer. Anal. Meth. Geomech. 2004; 28 (9): 947–962). In the present paper, that formulation is combined with an existing mechanical formulation in order to obtain a fully coupled hydro‐mechanical (or HM) model applicable to fractured/fracturing geomaterials. Each element (continuum or interface) is formulated in terms of the displacements (u) and the fluid pressure (p) at the nodes. After assembly, a particular expression of the traditional ‘u–p’ system of coupled equations is obtained, which is highly non‐linear due to the strong dependence between the permeability and the aperture of discontinuities. The formulation is valid for both pre‐existing and developing discontinuities by using the appropriate constitutive model that relates effective stresses to relative displacements in the interface. The system of coupled equations is solved following two different numerical approaches: staggered and fully coupled. In the latter, the Newton–Raphson method is used, and it is shown that the Jacobian matrix becomes non‐symmetric due to the dependence of the discontinuity permeability on the aperture. In the part II companion paper (Int. J. Numer. Anal. Meth. Geomech. 2008; DOI: 10.1002/nag.730 ), the formulation proposed is verified and illustrated with some application examples. Copyright © 2008 John Wiley & Sons, Ltd.     

Dynamics of unsaturated poroelastic solids at finite strain

We derive the governing equations for the dynamic response of unsaturated poroelastic solids at finite strain. We obtain simplified governing equations from the complete coupled formulation by neglecting the material time derivative of the relative velocities and the advection terms of the pore fluids relative to the solid skeleton, leading to a so‐called u^s ? p^w ? p^a formulation. We impose the weak forms of the momentum and mass balance equations at the current configuration and implement the framework numerically using a mixed finite element formulation. We verify the proposed method through comparison with analytical solutions and experiments of quasi‐static processes. We use a neo‐Hookean hyperelastic constitutive model for the solid matrix and demonstrate, through numerical examples, the impact of large deformation on the dynamic response of unsaturated poroelastic solids under a variety of loading conditions. Copyright © 2011 John Wiley & Sons, Ltd.     

A consistent formulation of a dilatant interface element

The paper considers a plane joint or interface element suitable for implementation into a standard non-linear finite element code. The element is intended to model discontinuities with rough contact surfaces, such as rock joints, where dilatant behaviour is present. Of particular concern is the formulation of a constitutive model which fully caters for all possible histories of opening, closing and sliding (accompained by dilation or contraction) in any direction. The non-linear incremental constitutive equations are formulated in a manner appropriate for a back-ward difference discretization in time along the path of loading. The advantage of such an approach is that no essential distinction need be drawn between opening, closing and sliding. Further, a convenient formulation of the constitutive equations is facilitated by representing the different contact conditions in relative displacement space. The state diagram in relative displacement space, however, changes from one time step to the next, and evolution equations for the updating must be formulated. These concepts are illustrated for two rock-joint models: a sawtooth asperity model and a limited dilation model. The models are based on a penalty formulation to enforce the contact constraints, and explicit equations for the tangent stiffness matrix and for the corrector step of the standard Newton–Raphson iterative algorithm are derived. These equations have been implemented as an user element into the finite element code ABAQUS⁷. Three examples are presented to illustrate the predictions of the formulation.     

Finite element modelling of transport of organic chemicals through soils

Aquifer contamination by organic chemicals in subsurface flow through soils due to leaking underground storage tanks filled with organic fluids is an important groundwater pollution problem. The problem involves transport of a chemical pollutant through soils via flow of three immiscible fluid phases: namely air, water and an organic fluid. In this paper, assuming the air phase is under constant atmospheric pressure, the flow field is described by two coupled equations for the water and the organic fluid flow taking interphase mass transfer into account. The transport equations for the contaminant in all the three phases are derived and assuming partition equilibrium coefficients, a single convective – dispersive mass transport equation is obtained. A finite element formulation corresponding to the coupled differential equations governing flow and mass transport in the three fluid phase porous medium system with constant air phase pressure is presented. Relevant constitutive relationships for fluid conductivities and saturations as function of fluid pressures lead to non-linear material coefficients in the formulation. A general time-integration scheme and iteration by a modified Picard method to handle the non-linear properties are used to solve the resulting finite element equations. Laboratory tests were conducted on a soil column initially saturated with water and displaced by p-cymene (a benzene-derivative hydrocarbon) under constant pressure. The same experimental procedure is simulated by the finite element programme to observe the numerical model behaviour and compare the results with those obtained in the tests. The numerical data agreed well with the observed outflow data, and thus validating the formulation. A hypothetical field case involving leakage of organic fluid in a buried underground storage tank and the subsequent transport of an organic compound (benzene) is analysed and the nature of the plume spread is discussed.     

Dynamics of unsaturated soils using various finite element formulations

Unsaturated soils are three‐phase porous media consisting of a solid skeleton, pore liquid, and pore gas. The coupled mathematical equations representing the dynamics of unsaturated soils can be derived based on the theory of mixtures. Solution of these fully coupled governing equations for unsaturated soils requires tremendous computational resources because three individual phases and interactions between them have to be taken into account. The fully coupled equations governing the dynamics of unsaturated soils are first presented and then two finite element formulations of the governing equations are presented and implemented within a finite element framework. The finite element implementation of all the terms in the governing equations results in the complete formulation and is solved for the first time in this paper. A computationally efficient reduced formulation is obtained by neglecting the relative accelerations and velocities of liquid and gas in the governing equations to investigate the effects of fluid flow in the overall behavior. These two formulations are used to simulate the behavior of an unsaturated silty soil embankment subjected to base shaking and compared with the results from another commonly used partially reduced formulation that neglects the relative accelerations, but takes into account the relative velocities. The stress–strain response of the solid skeleton is modeled as both elastic and elastoplastic in all three analyses. In the elastic analyses no permanent deformations are predicted and the displacements of the partially reduced formulation are in between those of the reduced and complete formulations. The frequency of vibration of the complete formulation in the elastic analysis is closer to the predominant frequency of the base motion and smaller than the frequencies of vibration of the other two analyses. Proper consideration of damping due to fluid flows in the complete formulation is the likely reason for this difference. Permanent deformations are predicted by all three formulations for the elastoplastic analyses. The complete formulation, however, predicts reductions in pore fluid pressures following strong shaking resulting in somewhat smaller displacements than the reduced formulation. The results from complete and reduced formulations are otherwise comparable for elastoplastic analyses. For the elastoplastic analysis, the partially reduced formulation leads to stiffer response than the other two formulations. The likely reason for this stiffer response in the elastoplastic analysis is the interpolation scheme (linear displacement and linear pore fluid pressures) used in the finite element implementation of the partially reduced formulation. Copyright © 2008 John Wiley & Sons, Ltd.     

Accelerated initial stiffness schemes for elastoplasticity

Iterative methods for the solution of non‐linear finite element equations are generally based on variants of the Newton–Raphson method. When they are stable, full Newton–Raphson schemes usually converge rapidly but may be expensive for some types of problems (for example, when the tangent stiffness matrix is unsymmetric). Initial stiffness schemes, on the other hand, are extremely robust but may require large numbers of iterations for cases where the plastic zone is extensive. In most geomechanics applications it is generally preferable to use a tangent stiffness scheme, but there are situations in which initial stiffness schemes are very useful. These situations include problems where a nonassociated flow rule is used or where the zone of plastic yielding is highly localized. This paper surveys the performance of several single‐parameter techniques for accelerating the convergence of the initial stiffness scheme. Some simple but effective modifications to these procedures are also proposed. In particular, a modified version of Thomas' acceleration scheme is developed which has a good rate of convergence. Previously published results on the performance of various acceleration algorithms for initial stiffness iteration are rare and have been restricted to relatively simple yield criteria and simple problems. In this study, detailed numerical results are presented for the expansion of a thick cylinder, the collapse of a rigid strip footing, and the failure of a vertical cut. These analyses use the Mohr–Coulomb and Tresca yield criteria which are popular in soil mechanics. Copyright © 2000 John Wiley & Sons, Ltd.     

Numerical simulation of a stratigraphic model

A multi-lithology diffusive stratigraphic model is considered, which simulates at large scales in space and time the infill of sedimentary basins governed by the interaction between tectonics displacements, eustatic variations, sediment supply, and sediment transport laws. The model accounts for the mass conservation of each sediment lithology resulting in a mixed parabolic, hyperbolic system of partial differential equations (PDEs) for the lithology concentrations and the sediment thickness. It also takes into account a limit on the rock alteration velocity modeled as a unilaterality constraint. To obtain a robust, fast, and accurate simulation, fully and semi-implicit finite volume discre tization schemes are derived for which the existence of stable solutions is proved. Then, the set of nonlinear equations is solved using a Newton algorithm adapted to the unilaterality constraint, and preconditioning strategies are defined for the solution of the linear system at each Newton iteration. They are based on an algebraic approximate decoupling of the sediment thickness and the concentration variables as well as on a proper preconditioning of each variable. These algorithms are studied and compared in terms of robustness, scalability, and efficiency on two real basin test cases.     

Seepage flow problems by variational inequalities: Theory and approximation

The theory of variational inequalities enables us to formulate and solve free boundary problems in fixed domains, while most other methods assume the position of the unknown domain in solving the problem. Here the problem of seepage flow through a rectangular dam with a free boundary is formulated as a vertical inequality following the ideas of Baiocchi. In order to demonstrate the essential ideas of extending the domain of the solution of problems with free boundaries, the problem of the deflection, of a string on a rigid support is first examined. Next, variational inequalities are derived which are associated with several cases of seepage problems. An approximation theory, including a priori error estimates, is developed using finite element methods, and an associated numerical scheme is given. It is shown that for linear and quadratic finite element methods, the rates of convergence are 0(h) and 0(h^1.25-δ), 0 < δ < 0.25, respectively, if the permeability is constant.     

An Empirical Method for Analysis of Load Transfer and Settlement of Single Piles

An empirical method is developed for estimating the load transfer and deformation of drilled, in situ formed piles subjected to axial loading. Firstly, governing equations for soil–pile interaction are developed theoretically, taking into account spatial variations in: (a) shaft resistance distribution and (b) ratio of load sharing between the shaft and base. Then generic load transfer models are formulated based on examination of data from 10 instrumented test piles found in the literature. The governing equations and load transfer models are then combined and appropriate boundary conditions defined. Using an incremental-iterative algorithm whereby all the boundary conditions are satisfied simultaneously, a numerical scheme for solving the combined set of equations is developed. The algorithm is then developed into an interactive computer program, which can be used to predict the load-settlement and axial force distribution in piles. To demonstrate its validity, the program is used to analyse four published case records of test piles, which other researchers had analysed using the following three computationally demanding tools: (a) load transfer (t–z), (b) finite difference and (c) finite element methods. It is shown that the proposed method which is much less resource-intensive, predicts both the load-settlement variation and axial force distribution more accurately than methods: (a–c) above.     

Finite element analysis of consolidation of layered clay soils using an elastic visco-plastic model

This paper presents a general one-dimensional (1-D) finite element (FE) procedure for a highly non-linear 1-D elastic visco-plastic (1-D EVP) model proposed by Yin and Graham for consolidation analysis of layered clay soils. In formulating the 1-D FE procedure, a trapezoidal formula is used to avoid the unsymmetry of the stiffness matrix for a Newton (modified Newton) iteration scheme. Unlike many other 1-D FE approaches in which the initial in situ stresses (or stress/strain states) are considered indirectly or even not considered, the initial in situ stress/strain states are taken into account directly in this paper. The proposed FE procedure is used for analysis of 1-D consolidation of a clay with published test results in the literature. The FE modelling results are in good agreement with the measured results. The FE model and procedure is then used to analyse the consolidation of a multi-layered clay soils with a parametric study on the effects of the variations of creep parameters in Yin and Graham's 1-D EVP model. It is found that the creep parameters ψ/V and t₀ have significant influence on the compression and porewater pressure dissipation. For some boundary conditions, changes of parameters in one layer will have some effects on the consolidation behaviour of another layer due to the different consolidation rates. Finally, the importance of initial stress/strain states is illustrated and discussed. Copyright © 1999 John Wiley & Sons, Ltd.     

Steady‐state solution for cylindrical penetrometers

This paper suggests a new method for obtaining steady‐state solutions for ‘full‐flow’ penetrometers. The method is based on the numerical solution of the small strain plastic‐flow problem (i.e. rigid plastic material) with an inhomogeneous strength field, which is determined by converting changes of material properties over time in a stationary frame of reference into spatial distribution of strength in a moving frame of reference. Rather than building streamlines from back integration of soil element distortion, as previous methods have suggested, the method treats the domain as continuous with the associated field equations. The method employs an upstream weighting technique for the determination of information flow within the domain. The execution order for the calculation is based on topological ordering. This results in the calculation having a complexity of O(N), as compared with O(N^1.5) for the strain path or streamline methods (N is the number of discretized points), which significantly reduces the calculation time. The formulation is presented for the cylindrical (T‐bar) penetrometer, and includes aspects of soil strength degradation, strain rate effects, strength anisotropy, and interface strength law. Comparison to previously published values, based on large displacement finite element simulations with remeshing, showed good agreement, indicating on the correctness of the suggested approach. Investigation into the soil rigid‐body rotation and the remolding effect on anisotropy characteristics showed an interesting behavior, where the decrease of strength anisotropy due to remolding has a greater influence when the soil strength is higher in the vertical direction. Copyright © 2009 John Wiley & Sons, Ltd.     

A study on iterative methods for solving Richards’ equation

This work concerns linearization methods for efficiently solving the Richards equation, a degenerate elliptic-parabolic equation which models flow in saturated/unsaturated porous media. The discretization of Richards’ equation is based on backward Euler in time and Galerkin finite elements in space. The most valuable linearization schemes for Richards’ equation, i.e. the Newton method, the Picard method, the Picard/Newton method and the L-scheme are presented and their performance is comparatively studied. The convergence, the computational time and the condition numbers for the underlying linear systems are recorded. The convergence of the L-scheme is theoretically proved and the convergence of the other methods is discussed. A new scheme is proposed, the L-scheme/Newton method which is more robust and quadratically convergent. The linearization methods are tested on illustrative numerical examples.     

A coupled finite element–rigid block method for transient analysis of rock caverns

A coupled finite element–rigid block model for the transient analysis of caverns in jointed media is presented. This coupling permits the modelling of lined openings in a jointed rock mass as well as the propagation of stress waves to the cavern. Both the finite element and the rigid block algorithms employ explicit time integration; an efficient, stable scheme is developed for coupling the two algorithms. Two numerical examples are given: one is a simple validation, the second is a representation of a lined cavern in a sparsely jointed medium.     

Finite Element Solvers for Biot’s Poroelasticity Equations in Porous Media

We study and compare five different combinations of finite element spaces for approximating the coupled flow and solid deformation system, so-called Biot’s equations. The permeability and porosity fields are heterogeneous and depend on solid displacement and fluid pressure. We provide detailed comparisons among the continuous Galerkin, discontinuous Galerkin, enriched Galerkin, and two types of mixed finite element methods. Several advantages and disadvantages for each of the above techniques are investigated by comparing local mass conservation properties, the accuracy of the flux approximation, number of degrees of freedom (DOF), and wall and CPU times. Three-field formulation methods with fluid velocity as an additional primary variable generally require a larger number of DOF, longer wall and CPU times, and a greater number of iterations in the linear solver in order to converge. The two-field formulation, a combination of continuous and enriched Galerkin function space, requires the fewest DOF among the methods that conserve local mass. Moreover, our results illustrate that three out of the five methods conserve local mass and produce similar flux approximations when conductivity alteration is included. These comparisons of the key performance indicators of different combinations of finite element methods can be utilized to choose the preferred method based on the required accuracy and the available computational resources.

     

Finite element consolidation analysis of soils with vertical drain

This paper presents a finite element procedure for the analysis of consolidation of layered soils with vertical drain using general one‐dimensional (1‐D) constitutive models. In formulating the finite element procedure, a Newton–Cotes‐type integration formula is used to avoid the unsymmetry of the stiffness matrix for a Newton (Modified Newton) iteration scheme. The proposed procedure is then applied for the consolidation analysis of a number of typical problems using both linear and non‐linear soil models. Results from this simplified method are compared with those from a fully coupled consolidation analysis using a well‐known finite element package. The average degree of consolidation, excess porewater pressure and average vertical effective stress are almost the same as those from the fully coupled analysis for both the linear and non‐linear cases studied. The differences in vertical effective stresses are tolerable except for the values near the vertical drain boundaries. The consolidation behaviour of soils below a certain depth of the bottom of vertical drain is actually one‐dimensional for the partially penetrating case. Therefore, there are not much differences in whether one uses a one‐dimensional model or a three‐dimensional model in this region. The average degree of consolidation has good normalized feature with respect to the ratio of well radius to external drainage boundary for the cases of fully penetrating vertical drain using a normalized time even in the non‐linear case. Numerical results clearly demonstrate that the proposed simplified finite element procedure is efficient for the consolidation analysis of soils with vertical drain and it has better numerical stability characteristics. This simplified method can easily account for layered systems, time‐dependent loading, well‐resistance, smear effects and inelastic stress–strain behaviour. This method is also very suitable for the design of vertical drain, since it greatly reduces the unknown variables in the calculation and the 1‐D soil model parameters can be more easily determined. Copyright © 2000 John Wiley & Sons, Ltd.