Combining the Modified Discrete Element Method with the Virtual Element Method for Fracturing of Porous Media

TL;DR

This paper introduces a novel coupling of the Modified Discrete Element Method with the Virtual Element Method to improve hydraulic fracturing simulations in porous media, enabling better boundary condition handling and multi-scale modeling.

Contribution

It establishes a new connection between MDEM and VEM, providing a basis for coupling discrete and continuum models on polyhedral grids for hydraulic fracturing.

Findings

Coupled MDEM-VEM method successfully simulates hydraulic fracturing.

The approach allows multi-scale modeling with complex reservoir geometries.

Implementation in MRST demonstrates practical applicability.

Abstract

Simulation of fracturing processes in porous rocks can be divided into two main branches: (i) modeling the rock as a continuum which is enhanced with special features to account for fractures, or (ii) modeling the rock by a discrete (or discontinuous) approach that describes the material directly as a collection of separate blocks or particles, e.g., as in the discrete element method (DEM). In the modified discrete element (MDEM) method, the effective forces between virtual particles are modified in all regions, without failing elements, so that they reproduce the discretization of a first order finite element method (FEM) for linear elasticity. This provides an expression of the virtual forces in terms of general Hook's macro-parameters. Previously, MDEM has been formulated through an analogy with linear elements for FEM. We show the connection between MDEM and the virtual element method (VEM), which is a generalization of FEM to polyhedral grids. Unlike standard FEM, which computes strain-states in a reference space, MDEM and VEM compute stress-states directly in real space. This connection leads us to a new derivation of the MDEM method. Moreover, it gives the basis for coupling (M)DEM to domain with linear elasticity described by polyhedral grids, which makes it easier to apply realistic boundary conditions in hydraulic-fracturing simulations. This approach also makes it possible to combine fine-scale (M)DEM behavior near the fracturing region with linear elasticity on complex reservoir grids in the far-field region without regridding. To demonstrate the simulation of hydraulic fracturing, the coupled (M)DEM-VEM method is implemented using the Matlab Reservoir Simulation Toolbox (MRST) and linked to an industry-standard reservoir simulator.

Equations

$$\nabla \cdot \sigma$$

$$\epsilon$$

$$\sigma$$

$$[{\epsilon }_{11},{\epsilon }_{22},{\epsilon }_{33},\sqrt{2}{\epsilon }_{23},\sqrt{2}{\epsilon }_{13},\sqrt{2}{\epsilon }_{12}{]}^T$$

$$\sigma =2\mu \epsilon +\lambda \mathrm{tr}(\epsilon )I.$$

$$\alpha :\beta =\mathrm{tr}({\alpha }^T\beta )=\sum _{i,j=1}^3{\alpha }_{i,j}{\beta }_{i,j},$$

$$\sigma =2\mu \epsilon +\lambda \mathrm{tr}(\epsilon )I+\kappa (\tau -\varphi )$$

$$\tau =\frac{1}{2}(\nabla u-{\nabla }^{T}u).$$

$$E={\int }_{\Omega }(\mu \epsilon :\epsilon +\frac{\lambda }{2}\mathrm{tr}(\epsilon {)}^2+\frac{\kappa }{2}(\tau -\varphi ):(\tau -\varphi ))dx.$$

$$\delta E={\int }_{\Omega }(\nabla \cdot (2\mu \epsilon +\lambda \mathrm{tr}(\epsilon )I+\kappa (\tau -\varphi ))\cdot \delta u)dx+{\int }_{\Omega }\kappa (\tau -\varphi ):\delta \varphi dx,$$

$$\nabla \cdot (2\mu \epsilon +\lambda \mathrm{tr}(\epsilon )I+\kappa (\tau -\varphi ))=0$$

$$\tau -\varphi =0.$$

$$\Delta X^m=X^{p_2}-X^{p_1},\Delta U^m=U^{p_2}-U^{p_1},I^m=\frac{\Delta X^m}{\left|\Delta X^m\right|}.$$

$$F_n^m=k_n\Delta U_n^m{F}_s^m=k_s\Delta U_s^m.$$

$$\Delta U_n^m$$

$$\Delta U_s^m$$

$$F^m=F_n^m+F_s^m.$$

$$F^m=\sigma I^{m}ds,$$

$$\sigma =\frac{1}{V}\sum _{m=1}^{N_c}{d}_m{F}^m\otimes I^m,$$

$$\Delta U^m=\ell \Delta X^m=(e+r)\Delta X^m,$$

$$F_n^m=d_m{k}_n(I^m\cdot e{I}^m+I^m\cdot r{I}^m)I^m=d_m{k}_n(I^m\cdot e{I}^m)I^m,$$

$$F_s^m=d_m{k}_s(e{I}^m-I^m\cdot e{I}^m+(r-\theta )\times I^m).$$

$$\sigma=\frac{1}{V}\sum_{m=1}^{N_{c}}d_{m}^{2}\big{(}(k_{n}-k_{s})(I^{m}\cdot eI^{m})\,[I^{m}\otimes I^{m}]\\ +k_{s}\,\left([(eI^{m})\otimes I^{m}]+[((r-\theta)\times I^{m})\otimes I^{m}]\right)\big{)}$$

$$\sigma =2(k_n-k_s)\mathrm{tr}(e)I+2k_{s}e+2k_s(r-\theta ).$$

$$\lambda =2(k_n-k_s)\text{and}\mu =k_s.$$

$$\mu =k_n+k_s\lambda =k_n-k_s$$

$$\nu =\frac{\lambda }{2(\mu +\lambda )}=\frac{1}{4}(1-\frac{k_s}{k_n})<\frac{1}{4}$$

$$\nu =\frac{\lambda }{2\mu +\lambda }=\frac{k_n-k_s}{3k_n+k_s}=\frac{1}{3}\frac{1-\frac{k_s}{k_n}}{1+\frac{k_s}{k_n}},$$

$$\sum _{m=1}^{N_c}{F}^m\times X^m=\sum _{m=1}^{N_c}{d}_m^2{k}_s((r-\theta )\times I^m)\times I^m=\sum _{m=1}^{N_c}{d}_m^2{k}_s((r-\theta )-(r-\theta )\cdot I^m{I}^m)$$

$$\sum _{m=1}^{N_c}{F}^m\times X^m=4d^2{k}_s(r-\theta ).$$

$$U=M\epsilon$$

$$E_{mdem}=\frac{1}{2}{U}^{T}KU=\frac{1}{2}\epsilon M^{T}KM\epsilon ,$$

Full text

1. Introduction

Effective control of flows in geological formations is a key factor for exploiting resources that are highly important for the society, such as ground water, geothermal energy, geological storage of CO2, high quality fossil fuel (gas and oil) and potentially large scale storage of energy in terms or heat or gas. Today 60% of the world energy consumption is based on oil and natural gas resources [22]. In addition 19% is based on coal which needs large scale CO2 storage to be safely exploited without large scale impact on climate [33]. Geothermal energy is an important source of green energy, which would be even more valuable in the future as the supply of fossil fuel is expected to decrease. Gas storage is today an integrated part of the energy supply and provides reliable large scale storage of energy. It enables to both attenuate the volatility of energy prices and ensure energy security.

The use of all of the above resources will benefit from a reliable control of the flow properties around the wells that are used to exploit them. Increased injectivity is particularly important for exploiting resources in tight formation or where high flow rates are required. For enhanced geothermal applications, rock fracturing is a prerequisite for economical exploitation. For CO2 injection, where large volumes of fluid have to be injected, high injectivity limits the increase in pressure near the well and simplifies the operation. Tight formations contain much of the hydrocarbon reserves. The exploitation of these formations has been a driving force for the technology of fracking, which is a more drastic well stimulation than the traditional ones. When increasing the injectivity in a well, it is of vital importance to be able to predict and control fracturing to avoid unwanted fractures or even induced seismic events, which may cause environmental damage as well as the disruption of operations. The key for enabling high injectivity is to induce and control the fracturing process using the coupling between fluid flow, heat and rock mechanics. The failure of the rock and the propagation of fractures depend on both global and local effects, through the stress distribution which is intrinsically global and failure criterias which are local. It is therefore important to have flexible simulation tools that are able to cover both large scale features with complex geometry and include specific fracture dynamics where the fracturing processes occur.

Numerical methods for simulating fracturing can typically be classified either as continuum or discontinuum based methods [23]. The modeling of fracturing in brittle materials like rock is particularly difficult. It is determined by the stress field in the vicinity of the fracture tip. As [15] showed, failures happen when the global energy release is larger than the energy required by the fracturing process. The first depends on the global stress field while the latter is associated with the energy needed at small scale to create a fracture. For brittle materials where the failure happens at very small scales, linear elasticity governs the behavior in most of the domain but the solution of linear elasticity in the presence of fractures is singular near the fracture tip (see [25] for a general description). This introduces challenges for numerical calculations and often results in artificial grid dependence of the simulated dynamics. From a physical point of view, such effects are removed when plasticity is introduced, however the length scale of this region may be prohibitively small to be resolved numerically on the original model. Several techniques have been introduced to incorporate the singularity at the fracture tip into the numerical calculations explicitly, for example specific tip elements in the finite element FEM method. In general, the methods using global energy arguments are less sensitive to the choice of numerical methods than those that use estimates of the strength of the singularity [25]. In the case of the fracturing of natural rock, the uncertainty in the model is large, small scale heterogeneities are important and several different fracturing mechanisms complicate the structure. Discrete modeling techniques have been very successful in this area [28], in particular if complex behaviors should be simulated.

An essential component to model fracturing is therefore the ability to account in a flexible manner for both large and small scale behaviors. This is reflected in the widespread use of tools based on analytic models for hydraulic fracturing (for a review see [26]). However, it becomes a challenge to incorporate interaction of fractures and fine-scale features into those simulation tools. Because of their simplicity and flexibility for incorporating different fracturing mechanisms, discrete element methods (DEM), also called in their explicit variants, distinct element methods, have been one of the main techniques used for hydraulic fracturing in commercial simulators. Those methods exploit the ability of easily modifying the interactions and connections between the discrete particles or elements. For continuum models, such behavior is more difficult to account for. However, the parameters in the DEM model are not directly related to physical macro-scale parameters and there are restrictions on the range of the parameters that can be simulated. In particular, [4] showed that only Poisson’s ratios (in plain stress) smaller than $1/3$ can be considered. The modified discrete element method (MDEM) was introduced to get rid of this restriction and gives also easier relationships between the macro parameters in the linear elastic domain, while keeping the advantages of DEM in the treatment of fracture. In this work, we show the connection between MDEM and the recent development of Virtual Element Methods (VEM). Such approach provides a simple derivation of the MDEM framework, and also highlight the discrepancy of the original DEM model from linear elasticity. The linear version of the VEM methods for elasticity can be used to extend first-order FEM on simplex grids, which was the basis of the MDEM method, to general polyhedral grids. We use the fact that both DEM, MDEM and VEM share the same degrees of freedom in the case of simplex grids to derive smooth couplings between these methods. A similar approach has been followed for coupling FEM with DEM previously [34, 38]. The introduction of VEM opens the possibility for flexible gridding on general polyhedral grids in the far-field region while keeping the DEM/MDEM flexibility in the near fracture domains. Geological formations are typically the result of deposition and erosion processes, which lead to layered structures and faults. Geometrical models using polyhedral grids, such as Corner Point Grid [36], Skua Grid [16] and Cut-Cell [29] are natural in this context and correspond to grids used in the industry of flow modeling in reservoirs. Our proposed method therefore may simplify the incorporation of fracture simulation in realistic subsurface applications.

2. Methods

We study the methods for the standard equations of linear elasticity given by

[TABLE]

where $\sigma$ is the Cauchy stress tensor, $\varepsilon$ the infinitesimal strain tensor and $u$ the displacement field. The linear operator $C$ is the fourth-order Cauchy stiffness tensor. In Kelvin notation [19], a three-dimensional symmetric tensor $\{\varepsilon_{ij}\}$ is represented as an element of $\mathbb{R}^{6}$ with components

[TABLE]

while a two-dimensional symmetric tensor is represented by a vector in $\mathbb{R}^{3}$ given by $[\varepsilon_{11},\varepsilon_{22},\sqrt{2}\varepsilon_{12}]^{T}$. Using this notation $C$ can be represented by a $6\times 6$ matrix and inner products of tensors correspond to the normal inner-product of vectors. For isotropic materials, we have the constitutive equations

[TABLE]

where $\mu$ and $\lambda$ denote the Lamé constants. The elastic energy density is given by $\frac{1}{2}\sigma:\varepsilon$ where we use the standard scalar product for matrices defined as

[TABLE]

for any two matrices $\alpha,\beta\in\mathbb{R}^{3\times 3}$.

2.1. Discete element method

Discrete element methods consist of modeling the mechanical behaviour of a continuum material by representing it as a set of particles, or discrete elements. The forces in the material are then modeled as interaction forces between the particles. There are several variants of the discrete element method [28]. Here we will use the simple version introduced in [4] where the particles are discs in 2D and spheres in 3D. We will also restrict the treatment to the linear case to compare with linear elasticity, but this is not a restriction of the method. For more in-depth presentation of different variants see [28] or [10] and the references therein.

The starting point of the DEM methods has its background in the description of granular media. This has a long history starting from the description of the contact force by Hertz and [30]. In this field an important question was to study how the effective elastic modulus of the bulk was related to the microscopic description [12, 40]. In DEM the basic ideas is to use a microscopic description to simulate the behavior of the bulk modulus.

For a complete relation between general DEM method and linear elasticity using shear forces, it is necessary to introduce the material laws for micropolar materials, see [4]. This introduces an extra variable associated with local rotation, as illustrated in Figure 1. For an isotropic micropolar media the strain stress relation is

[TABLE]

where the extra variable $\phi$ describe the local rotation and $\tau$ represent the asymmetric part of the strain tensor, i.e, rigid body rotations. In terms of displacements it can be written

[TABLE]

The state variable are the displacement $u$ and the local rotation $\phi$ and the total elastic energy is given by

[TABLE]

By computing the variation of $E$,

[TABLE]

we obtain the governing equations of the system, that is, the linear momentum conservation equation,

[TABLE]

We will use the expression of the stress given by (3) to compare with the DEM model which include local rotations. Let us consider two particles $p_{1}$ and $p_{2}$ which are connected through a contact denoted $m$. For a particle $p_{i}$, $i=1,2$, we denote by $X^{p_{i}}$ the position of the particle. The degrees of freedom of the system are the displacement $U^{p_{i}}$ and the microrotation $\theta^{p_{i}}$ for each particle $p_{i}$. Furthermore we set

[TABLE]

We introduce also the distance between the particles, $d_{m}=\left|\Delta X^{m}\right|$. We use the cross-product to represent the action of a rotation so that the rotation given by a vector $\theta$ is the mapping given by $X\to\theta\times X$. Our description of DEM follows [4] with slight differences in the notation. We introduce the normal and shear forces,

[TABLE]

For a given contact $m$ the relative shear and normal displacement are given by

[TABLE]

where $\theta^{m}=\frac{1}{2}(\theta^{p_{1}}+\theta^{p_{2}})$. Note that in the case where two adjacent spheres roll one over the other without sliding, we have $\theta^{p_{1}}=-\theta^{p_{2}}$ so that the term $\Delta U_{s}^{m}$ accounts only for the sliding part of the tangential component. Let us define the total force over a contact $m$ as

[TABLE]

Using the definition of stress tensor $\sigma$, we have that, at the contact between the spheres of center $p_{2}$ and $p_{1}$, and assuming that there exists a non-zero contact surface $ds$, the force $F^{m}$ can be written as

[TABLE]

as $I^{m}$ points in the normal direction. The Cauchy’s formula for the stress matrix [4], which is meant to invert this relation, is given by

[TABLE]

where $N_{c}$ denotes the number of contact points, that is the number of spheres in contact. Let us consider a linear deformation $\ell$ and write $\Delta U^{m}$ as

[TABLE]

where the tensor $e$ and $r$ are respectively the symmetric and skew-symmetric parts of $\ell$. Since $r$ is a skew-symmetric matrix, it corresponds to a rotation and, abusing the notations, we will write indifferently $r\Delta X^{m}$ or $r\times\Delta X^{m}$ to denote the same rotation operator (here applied to $\Delta X^{m}$). To proceed with the identification of the stress tensor, we assume small displacement, that is, $e$ and $r$ are small compared with the identity, and we assume also $\theta^{m}=\theta$ for some constant $\theta$. We use (8) and obtain

[TABLE]

as $r$ is skew-symmetric. For $F_{s}^{m}$, we have

[TABLE]

Hence, we obtain the following expression for the stress tensor,

[TABLE]

To illustrate the restriction that this expression imposes on the parameters, we consider a square packing in 3D. In this case, there are $N_{c}=6$ contact points and, using $I^{1}=-I^{2}=(1,0,0)$, $I^{3}=-I^{4}=(0,1,0)$ and $I^{5}=-I^{6}=(0,0,1)$, we obtain

[TABLE]

Note that we do not take $V$ equal to the volume of the sphere but $V=1$, that is the effective volume. In the expression above, $r$ and $\theta$ must be seen as matrices and not as vectors as in (10). We can also identified the parameter $\phi$ of local rotation for micropolar media with the local rotation $\theta$ in the DEM model. This gives the following Lamé coefficients,

[TABLE]

Hence, as $\sigma:e=2\mu\sum_{ij}e_{ij}^{2}+\lambda\sum_{i}e_{ii}^{2}=2k_{s}\sum_{i\neq j}e_{ij}^{2}+2k_{n}\sum_{i}e_{ii}^{2}$, we can conclude that, for square lattices, this is only stable if $k_{s}>0$ (we assume $k_{n},k_{s}\geq 0$). However, this is not a restriction for simplex grids. Using the same approach as above but now for regular simplices, it is shown in [4] that

[TABLE]

Since $k_{s}$ and $k_{n}$ are naturally positive, this restricts the Poisson’s ratio to

[TABLE]

in the 3D case. For the 2D case, we obtain the same expression in the case of plane strain boundary conditions and, in the case of plane stress, we get

[TABLE]

which implies $-1<\nu<\frac{1}{3}$. These limitations on the physical parameters have been the main motivation for introducing MDEM. Comparing the expression in equation (11) to the governing equations for a micropolar medium (6), we see that the conservation of torque is equivalent to the conservation of angular momentum. Indeed, we get from (9) and (10) that

[TABLE]

so that for a square lattice ($I^{1}=-I^{2}=(1,0,0)$, $I^{3}=-I^{4}=(0,1,0)$ and $I^{5}=-I^{6}=(0,0,1)$), we get

[TABLE]

and the requirement that the torque is zero yields $r-\theta=0$, which corresponds to the conservation of angular momentum equation (6b). This also highlights the need for introducing rotational degrees of freedom for the DEM method if shear forces are used. If not, one gets the non physical effects that rigid rotations introduce forces. Notice that the method which here is referred as DEM is a specific version of a lattice model where the edges of a simplex grid are used to calculate force and the normal force is independent of the rotation of the particles. The last statement could be understood as neglecting rolling resistance.

We also notice that the introduction of angles has been made in finite element literature for membrane problems [5, 9, 20]. In this context, it is called the ”drilling degree of freedom”, see [13] for review. The motivation has been to remove the singularity of the stiffness matrix and the angular degree of freedom as a stiffening effect on the structure. In fact in [21] the value for the free parameter associated with the non symmetric part in the variation principle is recommended to be $\mu$, the shear modulo. The degrees of freedom are completely the same as in DEM.

2.2. Modified discrete element method MDEM

The motivation for the introduction of the MDEM method is twofold. First, in DEM, the relation between macro parameters and the parameters is not simple. Secondly, given a configuration of particles, it is not possible to reproduce all the parameters associated with isotropic materials as discussed in the paragraph above. The same type of restrictions also holds for hexahedral and square grids, see [39, 35]. In [31], thermodynamical considerations are used to show that, for isotropic materials, the value of the Poisson’s ratio should satisfy $\nu>1/5$. In this perspective, the restriction $\nu<1/4$ established above for DEM appears very restrictive. The ability to vary the mechanical properties even for this configuration introduces non central forces between the particles, in this context called shear forces. As discussed above this can only be done if extra local rotation variables are introduced. This has two disadvantages, first it is more complicated, and secondly the final system is equivalent to a micopolar medium and not a purely elastic medium. Restricting oneself to central forces may therefore be in some cases preferable but one should remain aware that such assumption comes with very strong restriction on the material parameters. In [18], the authors consider the Cauchy relations which are known to be necessary for an elastic material where only central forces are present, each atom or molecule is a center of symmetry and the interaction forces are well approximated by an harmonic potential. They show that, for an isotropic material, the Cauchy relations imply that $\nu=1/3$. This very strong restriction makes it difficult to consider models only based on central forces.

The basic idea of MDEM is to use an interaction region, instead of looking at the forces on each particle as a result of interaction with neighboring particles like in DEM, see Figure 2. Then, the force at a particle is given by the sum of the forces computed at the particle for each interactive region the particle belongs to. In the finite element setting, the interaction region corresponds to an element and a particle to a node. The calculation of the forces is equivalent to the case of linear finite element. The original derivation [3, 2], was based on explicit representation of the geometry and calculation of forces. Here we will base our derivation on the variational form of linear elasticity. To simplify the derivation we will use the fact that for simplex grids there exists a one to one mapping between non rigid-body linear deformations and the length of the edges. By non rigid-body linear deformations, we mean the quotient space of the space of linear deformations with the space of rigid-body deformations (translation and rotations). Such space is in bijection with symmetric matrices, the strain tensors. Using the notation of [2] but with all tensors represented in the Kelvin notation where the tensor inner product reduces to normal inner-product. For simplices one can relate the non zero strain states to edge length $U$ by

[TABLE]

Note that, to simplify the expressions, we use different notations in the previous section where $U$ was denoted by $\Delta U$. We write the energy of the element as

[TABLE]

where $K$ is a symmetric definite positive matrix to be determined. The tensor $K$, which we will call in the paper the MDEM stiffness tensor, depends on the material parameter. This fulfills the requirement of linear elasticity that rigid motion does not contribute to the energy. The normal forces can be calculate as the generalized forces associated with the variable U, that is

[TABLE]

From (17), we can see that assuming that only central forces are present and the shear forces are negligible is equivalent to the requirement that $K$ is diagonal. Using the analogue definition of stress where we exploit the kelvin notation

[TABLE]

If we consider the energy of the same system for a linear elastic media assuming constant stress and strain, which is the case for linear elements on simplex grids, the result is

[TABLE]

where $D$ is the representation of the forth order stiffness tensor $C$ in Kelvin notation. Note also that $epsi$ in (19) is meant either as a tensor (in the first equality) or as a vector written in Kelvin notations (in the second). For the sake of simplicity, we will continue to do the same abuse of notations in the following. We see that one reproduces the energy of a linear elastic media if

[TABLE]

which gives that

[TABLE]

The difference between the matrix $K$ used in DEM and the matrix needed to reproduce linear elasticity used in MDEM is that the latter case normally is a full matrix. Since DEM methods solve Newtons’s equation with a dissipation term it will minimize this energy. The same is the case of standard Galerkin discretization of linear FEM on simplices, which by construction have the same energy functional as MDEM. Consequently, the only difference, if no fracture mechanism is present will be the method for computing the solution to the whole system of equations. The DEM methods rewrite the equations in the form of Newton laws with an artificial damping term and let time evolve to converge to the solution, see [10]. For FEM, the linear equations are usually solved directly.

The advantage of using the MDEM formulation compared to FEM is that it offers the flexibility to choose independently on each element if a force should be computed using linear elasticity or if a more traditional DEM calculation should be used.

The ability to associate the edge lengths to the non rigid body motions is only possible for simplex grids. An other important aspect to this derivation is that the degrees of freedom uniquely define all linear motions and no others. The importance of the last part will be more evident after comparison with the VEM method.

2.3. The Virtual Element Method

In contrast to FEM, the Virtual element method seeks to provide consistency up to the right polynomial order of the equation in the physical space. This is done by approximating the bilinear form only using the definition of the degrees of freedom, as described below. The FEM framework on the other hand defines the assembly on reference elements, using a set of specific basis functions. This however has disadvantages for general grids where the mappings may be ill defined or complicated. VEM avoids this problem by only working in physical space using virtual elements and not computing the Galerkin approximation of the bilinear form exactly. This comes with a freedom in the definition of the method and a cost in accuracy measured in term of the energy norm.

As the classical finite element method, the VE method starts from the linear elasticity equations written in the weak form of Equation 1,

[TABLE]

We have also introduced the symmetric gradient $\varepsilon$ given by

[TABLE]

for any displacement $\bm{u}$. The fundamental idea in the VE method is to compute on each element an approximation $a_{K}^{h}$ of the bilinear form

[TABLE]

that, in addition of being symmetric, positive definite and coercive with respect to the non rigid-body motions, it is also exact for linear functions. The correspondence between MDEM and VEM we study here holds only for a first-order VEM method. When higher order methods are used, the exactness must hold for polynomials of a given degree where the degree determines the order of the method. These methods were first introduced as mimetic finite element methods but later developed further under the name of virtual element methods (see [11] for discussions). The degrees of freedom are chosen as in the standard finite element methods to ensure the continuity at the boundaries and an element-wise assembly of the bilinear forms $a_{K}^{h}$. We have followed the implementation described in [14]. In a first-order VE method, the projection operator $\mathcal{P}$ into the space of linear displacement with respect to the energy norm has to be computed locally for each cell. The VE approach ensures that the projection operator can be computed exactly for each basis element. The projection operator is defined with respect to the metric induced by the bilinear form $a_{K}$. The projection is self-adjoint so that we have the following Pythagoras identity,

[TABLE]

for all displacement field $\bm{u}$ and $\bm{v}$ (in order to keep this introduction simple, we do not state the requirements on regularity which is needed for the displacement fields). In [14], an explicit expression for $\mathcal{P}$ is given so that we do not even have to compute the projection. Indeed, we have $\mathcal{P}=\mathcal{P}_{R}+\mathcal{P}_{C}$ where $\mathcal{P}_{R}$ is the projection on the space $R$ of translations and pure rotations and $\mathcal{P}_{C}$ the projection on the space $C$ of linear strain displacement. The spaces $R$ and $C$ are defined as

[TABLE]

Then, the discrete bilinear form $a_{K}^{h}$ is defined as

[TABLE]

where $s_{K}$ is a symmetric positive matrix which is chosen such that $a_{K}^{h}$ remains coercive. Note the similarities between (25) and (24). Since $\mathcal{P}_{R}$ and $\mathcal{P}_{C}$ are orthogonal and $\mathcal{P}_{R}$ maps into the null space of $a_{K}$ (rotations do not produce any change in the energy), we have that the first term on the right-hand side of (24) and (25) can be simplified to

[TABLE]

The expression (25) immediately guarantees the consistency of the method, as we get from (25) that, for linear displacements, the discrete energy coincides with the exact energy. Since the projection operator can be computed exactly for all elements in the basis - and in particular for the virtual basis elements for which we do not have explicit expressions - the local matrix can be written only in terms of the degrees of freedom of the method. In our case the degrees of freedom of the method are the value of displacement at the node. Let us denote $\bm{\varphi}_{i}$ a basis for these degrees of freedom. The matrix $(A_{K})_{i,j}=a_{K}^{h}(\bm{\varphi}_{i},\bm{\varphi}_{j})$ is given by

[TABLE]

In (26), $W_{C}$ is the projection operator from the values of node displacements to the space of constant shear strain and $S_{K}$, which corresponds to a discretization of $S_{K}$ in (25), is a symmetric positive matrix which guarantees the positivity of $A_{K}$. There is a large amount of freedom in the choice of $S_{K}$ but it has to scale correctly. We choose the same $S_{K}$ as in [14]. The matrix $D$ in (26) corresponds to the tensor $C$ rewritten in Kelvin notations so that, in three dimensions, we have

[TABLE]

Finally, the matrices $A_{K}$ are used to assemble the global matrix $A$ corresponding to $a^{h}$. In this paper, we use the implementation available as open source through the Matlab Reservoir Simulation Toolbox (MRST) [32]. The approach of splitting the calculation of the energy in terms of a consistent part block on one side and a higher order block one the other side was also used in the free formulation of finite elements [6]. In this case the motivation was to find an alternative element formulation,

2.4. Correspondences between VEM and MDEM

For simplex grids the regularization term $S_{K}$ in the expression for the local stiffness matrix in Equation 26 is zero because in this case the projection operator is equal to the identity. If we introduce the operator from the degrees of freedom for the element to the edge expansions $R$, we can compare the two expression for the local energy,

[TABLE]

One easily identifies the operators $W_{C}$ and $M^{-1}R^{T}$ as the projection operator$\mathcal{P}_{c}$ to the non-rigid body motions represented in Kelvin notation type of symmetric strain. The degrees of freedom span exactly the space of linear displacement and do not excite any higher order modes with nonzero energy. An illustration of the different concepts is given in Figure 2. We point that both DEM and VEM calculate the basic stiffness matrix in real space, contrary to most FEM methods which do this on the reference element. When dealing with simplex grids the advantage of using the DEM method within an explicit solving strategy (often called distinct element method) is that the calculation of the edge length extensions $U$ can be calculated for each edge, and only the matrices $M^{-1}$ and the Cauchy stiffness tensor $C$ are needed locally. These matrices only operate on the small space of non rigid motion with dimension $((d(d+1))/2)$ while the operator $W_{c}$ works on the all the deformations which have dimension $((d+1)d)$. The edge length can thus be seen as an efficient compact representation of the non rigid motions, which holds only on simplices. As we have seen, both MDEM and VEM can be derived from the calculation of the energy in each element. For the linear elastic part it is not necessary to introduce extra angular degrees of freedom. However, this may be needed for certain DEM methods. In this case, we refer to the use of drilling elements in combination with the use of the free formulation of FEM [13, 7], which, as discussed earlier, shares some fundamental ideas with VEM, such as energy orthogonality (which corresponds (24)), and rc-modes exactness (whic) and the freedom in choosing the stabilization term.

2.5. Fluid mechanics coupling

We introduce a coupling with a fluid flow through the Biot’s equations [8]. The Biot’s equations are given by

[TABLE]

where $S_{c}p+\alpha\nabla\cdot\bm{u}$ denotes the fluid content. The fluid content depends on the storativity $S_{c}$, the fluid pressure $p$ and on the rock volume change given by $\nabla\cdot\bm{u}$ which is weighted by the Biot-Willis constant $\alpha$. In (28), $K$ denotes the permeability and $\mu_{v}$ the fluid viscosity. For flow, and in particular if multiphase behaviors are considered, the most successful methods have been based on finite volume methods. The basic time discretization using the two point flux method or multi point flux methods [1] can be written as

[TABLE]

Here $\text{grad}_{p}$ is a discrete gradient operator from cell pressures to face, $\text{div}_{f}$ is the corresponding discrete divergence acting on face fluxes. The source term $Q$ represents the injection of fluids, see [24] for more details on those discrete operators.

Given an implicit time discretization the coupling term in the Biot case requires a discrete divergence operator $\text{div}_{d}$ for displacement field. Note that this discrete operator can be implemented exactly for first order VEM, see [37] for more details. The semi-discrete equations are

[TABLE]

Here $A_{s}$ system matrix of the mechanical system, $\text{div}_{d}$ is the divergence operator acting on the nodal displacement and gives a volume expansion of a cell and $\alpha$ is the Biot parameter depending on the ratio between the rock and fluid compressibility. In the context of MDEM when the simulation of fracturing is the main purpose, we normally approximate only the volume expansion term in the transport equation for the fractured cells, where the expansion is also the largest. Except for this term an explicit update of pressure is used. This approximation also avoids problems due to small permeabilities which can cause numerical locking and artificial oscillations in the fluid pressure, see [17].

2.6. Fracturing criteria

In the MDEM method, before an element is fractured, it behaves as in FEM and the MDEM stiffness tensor $K$ is obtained from the Cauchy stiffness tensor through the relation established in (20). Depending on the physical situation a fracturing criteria based on stress is used, for example Mohr-Coulomb. In the examples in this paper we will use the simple tensile failure criteria, namely

[TABLE]

where $\sigma_{tens}$ is the tensile stress. After failure, we use a central force model, where the forces are calculated individually for each edge as

[TABLE]

where $K_{d}(\Delta U)$ denotes the diagonal matrix such that

[TABLE]

If a fracture is closing, then the effective force will in this case be as for DEM using only central forces. As for all methods trying to simulate fracturing, the critical point is how to avoid grid dependent fracturing, due to the singularity of the stress field near the fracture front. In this work however, the main aim is to see how the far-field solution can be simulated using general grid, independently of the fracturing modeling.

2.7. Solution method

The solution method in MDEM is chosen to be similar to the one used in DEM. That is explicit time integration of Newton’s laws. To get fast convergence to the physical stationary state, the local artificial damping term that can be found [10] is often preferable. This is not a physical damping mechanism, but it avoids large differences in local time steps restrictions. The advantage of this approach is that it is less sensitive to global changes than the traditional FEM approach which solves directly the stationary state by solving the linearized equations. This is particularly important when discontinuous changes of the forces due to changes in the medium is present. For MDEM this is the case for the situation of initial fracturing, equation (31) or in contact properties for fracture cells, equation (32). The result in all cases is that the forces are discontinuous with respect to the degrees of freedom. Explicit methods have been shown to have advantages for such problems even if the main dynamics is globally elliptic, because the non-linearities in the problem impose stronger time-step size requirement that those needed for the explicit integration of the elliptic part. As the damping criteria depends on the concept of total nodal forces, it can also be used on the nodes connected with VEM type of force calculations. No other modification apart from the force calculations are needed.

3. Examples

We demonstrate the features of the presented framework with two examples. First we show how the effective parameters of linear elasticity in simple DEM with only normal forces depend on the particular choice of the grid cells. Second, we use VEM, MDEM and DEM on a general a polyhedral grid to demonstrate how this can be combined within a uniform framework.

When a fracture has occurred in a cell, but the whole system evolves in such a way that the fracture closes again, then we should have forces normal to the fracture faces and, depending on the fracture model, forces along the fracture. Here, we choose to model this by an effective stiffness tensor. Indeed, we keep using the DEM model (and solver) after the fracture closes, meaning that the materiel parameters for the cell are given by the diagonal MDEM stiffness tensor $K_{d}$ as defined in (33). From Section 2.4, we know that it also corresponds to a unique Cauchy stiffness tensor. Let us study the effect of such choice and measure the difference between the original and this post-fracturing stiffness tensor. If we denote by $\mathcal{K}$, the one-to-one transformation from the MDEM stiffness tensor $K$ to the Cauchy stiffness tensor $C$, we compute, for a given $C$, the difference between $C$ and $\mathcal{K}^{-1}(\mathcal{K}(C)_{d})$. We consider a equilateral triangle and an isotropic material with Young’s modulus $E=1$ and Poison ratio $\nu=1/4$. For this value of $\nu$ and this shape, the matrix $K$ is diagonal, so that $C=\mathcal{K}^{-1}(\mathcal{K}(C)_{d})$. This reference triangle is plotted in yellow in Figure 3. We keep the same Cauchy stiffness tensor but modify the shape of the triangle by translating one of the corners. For each configuration that is obtained, we get a different post-fracturing MDEM stiffness tensor given by $\mathcal{K}(C)_{d}$ and we plot the six component of the tensor $T$ defined as the difference $T=C-\mathcal{K}^{-1}(\mathcal{K}(C)_{d})$. We notice that the changes in $C_{2,3}$ and $C_{1,3}$ is zero when on the line $x=0$ in the figure. This show that in this case as expected the effective model has biaxial symmetry. We also notice that there is quite strong changes in the effective parameters even for relatively small changes in the triangles. It should also be noted that a break of the edges along the x-axis, which in the MDEM fracture model result in putting one of the corresponding diagonal element to zero only changes the value of $C_{1,1}$. This is because this only acts in the x direction.

We use MRST [32, 27] to generate the unstructured grid presented in Figure 4 and set up an example which combines the use of VEM for the general cell shapes and the use of MDEM for the triangular cells that can easily be switched to a DEM model when a fracture is created. The total grid size is 30\text{\,}\mathrm{m}$\times$30\text{\,}\mathrm{m}. In the middle within a diameter $0.5\text{\,}\mathrm{m}$ we have placed cells associated with a well. The permeability used was $10\text{\,}\mathrm{nd}$, porosity of $0.3$, the compressibility of the fluid is similar to water, $1\text{\times}{10}^{-10}\text{\,}{\mathrm{Pa}}^{-1}$, and it is injected fluid corresponding to the pore volume of all well cells in an hour. The solution is shown after $40$ minutes. The initial condition was given by the mechanical solution with a force of $1\text{\times}{10}^{7}\text{\,}\mathrm{Pa}$ at the top and rolling all other places. The initial condition for pressure is constant pressure equal to $1\text{\times}{10}^{7}\text{\,}\mathrm{Pa}$. The mechanical parameters are given by $E=$1\text{\times}{10}^{9}\text{,}\mathrm{Pa} and $\nu=0.3$. The well cells are set to have Young’s modulo $E=$1\text{\times}{10}^{4}\text{\,}\mathrm{Pa} and finally the tensile strength is $2\text{\times}{10}^{5}\text{\,}\mathrm{Pa}$. We observe that the fracture propagates in the direction so that the fracture plane (or line in 2D) is aligned with the maximum stress plane (or line in 2D). We get slight grid orientation effect since there is no way a planar fault in the y direction can be obtain using the given triangular grid. The interface between the grids has large steps in cell sizes and include hanging nodes, but no effects due to these features are observed as long as the fracture does not reach the interface. Near the tip of the fracture we observe oscillation of the stress on cells, which is a well known problem for first order triangular elements. However the values associate with the nodes is better approximated and patch recovery techniques [41] can be used to get better stress fields as seen in Figure 5. A note is that the dynamics of DEM or MDEM, is associated with the sum of all forces from all elements around a node, not individual stresses for cells.

4. Conclusions

In this paper we have shown how MDEM and VEM for linear elasticity share the same basic idea of projection to the states of linear non-rigid motions, although with different representations, length extension for MDEM and polynomial basis for VEM. Both are equivalent to linear FEM on simplices, but the viewpoint presented here gives a more direct way on how they relate. Since both share the same degrees of freedom, except possibly the angular degree of freedom of MDEM/DEM, we combine these methods with minimal implementation issues. This is used to simulate fracture growth, where the near field regions is described by a simplex grid which is suited for DEM and MDEM, while the general polyhedral grids is used in the far-field region. The coupling between the grids which can contain hanging nodes and significant changes in cell shapes and sizes, can be done without introducing large errors. We see this method as a valuable contribution to flexible coupling of MDEM/DEM methods with traditional reservoir modeling grids.

5. Acknowledgments

This publication has been produced with support from the KPN project Controlled Fracturing for Increased Recovery. The authors acknowledge the following partners for their contributions: Lundin and the Research Council of Norway (244506/E30).