Virtual Element Method for geomechanics on reservoir grids
VVIRTUAL ELEMENT METHOD FOR GEOMECHANICALSIMULATIONS OF RESERVOIR MODELS
ODD ANDERSEN, HALVOR M. NILSEN, AND XAVIER RAYNAUD
Abstract.
In this paper we study the use of Virtual Element Method (VEM)for geomechanics. Our emphasis is on applications to reservoir simulations.The physical processes behind the formation of the reservoirs, such as sedi-mentation, erosion and faulting, lead to complex geometrical structures. Aminimal representation, with respect to the physical parameters of the sys-tem, then naturally leads to general polyhedral grids. Numerical methodswhich can directly handle this representation will be highly favorable, in par-ticular in the setting of advanced work-flows. The virtual element method isa promising candidate to solve the linear elasticity equations on such models.In this paper, we investigate some of the limits of the VEM method whenused on reservoir models. First, we demonstrate that care must be taken tomake the method robust for highly elongated cells, which is common in theseapplications, and show the importance of calculating forces in terms of trac-tion on the boundary of the elements for elongated distorted cells. Second, westudy the effect of triangulations on the surfaces of curved faces, which alsonaturally occur in subsurface models. We also demonstrate how a more stablestabilization term for reservoir application can be derived. a r X i v : . [ m a t h . NA ] F e b ODD ANDERSEN, HALVOR M. NILSEN, AND XAVIER RAYNAUD Introduction
Sedimentary formations are the result of long and complex geological processes.Sedimentation creates thin layers, faulting creates nontrivial connections betweenthe layers and erosion creates degenerate layers. The formation retains an overallstratigraphic structure, in the sense that very different spatial correlations in thematerial properties can be observed between the horizontal and vertical directions,and long and thin cells are specific to reservoir simulations, see the section rep-resent in Figure 1. The geometric modeling of sedimentary formations requires
Figure 1.
Section of the Gullfaks reservoir model (Norway). Eachcolor represents a different material property. We observe how thelarge aspect ratio in the cells.the parameterization of a very large number of complicated interfaces. Each in-terface then separates regions with material properties that may differ of severalorder of magnitude and must be captured with maximum accuracy. Because ofthese difficulties, computational considerations are often not prioritized in the de-sign of geological grids, which will typically contain highly irregular cell shapes.The grid and material properties are strongly related, which cause severe limita-tions on remeshing. The industry standard for reservoir grids is the corner-pointformat [15]. In a corner-point grid, pillars which have a dominant vertical direc-tion are first defined from a two-dimensional Cartesian partition. Then, for eachset of four adjacent pillars, hexahedron cells are constructed by choosing 2 pointson each pillars and connecting all these points (see the detailed in Section 6.2).Many geometrical grid formats have been proposed to improve on this format, forexample Skua Grid [10], S-Grid, Faulted S-Grid and Cut-Cell [12]. By refining themesh, it is of course always possible to improve the quality of the mesh from thepoint of view of numerical computation, but all compact representation of the un-derlying geology, that is a representation where the data (the material properties)is represented by the minimum number of cells, will lead to cells with high aspectratio, distorted cells, faces or cells of very different sizes, cells or faces with differentshapes. Methods which are robust for such grid will greatly simplify the modelingof subsurface physics.In recent years, the coupling of geomechanical effects with subsurface flow hasbecome more and more important in many areas including: oil and gas produc-tion from mature fields, oil and gas production from fractured tight reservoirs,fractured rock for geothermal application and risk assessment of CO injection.Realistic modeling of these applications is hampered by the differences in the waygeomechanics and flow models are build and discretized. Traditionally, the me-chanic problems are solved using finite element methods but they are difficultto adapt to the standard geometrical representation of reservoir models, such as IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS3 corner-point grids. In contrast, the Virtual Element Method (VEM) can operateon general polyhedral grids. As such, the ability of the method to handle irregulargrids makes it very attractive for geomechanical applications. In this paper, weinvestigate if the method can effectively be applied on realistic reservoir grids. Ourmain result concerns the treatment of the load term. We observe that standardstabilization terms presented in the literature are not adapted to elongated cellswith large aspect ratios, which are standard in geological models. We propose amodification of the stabilization constant which can be used in the 2D case and adiscrete gradient approach to compute the load term which turns out to be littlesensitive to the choice of stabilization and can be easily extended to 3D. In a firstpart of this paper, we present the VEM method following mostly [9] but we alsotry to clarify the connection with the construction of the projection operators, asintroduced in the basic principles of [3]. In our numerical experiments, we willfocus on the performance of the VEM method on geological grids. The emphasiswill be on corner-point grids and complex small scale sedimentary models. To beable to demonstrate the VEM method in the reservoir setting we used the MRSTframework [11] to simplify the grid handling. The numerical implementation of theVEM method used in this paper can be downloaded from [13], in particular thetest case concerning the Norne reservoir model (see Section 6.2) is readily availablefrom there. 2.
The equations of linear elasticity
We consider the equation of linear elasticity for small deformations. The dis-placement is given by u ( x ) ∈ R d for d = 2 , x ∈ Ω ⊂ R d .The equations are given by(1) ∇ · σ = f with(2) σ = Cε and ε = 12 ( ∇ + ∇ T ) u Here, σ , ε and u denote the Cauchy stress tensor, the infinitesimal strain tensorsand the displacement field, respectively. The vector function f : R d → R d is anexternal volumetric force that we will refer to as the load term . The linear operator C is a fourth order stiffness tensor which satisfies, for some constant c >
0, theellipticity condition(3) c S : S ≤ S : CS , for any symmetric matrix S ∈ R d × d . The symbol : denotes the scalar product in R d × d defined as(4) A : B = tr( A T B ) , for any A , B ∈ R d × d .3. Presentation of the VEM method for linear elasticity
The VEM method was first introduced in the framework of mimetic discretizationmethods but later rephrased in the language of finite element methods (see [6] fordiscussions). A general presentation of VEM is given in [3]. The same authorspresent convergence results for linear elasticity in [5]. Practical details on theimplementation of VEM are given in [4]. Our implementation of the VEM follows
ODD ANDERSEN, HALVOR M. NILSEN, AND XAVIER RAYNAUD the presentation done in [9] where the specific case of linear elasticity is considered.We rewrite the equation of linear elasticity in the weak form,(5) (cid:90) Ω ε ( v ) : Cε ( u ) d x = (cid:90) Ω f · v d x , which must hold for all displacement field v : R d → R d . Let N c denote the numberof cells and { E i } N c i =1 the grid cells. We define the bilinear form a E i as a E i ( u , v ) = (cid:90) E i ε ( u ) : Cε ( v ) d x and decompose the global bilinear energy form a ( · , · ) in cell contributions,(6) a ( u , v ) := (cid:90) Ω ε ( u ) : Cε ( v ) d x = N c (cid:88) i =1 a E i ( u , v ) . In the rest of this section, we consider a given cell E and will denote by V E thefinite dimensional approximating function space in E . In the VEM approach, thebasis functions of V E are not known explicitly but, for a first-order VEM method,the requirements on V E are that it contains the space of polynomials of order 1,denoted P ( E ), and that the bilinear form a E ( u , v ) can be computed exactly forany u ∈ P ( E ) and any v ∈ V E , using only the degrees of freedom of v . As in thestandard finite element method, the degrees of freedom are the nodal displacements,so that the continuity at the boundaries of each element is ensured by requiringlinearity on the edges and a local reconstruction on the faces, which depends onlyon the values at the edges of the face where the reconstruction is done. The systemmatrix can be assembled element-wise. Let us denote(7) V scalar E = { v ∈ H ( E ) | v | e ∈ P ( e ) for all edges e } , for i = 1 , . . . , d . For a given node η of E , we can construct a function φ η in V scalar E such that φ η (¯ η ) = 1 if ¯ η = η and zero if ¯ η (cid:54) = η . The virtual basis functions of V E are then given by(8) φ kη ( x ) = φ η ( x ) e k for η ∈ N ( E ) and k = 1 , . . . , d , where N ( E ) denotes the set of nodes of the cell E and e k is the unit vector in the direction given by the index k . After havingintroduced the projection operator, we will add some extra requirements for φ η concerning its first and second order moment. But beside that, no more explicitproperties for φ η are needed and this is one of the important point of the method,which also makes it so flexible. The projection operator, which we denote π ∇ , isdefined with respect to the energy norm a E . We consider first order approximationsand, for any displacement field u in the Hilbert space [ H ( E )] , the projection π ∇ ( u ) of u is defined as the element p ∈ [ P ( E )] such that a E ( p , q ) = a E ( u , q ) , for all q ∈ [ P ( E )] . Since the bilinear form a E is degenerate, additional conditionsmust be imposed to define completely π ∇ , see (22) in Section 3.2 for the rigorousdefinition. For any displacement field u , the energy a E ( u , u ) can be decomposedusing Pythagoras’ identity,(9) a E ( u , u ) = a E ( π ∇ u , π ∇ u ) + a E (( I − π ∇ ) u , ( I − π ∇ ) u ) . IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS5
The first term on the right-hand side can be computed exactly from the degree offreedom, for any u ∈ V E . The last term can not be handled generically and istherefore replace by a stabilization term which takes the form of a bilinear form s E ,whose role is to ensure that the ellipticity of a E is retained. Hence, the energy isfinally approximated by(10) a h,E ( u , u ) = a E ( π ∇ u , π ∇ u ) + s E (( I − π ∇ ) u , ( I − π ∇ ) u ) . In the general framework of VEM, as introduced in [3], the computation of theprojection operator typically requires the computation of an inverse, locally foreach cell. The formulation in [9] has the advantage of giving an explicit expressionof the projection operator. In the presentation that follows, we will try to clarifythe connection between the two approaches.3.1.
The kinematics of affine displacement.
The physics of linear elasticityis associated with linear deformations, in particular the rigid body motions play acrucial role. Let us recall some simple facts on the kinematics of affine displace-ments. The linear space of affine displacements, which we denote by P , correspondsto the sum of the translations and linear transformation so that any l ∈ P can bewritten as l ( x ) = u + Lx , for u ∈ R d and L ∈ R d × d . The dimension of the spaceof P is d + d . The subspace of rigid body motion, which we denote P r , containsthe rotation and the translation. Any l ∈ P r can be written as(11) l ( x ) = u + Ω ( x − x ) , for any u , x ∈ R d and Ω that belongs to the space of skew-symmetric matrices,denoted aSym( R d ). There is a redundancy in the choice of x and u so that aunique decomposition of l ∈ P r is given by l ( x ) = u + Ω x for u ∈ R d and Ω skew-symmetric. Hence, the space P r is isomorphic to the sum of the linear spaceof translation and the linear space of skew-symmetric matrices, and its dimensionis therefore d ( d + 1) /
2. The space of non rigid body motion is the quotient of P with respect to P r , which we denote P / P r . We introduce the projection operator π c in P defined as(12) π c ( l ) = 12 ( L + L T )( x − ¯ x E ) , for any l ( x ) = u + Lx ∈ P . Here, ¯ x E denotes the arithmetic average of thepositions x i of all the nodes of the cell E , that is¯ x E = 1 n n (cid:88) i =1 x i , where n corresponds to the number of nodes in E . We can check that π c is aprojection and π c ( l ) = 0 if and only if l ∈ P r . Hence, the image of π c is in bijectionwith the space of linear strain P / P r which we therefore identify to P c = π c ( P ).Note that P c can also be defined as(13) P c = { l ∈ P | ∇ l = ∇ l T and l ( ¯ x E ) = 0 } . Then, we introduce the projection π r from P to P r as(14) π r ( l ) = l − π c ( l )so that π r ( l ) = l ( ¯ x E ) + 12 ( L − L T )( x − ¯ x E ) . ODD ANDERSEN, HALVOR M. NILSEN, AND XAVIER RAYNAUD
We can check that P r = π r ( P ), π c π r = π r π c = 0, π c + π r = I P . The space P c is isomorphic to the space of symmetric matrices, denoted Sym. We consider thecase d = 3 and use Kelvin’s notation to represent Sym so that, for any a ∈ Sym,its Kelvin representation in R , which we denote ˆ a ∈ R , is given by(15) ˆ a T = [ a , a , a , √ a , √ a , √ a ] . The square root in the definition above has the advantage to lead to the followingcorrespondence between the scalar products in Sym and R ,(16) a : b = ˆ a · ˆ b , for any a , b ∈ Sym. Note that the authors in [9] use Voigt instead of Kelvinnotations, which explains why the expressions given in the present paper differ upto a coefficient to those in [9]. We define the symmetric tensor ˆ C ∈ R × by theidentity(17) a : Cb = ˆ a T ˆ C ˆ b , for all a , b ∈ Sym. Then, we obtain that, for any l , m ∈ P , we have(18) a E ( l , m ) = (cid:90) E ε ( l ) : Cε ( m ) d x = | E | (cid:91) π c ( l ) T ˆ C (cid:92) π c ( m )The projection π c ( l ) belongs to P c and can therefore be written as π c ( l ) = Q ( x − ˆ x E ) for Q ∈ Sym. In (18), we wrote (cid:91) π c ( l ) for the Kelvin representation ˆ Q of Q ,and similarly for (cid:92) π c ( m ). For the space P r , we use the mapping between aSymand R given by the cross-product operation. For any a ∈ aSym, we can define arotation vector ˆ a ∈ R as(19) ˆ a = √ − a , , a , , − a , ] T . Then, we have(20) ax = 1 √ a × x and(21) a : b = ˆ a · ˆ b for any a , b ∈ aSym. The normalization using √ in (20) is used in order to getan exact correspondence between the scalar products in aSym and R , as in (16)in the case of Sym. The basis of P r is given by the canonical basis of R , thefirst three components corresponding to the translation vector while the three lastcomponents correspond to a rotation vector.3.2. The projection operator.
The projection operator on the space of affinedisplacement, with respect to the bilinear form a , plays an essential role in theformulation of a first order VEM method. We follow the notation of [3] and denotethis projection by π ∇ , even if the bilinear form we consider is not the H semi-norm, which is used as example in [3]. For a given displacement function ν ∈ V E ,the projection p = π ∇ ( ν ) is defined as the unique element p ∈ P which satisfies(22a) a E ( p , q ) = a E ( ν , q ) , for all q ∈ P and such that(22b) π r p = π r ν , IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS7 which means that the projection of p and ν on P r coincide. The condition (22b) isnecessary to determine a unique solution. Indeed, the bilinear form a is degenerateas it is invariant with respect to the space of rigid body motion and the condition(22b) eliminates this underteminancy by imposing a rigid body motion for p . It isimportant to note that the definition of the projection operator is in fact indepen-dent of V E and could be extended to [ H ( E )] . However, for V E , the projectioncan be computed exactly, using directly the degrees of freedom and without furtherintegration. Up to now, we have only introduced π r and π c on P and we willnow extend these definitions to V E so that the definition in (22b) makes actuallysense. We let the reader check that the new definitions of π c and π r on V E , whenrestricted to P , coincide with those introduced previously. We define the projection π r : V E → P r as(23) π r ( ν ) = ¯ ν E + 12 (cid:10) ∇ ν − ∇ ν T (cid:11) ( x − ¯ x E ) , where the bracket denote the cell average, i.e., (cid:104) w (cid:105) = 1 | E | (cid:90) E w d x and ¯ ν = 1 n (cid:88) E ν i . We define the projection π c : V E → P c as π c ( ν ) = 12 (cid:10) ∇ ν + ∇ ν T (cid:11) ( x − ¯ x E ) . In a moment, we are going to check that both projections can be computed directlyfrom the degree of freedoms. First, we use these definitions to compute π ∇ . Westart by considering a solution p = π ∇ ( ν ) to (22) and show that (22a) yields(24) a E ( p c , q c ) = a E ( ν , q c ) , for any q ∈ P ( E ), where p c = π c ( p ) and q c = π c ( q ). The symmetric gradient iszero for any element in P r , that is ε ◦ π r = 0. Hence, a E ( π r ( p ) , ν ) = 0 for any p ∈ P ( E ) and ν . It implies that a E ( p , q ) = a E ( p c , q c ) and a E ( ν , q ) = a E ( ν , q c ) sothat Equation (22a) indeed implies (24). Let us now determine the element p thatsatisfies (22) for a given ν ∈ V E . The coercivity of the the form a E on P c followsfrom the definition of P c and the coercivity of the tensor C , see (3). Therefore,there exists a unique solution p c ∈ P c such that (24) holds for all q c ∈ P c . For any q c ∈ P c , we have(25) a E ( p c , q c ) = (cid:90) E ∇ p c : C ∇ q c d x = | E | ∇ p c : C ∇ q c and(26) a E ( ν , q c ) = (cid:90) E
12 ( ∇ ν + ∇ ν T ) : C ∇ q c d x = (cid:18) (cid:90) E ( ∇ ν + ∇ ν T ) d x (cid:19) : C ∇ q c . Hence, ∇ p c = (cid:82) E ( ∇ ν + ∇ ν T ) d x which implies that p c is uniquely defined as p c = π c ( ν ). We can conclude that p defined as(27) p = p c + π r ( ν )is the unique solution to (22). Indeed, π c ( p ) = p c and π r ( p ) = π r ( ν ) are bothuniquely defined by (24) and (22b). ODD ANDERSEN, HALVOR M. NILSEN, AND XAVIER RAYNAUD
Let us now give more details on the assembly. To do so, we consider a basisfunction ν i ∈ V E for which the only non-zero displacement can only occur at thenode i , that is ν i ( x j ) = 0 if i (cid:54) = j . Such function can be written as ν i ( x ) = (cid:80) dj =1 ν j φ i ( x ) e j , where { e } dj =1 is the basis for Cartesian coordinates. We have (cid:10) ∇ ν i (cid:11) = (cid:80) dj =1 ν ij e j (cid:104)∇ φ i (cid:105) T , and we have to compute q i = (cid:104)∇ φ i (cid:105) . The expression above simplifies to (cid:10) ∇ ν i (cid:11) = ν i q iT . For q i , using Stokes’ theorem, we have(28) q i = (cid:90) E ∇ φ i d x = (cid:90) ∂E φ i n d x = (cid:88) f ∈ F ( E ) ( (cid:90) f φ i d x ) n f , where F ( E ) denotes the set of faces that belong to E . The integral in (28) canbe computed exactly. For the 3D, we use a virtual space such that the first twomoments of the virtual basis elements coincide with those of their projection, see[1]. The integral is zero if the node i does not belong to the face f and, otherwise,(29) (cid:90) f φ i d x = (cid:40) | f | m + ( n e,i − + n e,i + ) · ( x f − ¯ x f ) in 3D , | f | in 2D , where x f is the centroid of the face f and ¯ x f the arithmetic average of the nodecoordinates, i.e. ¯ x f = m (cid:80) mj =1 x fj . We denote by W ic ∈ R × the matrix represen-tation of π c written in the basis of displacement for the node i (that is R ) and thebasis of P c (that is R , using the Kelvin notation). For l, m = { , , } , we have12 (cid:10) ∇ ν i + ∇ ν iT (cid:11) l,m = 12 ( ν il q im + ν im q il )so that ( W ic ) T = q i √ q i √ q i q i √ q i √ q i q i √ q i √ q i . We have 12 (cid:10) ∇ ν i − ∇ ν iT (cid:11) = 12 ( ν i q iT − q i ν iT ) . Using the general identity ( q i × ν i ) × x = ( q i · x ) ν i − ( ν i · x ) q i , we get that the R representation of the matrix in aSym above is given by √ q i × ν i . Hence, thematrix W ir ∈ R × that represents π r written in the basis of displacement for thenode i (that is R ) and the basis of P r (that is R for the translation and therotation vector) is given by( W ir ) T = n − √ q i √ q i n √ q i − √ q i n − √ q i √ q i . The matrices W c from the space of all the degrees of freedom (that is R n ) to P c isobtained by concatenating W ic and similarly for W r . To obtain, from W c and W r ,the matrix representations of π c and π r in terms only of the degrees of freedom, IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS9 we have to find the decomposition of P c and P r in terms of the degrees of freedom.To do so, we introduce the vectors r i , for i = { , . . . , n } as r i = x i − ¯ x E . We define N ic , N ic ∈ R × as(30) N ic = r i √ r i √ r i r i √ r i √ r i r i √ r i √ r i and N ir = − √ r i √ r i √ r i − √ r i − √ r i √ r i . Then, the matrices N c , N r ∈ R n × are obtained by concatenating N ic , N ir , respec-tively.The projections can be then written in terms as a mapping from degrees offreedom to degrees of freedom,(31) P r = N r W r and P c = N c W c and the projection on affine displacement is given by P = P c + P r . For any ν , η ∈ V E , we have that a E ( π ∇ ν , π ∇ η ) = | E | (cid:92) π c ( ν ) T ˆ C (cid:92) π c ( η ) = | E | ν T W c T ˆ CW c η , where in the last term, slightly abusing the notations, we denote by ν , η ∈ R n the vector composed of the degrees of freedoms of ν , η ∈ V E . Using the sameconvention, we can write the bilinear form as(32) a h,E ( ν , η ) = ν T (cid:16) | E | W c T ˆ CW c + ( I − P ) T S ( I − P ) (cid:17) η , where S ∈ R n × n is a stabilization term. For the VEM method to be well-defined,the matrix S must be chosen such that it is positive, symmetric and definite onthe kernel of P . We note that the decomposition of the energy in two orthogonalparts, the linear part which ensures consistency and the higher order part whichare handled so that stability is preserved, is analog to the decomposition used in[7], even if it was introduced there to add some freedom in the choice of the basisfunctions. 4. Implementation of the load term
The load term can be calculated in several different ways which are equivalent upto the order of accuracy of the methods. We have investigated the three followingalternatives,(1) Computation using the projection operator π ∇ ,(2) Integration using nodal quadrature,(3) Computation based on a discrete gradient operator.Alternative 1 is the choice that naturally follows from the VEM approach and whichis proposed in [3]. Alternative 2 was argued to be simpler and with similar accuracyin [9]. Alternative 3 is possible when the force is equal to the gradient of a potential.This last alternative actually came to the mind of the authors when they consideredthe poro-elasticity equation, where the divergence operator naturally arises. As wewill see below, the discrete gradient is in fact derived from the discrete divergenceoperator by duality. The two first alternatives give similar results. We show thatthe last one has significantly less errors than the others for elongated grid cells. Standard assembly of the load term (alternative 1 and 2).
For a givenforce f , we consider the work done by the force for a given displacement field u ,(33) (cid:90) Ω f · u d x . This expression defines a linear form on the space of displacement. We denote by V the global discrete function space of displacement, which is constructed by takingthe product of the V E for all the cells E of the grid and requiring continuity at thecell boundaries and correspond to the nodal displacement in terms of the degrees offreedom. We want to find a discrete linear form on V that approximates (33). Wecan equip V with the standard scalar product in R n N , that is (cid:80) η u η · ν η . Here, n N denotes the total number of nodes. Any linear form on V can be represented by anelement in V , using this scalar product. Hence, we end up looking for an element (cid:98) f ∈ V such that(34) (cid:90) Ω f · u d x ≈ (cid:88) η (cid:98) f η · u η . The vector (cid:98) f ∈ R N n can be interpreted as a vector of nodal forces . We present sev-eral expressions for (cid:98) f corresponding to the three alternatives presented previously.First, we can use weights which are obtained using a first-order quadrature. For anode η , let us denote by E ( η ) the set of cells to which the node η belongs. Usingquadrature rules to integrate f on each cell, we obtain(35) (cid:98) f η = (cid:88) i ∈ E ( η ) w ηi f ( η ) , see [9] for the definitions of the weights w ηi . This corresponds to alternative 2. Foralternative 1, let u η ∈ V be a displacement for which the only non-zero degrees offreedom are those corresponding to the node η . Then, we have (cid:90) Ω f · u η d x = (cid:88) i ∈ E ( η ) (cid:90) E i f · u η d x ≈ (cid:88) i ∈ E ( η ) (cid:90) E i π i ( f ) · u η d x = (cid:88) i ∈ E ( η ) π i ( f ) · (cid:90) E i u η d x = (cid:88) i ∈ E ( η ) π i ( f ) · (cid:90) E i π ∇ i ( u η ) d x . (36)Here, π i denotes the L projection to the space of constant functions (polynomialsof degree zero) in the element E i . To obtain the last integral, we use that fact thatthe virtual basis functions can be chosen such that the zero and first moment of afunction ν in V E coincide with those of its projection π ∇ i ν , that is (cid:90) E p · ν d x = (cid:90) E p · π ∇ i ( ν ) d x , for any p ∈ P and ν ∈ V E . See [1] for more details. The choice of such basisimplies that, for an element E , the modes that belong to ker π ∇ , typically higher IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS11 nonlinear modes, will not be excited directly by the force. From (36), we infer that (cid:98) f η is defined as a linear combination of cell averages of f ,(37) (cid:98) f η = (cid:88) i ∈ E ( η ) m ηi π i ( f )where m ηi are the weights given m ηi = e k · (cid:90) E i π ∇ i ( φ ηk ) d x . Note that the expression on the right above do not depend on k , as the same basisfunction is used in all directions.4.2. The discrete gradient approach.
Let us now turn to alternative 3. Weassume that the force can be written as the gradient of a potential, f = ∇ ψ . Wehave, for a node η and a dimension k ∈ { , , } ,(38) (cid:90) Ω f · u d x = (cid:90) Ω ∇ ψ · u d x = − (cid:90) Ω ψ ∇ · u d x + (cid:90) ∂ Ω ψ u d x = − (cid:90) Ω ψ ∇ · u d x . The boundary integral vanishes because we assume Dirichlet boundary condition, u = 0 for x ∈ ∂ Ω. In the VEM space, there exists a natural discretization of thedivergence operator as an operator from V to cell-wise constant functions, denoted T , which is isomorphic to R N c , where N c denotes the number of cells. Indeed, forany discretized potential ˆ ψ ∈ T and ν ∈ V , we have(39) (cid:90) Ω ˆ ψ ∇ · ν = (cid:88) E (cid:90) E ˆ ψ E ∇ · ν d x = (cid:88) E (cid:88) f ∈ F ( E ) ˆ ψ E (cid:90) f ν · n d x , where F ( E ) as before denotes the set of faces that belong to E . The last integralcan be computed exactly as shown in (29). Then, using partial integration, we get(40) (cid:90) Ω ˆ ψ ∇ · φ kη d x = (cid:88) j ∈ E ( η ) (cid:88) f j,l ∈ E j ∩ E l ˆ ψ E j ( e k · n j.l ) (cid:90) f j,l φ η d x , with the convention that we only get contribution in the integral when the face f j,l exists, that is when E j and E l share a common face. Note that by definition ofthe exterior normal, we have n j.l = − n j.l . We use (29) to compute the integraland therefore the divergence operator div : V → T is defined and can be computedexactly in the sense that div ( ν ) = π ∇ · ν for any ν ∈ V , where π denotes the L projection to T . The transpose of thediscrete divergence operator will give us a discrete approximation of the gradient.We can obtain an expression of the discrete gradient by reverting the order of thesum in (40). Let us denote by F ( η ) the set of faces to which the node η belongsand, for a face f k , we denote the neighboring cells of f k by E + k and E − k . From (40),we can rewrite (cid:90) Ω ˆ ψ ∇ · φ kη d x = − (cid:88) f ∈ F ( η ) ( ˆ ψ E + f − ˆ ψ E − f )( e k · n f ) (cid:90) f φ η d x . where the normal n f of the face f is directed from E − f to E + f . This conventionimplies that n jl = − n lj = n f if E j = E − f and E l = E + f . Hence, the discrete gradient operator grad is the mapping from scalar cell values to vector node valuegiven by(41) [ grad ( ˆ ψ )] η,k = (cid:88) f ∈ F ( η ) ( ˆ ψ E + f − ˆ ψ E − f )( e k · n f ) (cid:90) f φ η d x . Hence, gathering (38), (39) and (41), in this formulation, we obtain the followingexpression for (cid:98) f , as the discrete gradient of the discretized potential, that is (cid:98) f = grad ( ˆ ψ ) . In (41), the expression only depends on differences of the potential, which can beestimated locally without knowledge of the global potential, i.e.(42) ˆ ψ E + f − ˆ ψ E − f = ˆ f f · dr f , where ˆ f f is an approximation of the force on the face f and dr f is the vectorjoining the centroids of E − f and E + f . In practice, it means that the method canbe applied even if the force is not derived from a potential, as we can see that thepotential ψ does not have to be computed. Note that, in the numerical tests thatfollow, we have not tested this case.4.3. Interpretation of the discrete gradient approach using singular loadterm functions.
When we consider a cell-valued potential ψ , the correspondingforce f = ∇ ψ can be defined as a singular function with support on the cell faces.Let us define this class of function, which we will refer to as functions.Given an internal 2D surface S in Ω (or 1D line in 2D), we define the constant2D-Dirac function δ S ( x ) as the distribution given by < δ S , φ > = (cid:90) S φ ( x ) dx, for all φ ∈ C ∞ (Ω). The 2D-Dirac δ S is a measure which coincides with the Hausdorffmeasure on the d − S . Then, we can also define 2D-Dirac function h ( x ) δ S ( x ), for any h ∈ L ( S ) as < hδ S , φ > = (cid:82) S h ( x ) φ ( x ) dx, , for any φ ∈ C ∞ (Ω).If the surface S is Lipschitz, then a continuous trace operator from H (Ω) to H ( S )can be defined, see for example [8]. Therefore, at least if h ∈ L ( S ), we have that h ( x ) δ S ( x ) ∈ H − (Ω). Indeed, we have, for any φ ∈ H (Ω), < hδ S , φ > = (cid:90) S h ( x ) φ ( x ) dx ≤ (cid:107) h (cid:107) L ( S ) (cid:107) φ (cid:107) L ( S ) ≤ C (cid:107) h (cid:107) L ( S ) (cid:107) φ (cid:107) H ( S ) ≤ C C (cid:107) h (cid:107) L ( S ) (cid:107) φ (cid:107) H (Ω) , for two constants C and C . From this observation, we can infer that the originalsystem of equation (1) is well-posed for 2D Dirac vector functions f .Let us now consider a surface S that splits the domain Ω in two sub-domains,namely Ω − and Ω + , and a potential ψ which is piecewise constant and takes thevalues ψ ± in Ω ± . The gradient of ψ in the sense of distribution is defined as(43) < ∇ ψ, φ > = − (cid:90) ψ ∇ φ dx IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS13
Let us consider φ with compact support in Ω so that, after using integration bypart we obtain,(44) < ∇ ψ, φ > = − ( ψ − (cid:90) Ω − ∇ φ dx + ψ + (cid:90) Ω + ∇ φ dx ) == (cid:90) S (( ψ + − ψ − ) n ) φ dx, where n ( x ) denotes the normal to S at x ∈ S pointing from Ω − to Ω + . From thedefinition (43) and (44), we get that the gradient of ˆ ψ is a 2D Dirac vector functiongiven by(45) ∇ ψ = [ ψ + − ψ − ] n ( x ) δ S . Let us now consider again a cell-wise constant potential function ˆ ψ defined on amesh. Using the same notation as in the previous section, we infer from (45) thatthe gradient ˚ f of ˆ ψ in the sense of distribution is given by(46) ˚ f = ∇ ˆ ψ = (cid:88) f ∈ F int ( ˆ ψ E + f − ˆ ψ E − f ) n f δ f , where F int denotes the set of internal faces. Note that for the basis function φ kη asdefined in (8), we get(47) (cid:90) Ω ˚ f · φ kη = (cid:88) f ∈ F i ( ˆ ψ E + f − ˆ ψ E − f )( e k · n f ) (cid:90) f φ η ( x ) dx and we recover expression (41). Hence, the discrete gradient approach can beinterpreted in the following way. First, we approximate the volumetric load term f by a 2D Dirac function ˚ f with support on the cell faces and which is constanton each face, that is ˚ f has the form(48) ˚ f ( x ) = (cid:88) f ∈ F i ˚ f f δ f ( x )where ˚ f f is a constant vector, for each face f . In the case the force f is derivedfrom a potential, we can use the expression (46) to carry on this approximation.Otherwise, we propose to use expression (42) and consider˚ f f = ( ˆ f f · dr f ) n f . Once ˚ f is computed, we use the VEM method to solve the problem defined as ∇ · σ = ˚ f . Then, the assembly of the load term can be done exactly , as we can see from (47)in the case of a potential and otherwise(49) (cid:90) Ω ˚ f · φ kη = (cid:88) f ∈ F i ˚ f f · e k (cid:90) f φ η ( x ) dx. in the case where (48) is used. Note that the integrals in (47) and (48) can becomputed exactly we use the virtual basis proposed in [1]. Stability with respect to aspect ratio
Let us now discuss the choice of the stabilization matrix S in (32). In [5], theauthors propose S = I , which is the simplest choice. In [9], the authors look at several cell shapes andrecommend the stabilization term given by(50) S = α I where the constant α is chosen as(51) α G = | E | tr( ˆ C )tr( N Tc N c ) , as it gives an overall satisfactory approximation of the higher order nonlinear modes.This constant is stable with respect to isotropic scaling but it is not stable withrespect to the aspect ratio.5.1. Instability of α G with respect to aspect ratio. To demonstrate that, weconsider a rectangular element in 2D given by [ − h , h ] × [ − h , h ]. In this case,an explicit definition of the virtual element space is available, as it is spanned bythe four following functions(52) ϕ l ( x ) = 1 , ϕ l ( x ) = x h , ϕ l ( x ) = x h , ϕ ( x ) = x x h h , in each Cartesian direction. They coincide in this case to the standard finite ele-ments for quadrilaterals. The functions ϕ lj ( x ) e i for j = 1 , , i = 1 , P . Let ϕ i ( x ) = ϕ ( x ) e i . We check directly, using thesymmetry of the domain, that π ∇ ( ϕ i ) = 0 . Hence, for each basis functions in (52), we have that the zero and first order mo-ments correspond to those of their projection so that, indeed, they form a basis of V . Moreover { ϕ i } i =1 , constitutes a basis for ker π ∇ . In this two-dimensional case,the matrix N c is given by N ic = (cid:32) h √ h h √ h (cid:33) We collect the contributions of the four nodes of the cell and obtain the matrix N c given by(53) N Tc = h − h − h h h h − h − h √ h √ h √ h − √ h − √ h − √ h − √ h √ h which yields(54) N Tc N c = h h
00 0 2( h + h ) so that tr( N Tc N c ) = 6( h + h ). Hence, the scaling ratio α G is given by(55) α G = 4 h h tr( ˆ C )6( h + h ) = 23 tr( ˆ C )( ε + ε − ) , IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS15 where ε = h h denotes the aspect ratio. Let us now compute how this weight in thestabilization term compares with the actual energy for the functions that belong toker π ∇ . To do so, we consider an isotropic material where the stress is given as(56) σ = λ tr( ε ) + 2 µ ε , which implies(57) a ( u , ¯ u ) = (cid:90) Ω ( λ tr( ε ) tr( ¯ ε ) + 2 µ ε : ¯ ε ) d x For ϕ i we denote by ε i , the corresponding strain, which is given by(58) ε i = 12 ( e i ∇ φ T + ∇ φ e Ti ) . We get ε i : ε j = 12 (cid:18) δ i,j |∇ ϕ | + ∂ϕ∂x i ∂ϕ∂x j (cid:19) , where δ i,j = 1 if i = j and zero otherwise. Hence, using the symmetry of thedomain, we get (cid:90) E ε i : ε j d x = δ i,j (cid:90) E |∇ ϕ | d x . We have tr( ε i ) = ∂ϕ∂x i . Hence, using the symmetry of the domain we get (cid:90) E tr( ε i ) tr( ε j ) d x = δ i,j (cid:90) E (cid:12)(cid:12)(cid:12)(cid:12) ∂ϕ∂x i (cid:12)(cid:12)(cid:12)(cid:12) d x . Finally, the restriction of the bilinear form a to ker π ∇ takes the form a ( ϕ i , ϕ j ) = (cid:82) E ( λ (cid:12)(cid:12)(cid:12) ∂ϕ∂x (cid:12)(cid:12)(cid:12) + 2 µ |∇ ϕ | ) d x (cid:82) E ( λ (cid:12)(cid:12)(cid:12) ∂ϕ∂x (cid:12)(cid:12)(cid:12) + 2 µ |∇ ϕ | ) d x The integrals above can be computed exactly and we have (cid:90) E (cid:12)(cid:12)(cid:12)(cid:12) ∂ϕ∂x (cid:12)(cid:12)(cid:12)(cid:12) d x = 43 ε − , (cid:90) E (cid:12)(cid:12)(cid:12)(cid:12) ∂ϕ∂x (cid:12)(cid:12)(cid:12)(cid:12) d x = 43 ε, Hence, a ( ϕ i , ϕ j ) = (cid:18) λε − + µ ( ε + ε − ) 00 λε + µ ( ε + ε − ) (cid:19) . We denote by α and α the two eigenvalues of the matrix above. We obtainlim ε → , ∞ α α G = lim ε → , ∞ α α G = ∞ This enables us to conclude that, when the aspect ratio ε tends either to zero orinfinity, the ratios above tends to infinity so that we cannot find a constant c > ε , such that ca E ( u , u ) ≤ s E ( u , u ) , for all u ∈ ker π ∇ . It implies that the stabilization term is not stable with respectto the aspect ratio. An alternative choice of the stabilization scaling.
Instead of using α G ,let us use(59) α N = 19 | E | tr( ˆ C ) tr(inv( N Tc N c )) . Both α N and α G are invariant with respect to rotation. Because of the coefficient , we have that if N Tc N c were diagonal with constant coefficient, then α N and α G would be equal. But in general they differ and we have(60) α N = 2 λ + 6 µ (cid:18) ε + ε − + 2 ε + ε − (cid:19) It implies that lim ε → α α N = lim ε →∞ α α N = 163 λ + 2 µλ + 3 µ and lim ε →∞ α α N = lim ε → α α N = 163 µλ + 3 µ . Therefore, for this choice of α , there exist two constants c , c > independent of the aspect ration ε and such that c a E ( u , u ) ≤ s E ( u , u ) ≤ c a E ( u , u )for all u ∈ ker π ∇ . We can conclude that the stabilization provided by α N is stablewith respect to the aspect ratio, at least for quadrilaterals. Let us now try toexplain the motivation back the introduction of α N . We denote by λ i the singularvalues of N c and introduce the following averages λ arithm = ( d (cid:88) i =1 λ i ) and λ harm = ( d (cid:88) i =1 λ i ) , which, for simplicity, we refer to as arithmetic and harmonic averages. Note thatthe matrix N c , which is given in (26), accounts for the geometry and the unit ofeach coefficient is a unit length. We could therefore interpret the values of λ arihm and λ harm as characteristic lengths of the cell. Using these values, we can rewritethe scaling coefficients as α G = 1 λ | E | tr( ˆ C ) and α N = 1 λ | E | C ) , so that the difference between the two scalings is that they consider different typeof averages. Let us use eigenmodes to estimate the energy in each direction. Forsimplicity, we consider the Laplace equation and the normalized energy of the mode φ i ( x ) = cos( π h i x i ) in the i -th direction is given by (cid:82) K |∇ φ | dx (cid:82) K | φ | dx = π (2 h i ) , from which we infer that a typical scale for the energy in the direction x i is given by h i . If we consider a linear combination of such unidirectional functions and neglectthe interactions between them, then we are naturally led to consider the sum d (cid:88) i =1 h i IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS17 as a typical scale for the energy. To obtain a typical length, we end up by takingthe harmonic average as defined above.6.
Numerical test cases
The great advantage of VEM methods is that they are valid for very generalgrids including non-convex cells and more than one face between two cells, [3]. Thisproperty can be used to avoid curved faces on general cells, simply by triangulatingthe surface. The VEM theory does not cover curved surfaces and in the nextexamples we investigate the need for triangulation in 3D.6.1.
A two-dimensional compaction case. Case description : We considera rectangular domain made of an isotropic material with the following properties, ρ = 3 × kg m − , E = 3 × Pa and ν = 0 .
3. The vertical length of the grid is L y = 15m and the horizontal length will by determined by the aspect ratio L y /L x .Different values of the aspect ratios will be tested. The boundary conditions arezero displacement at the bottom, rolling boundary condition on the sides, that is nodisplacement in the normal direction and no force in the tangential direction. At thetop, we have no force and free displacement. Even if the model is two-dimensional,we have to set up boundary conditions for the third dimension, perpendicular tothe plane, as the material is going to expand or withdraw in this direction due tothe Poisson ratio. We impose zero displacement in the perpendicular direction, theother standard choice being no force in that direction. The load term is gravitation,that is, a constant vertical vector pointing downwards and we simulate the situationwhere the material is going to subside by the effect of its own weight, hence thename of compaction . An analytical solution is available for this case and given by(61) u = [0 , γ ( L y − ( y − L y ) )] γ = gρ C , , where C , is the second diagonal coefficient of the stiffness matrix C . We start witha Cartesian grid that we twist in order to avoid artifact effects from symmetries.We will refer to this grid as the twisted grid . We consider a variation of this gridwhere we add extra degrees of freedom in the form of extra nodes on the horizontaledges, see Figure 6.1. The motivation for introducing such extra nodes is explainedin the next paragraph. Results : We test the three different implementations of the load term, as de-scribed in Section 4. It is important to note that the first alternative, which usesthe projection operator, see the definition in (36), is exact in the case we are con-sidering. Indeed, since f = ρ g is constant, we have π i ( f ) = f so that no error is introduced by the assembly of the load term. In the remaining, wewill refer to this implementation of the load as the exact load term. In comparison,the third method is not exact, as the potential function, here given by ψ = ρgy , isapproximated by a cell-valued function. For the stabilization term, we test the twoscaling variables α G and α N presented in Section 5.We start with the scaling variable α G taken from [9] and the exact load imple-mentation. We use the grid with extra nodes. For such grid, each cell gets twoextra degrees of freedom. However, these extra degrees of freedom do not enrichthe approximation space as they do in the case of a finite element method. TheVEM method retains the same degree of accuracy, that is first order in our case. The extra basis functions introduced by the extra degrees of freedom are handledby the stabilization term. But the stabilization term only guarantees that theseextra functions do not break the ellipticity of the system but it is an artificial termwhich cannot add any accuracy. Therefore, by adding an extra node on the edges,we increase the relative importance of the stabilization term, so that its deficiencywill be more apparent. As predicted by the results of Section 5, we observe a se-vere dependence on the aspect ratio. When the aspect ratio is minimal, that is L x /L y = 1, then the solution is close to the analytical one but, when the aspectratio is increased to L x /L y = 10, by stretching the grid in the horizontal direction,the results deteriorate severely, see the top panels in Figure 3. We run the samesimulations but, instead of the exact load term, we use the load term computed bythe discrete gradient operator. Then, the results do not deteriorate as the aspectratio is increased.In Figure 4, we plot the error in displacement as a function of the aspect ratio(from 1 to 100) for the different grid cases and the three implementation of the loadterm. The left figure shows that the exact load and nodal load calculations fail forthe grid with extra nodes. The error apparently follows a second order growth,that is err ∼ ( L y /L x ) . The plot on the right shows the error for the twistedgrid without extra nodes for the exact load computation and the error for bothgrids for the discrete gradient approach. All the methods give reasonable results,but the exact load calculation seems to deteriorate more than the others. Thediscrete gradient approach is stable in both cases. Note that, if we had used a gridwithout disturbance, all the methods would give exact results for the grid withoutextra nodes on the faces while the extra node case will still fail for the exact loadcalculation. The reason is that, in the non disturbed case with no extra nodes, allthe implementations of the load term give the same result in the case of a constantvertical load term.Finally in Figure 5, we consider the scaling α N introduced in (59), which is stablewith respect to aspect ratio. The error does not grow as the aspect ratio is increased,as opposed to α G . The use of α N deteriorates the solution computed using thediscrete gradient approach, while it significantly improves the solution using theexact method. However, this conclusion is difficult to extend to more general cases.The value of α N has been derived from an analysis done on regular quadrilateralsand we observe that the stability properties extend to a twisted Cartesian grid.However, separate studies would have to be done for more complicated shapes andalso in 3D, where the situation is expected to be more complicated. Indeed, whilein 2D the aspect ratio is described by a scalar quantity namely ε = ∆ x ∆ y , in 3D weneed 2 values, say ∆ x ∆ z and ∆ y ∆ z , the third quantity ∆ x ∆ y being imposed by the fact thatwe will anyway require isotropic stability. It means that a scalar approximation ofthe stabilization term, as given in (50) and also in [5], will not be enough. Thisproblem was noticed in [2], and the exact stabilization term corresponding to finiteelement was used there to study a poro-elastic response function in the 3D case. Comment : We do not really understand why the discrete gradient approach(Alternative 2) performs significantly better than the projection approach which isexact in this case (Alternative 1). However, we note some fundamental differencesbetween the force-based methods (Alternatives 1 and 2) and the discrete gradientapproach, which may help to understand the differences in the results. As explainedin the previous section, see (34), the difference between the methods is in how they
IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS19 divide the weights between the nodes. All the force-based methods divide forcesaccording to a weight for each node associated with volume integrals. These weightsare equal for all
Cartesian directions. In contrast, the discrete gradient methoduses weights associated with surface integrals, so that the weights can depend onthe direction, and the corresponding degrees of freedom. These weights can beassociated with the projected area of the faces associated with a node divided bythe projection of the cell in the same direction. This is most easily seen from theexpression in equation (42). In the case of the extra nodes on the edges, these nodeswill have associated weights in the horizontal direction only due to the tilt of thegrid and the weights will in the simple case be doubled of the corner nodes whilethe exact case will give all nodes the same weights. In [9], the method using nodequadrature (Alternative 2) is considered, this will in the above case give a smallerweight to the midpoint and behave worse for the case with extra node, as seen inthe left panel of Figure 4.
Figure 2.
A twisted Cartesian grid is obtained by starting froma regular Cartesian grid and moving the nodes, here by using asmooth given displacement field. We plot the grid that is obtainedafter adding one extra node on each horizontal edge. Such gridis used to demonstrate the failure of the stabilization term wherethe aspect ratio is increased. The grid plotted here is the referencegrid with aspect ratio, by definition, equal to L x /L y = 1. L x /L y = 1 L x /L y = 10 E x a c t l oa d D i s c r e t e g r a d i e n t a pp r oa c h Figure 3.
We plot the computed displacement in the vertical di-rection for the 2D compaction example. The result for every nodeof the grid is represented as a dot, where the x -coordinate of thenode corresponds to the vertical position of the node and the y -coordinate corresponds to the value of the vertical displacementcomputed at the node. The analytical solution is plotted as a con-tinuous line. For these plots, the twisted Cartesian grid with anextra node on each horizontal edges, see Figure 6.1, has been used.The left column is for aspect ratio 1 and the right is for aspectratio 10. For the first row, the exact load calculation based on theexact integration of the VEM basis function has been used, whilethe lower row corresponds to the discrete gradient approach. Weuse the scaling factor α G as proposed in [9], see (51). We observethat, for the exact load computation, the solution quickly deteri-orates when the aspect ratio is increased while the results remaingood for the discrete gradient approach. IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS21
Aspect ratio L x /L y e rr ( d z ) / m a x ( d z ) exact forcenodal force Aspect ratio L x /L y e rr ( d z ) / m a x ( d z ) exact force cartesiandual div extra nodesdual div cartesian Figure 4.
Plots of the maximum error in the vertical displace-ment as a function of the grid aspect ratio. The left figure showsresults for the exact integration method in the case with extranodes on horizontal faces. We observe that the method fails asthe error blows up. The extra points in this plot are referencepoints that indicate a quadratic scaling of the error with respectaspect ratio. In the right figure, the results are shown for the ex-act method on the twisted grid without extra nodes and for thediscrete gradient approach on the same grid with and without theextra nodes.
Aspect ratio L x /L y e rr ( d z ) / m a x ( d z ) α G , exact α G , dual approach α N , exact α N , dual approach Figure 5.
Comparison between the two scaling constants α G and α N . We plot the error of the vertical displacement, as the aspectratios is increased. Here, we use the grid without the extra nodes.We observe that the scaling constant α N yields stability with re-spect to aspect ratio, independently of which method is used tocompute the load term. Compaction 3D.
In order to investigate the performance of the VEM methodon real reservoir geometries, we use two grids which includes standard features ofsubsurface models. The first one is based on a local sedimentary model called sbed . The model was used for upscaling permeability. Our version is 15 m × × m with logical Cartesian dimensions 15 × × p i,j and q i,j the bottom and top and the pillarthat is indexed by ( i, j ). For each region contained between the four pillars ( i, j ),( i + 1 .j ), ( i + 1 , j + 1) and ( i, j ), points are defined on each of this pillar in equalnumber. We denote those points by x ki (cid:48) ,j (cid:48) for i (cid:48) ∈ { i, i +1 } and j (cid:48) ∈ { j, j +1 } . Then,the region between the four pillars is meshed with hexahedrons with eight cornerpoints given such as x k (cid:48) i (cid:48) ,j (cid:48) for i (cid:48) ∈ { i, i + 1 } , j (cid:48) ∈ { j, j + 1 } and k (cid:48) ∈ { k, k + 1 } . Thisconstruction naturally leads to irregular cell shapes and faces that are not planar,see the illustrations given in Figure 6 . Therefore we end up outside the theoretical Figure 6.
On the right, two neighboring cells of a corner-pointgrid. On the left, examples of the irregular cell shapes that thecorner-point format can produce.framework of the VEM method, which only cover planar polygonal faces. However,the computation of the stiffness matrix for VEM relies on geometrical propertiesthat are all available, either as exact or approximated values (such as face areas,face normals, etc.), in the case of a corner-point grid, so that the stiffness matrixcan be assembled and a solution computed. To evaluate the error that is introducedby this geometrical approximation, we compare the solution obtained this way withthe solution that is obtained after triangulating the non planar surfaces, by addinga point in the middle of the faces. For such grid, the faces will be planar and thetheoretical framework of the VEM method applies.In Figure 7, we show the effect of compression with two types of load given by aconstant gravitational force and a constant load applied on the top surface. For both
IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS23 loads, the analytical solutions can be computed and they are respectively, quadraticand linear in z . We consider both the original corner-point grid and the triangulatedgrid. By triangulated grid, we mean a grid where the faces ares triangulated, aswe just explained. For all these cases, the VEM method gives accurate results,given that we use the discrete gradient approach to compute the load term. Theother alternatives simply fail in this case, by errors that are larger than the spanof the exact solution. In the sbed model, the pillars are all vertical lines, whichimplies that the vertical faces are planar. For the linear case corresponding to aconstant load on the top surface, we see that the triangulated version gives exactresult, as predicted by the VEM theory, since in this case all the surfaces areplanar and the solution is linear. For the original grid, we get an error due to thecurved top and bottom faces in each cell. For the pure gravitational case, bothgrids give comparable results. Thus, we can conclude that in practice, it may notbe worth triangulating the faces because it introduces more degrees of freedomwithout significantly improving the accuracy of the solution. We consider the caseof a flipped model for Norne in Figure 8. In this way, we can investigate the effectof having non planar surface in the vertical direction. Typically, for the cells of theoriginal reservoir, we have ∆ x ∆ z ≈ ∆ y ∆ z (cid:29) sbed model and show that effects of curvature on the faces can be neglected.This indicates that, for many practical applications, the VEM method can be useddirectly on the original grid of reservoirs without deteriorating the accuracy of theresults. Original corner-point grid Using triangulated faces
Figure 7.
Effect of compression on the sbed model. The firstrow shows a plot of the vertical deformation on the grid (left),the original grid (middle) and the same grid where the surfaceare triangulated (right). Two types of loads are considered: puregravitational compression (second row), load at the top surface(third row). The first column shows the displacement obtained foreach loading case, which is very close to the analytical solution.The remaining plots show the errors for the original cornerpointgrid with curved faces (middle column) and the triangulated gridwith only planar faces (right column).
IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS25
Original corner-point grid Using triangulated faces error at the extra nodes
Figure 8.
Effect of compression on a fliped sbed model. The firstrow shows a plot of the vertical deformation on the grid (left),the original grid (middle) and the same grid where the surface aretriangulated (right). We consider only the case with gravitationload. The first column shows the displacement. The remainingplots show the errors for the original corner-point grid (middlecolumn) and the triangulated grid (right column). On the plot atthe lower right, we observe that the error splits clearly between thetype of nodes, the extra face node at the bottom and the other atthe top.
Original pillars Vertical pillars
Figure 9.
Effect of compression on a part of the Norne model.The first column shows the plot the vertical deformation on the grid(left), the original grid after removing the padding (middle) andthe straightened grid where the pillars are made vertical (right).The figure in the upper left corner shows the bounding box whichis used for the calculation, while the two other grids show theembedded Norne grid. The second row show the results for puregravitational compression. The first column shows the verticaldisplacement while the second and third show the errors in thevertical displacement for the original and triangulated grid.
IRTUAL ELEMENT METHOD FOR GEOMECHANICAL SIMULATIONS OF RESERVOIR MODELS27 Conclusion
We have demonstrated how geomechanical calculations can be done directly oncomplex geological models frequently encountered in reservoir modeling, by usingthe flexibility of the VEM method which can handle general geometries. In thismethod, the energy is not computed exactly for each basis element functions. Wedemonstrate that this approximation can come at the cost of large errors for de-formed grids, if not care is taken when defining the approximate bilinear form. Inparticular we study the effect of the load term calculation and show that, with sta-bilization terms and load term calculations presented earlier in the literature, evensimple 2D cases fails severely when the aspect ratio is increased. We found that boththe choices of discretization and of the load term calculation are in combinationresponsible for the failure. Using the exact equivalence with FEM on quadrilateralgrid, we presented a modification of the discretization that makes the method morerobust in the 2D case. In addition, we demonstrated that a calculation of the loadin term of a gradient of a potential was robust in 2D and the only approach whichgave sufficient accuracy in 3D. This holds in particular for grid cells that are outsidethe reach of FEM, such as those containing hanging nodes. The VEM theory doesnot cover curved faces, which are common in subsurface models. We saw that forour tests the error associated with this feature was negligible comparable with othererrors, with the natural exception of the case when VEM gives the exact solution(linear displacement). acknowledgements
This work has been partially funded by the Research Council of Norway throughgrants no. 215641 from the CLIMIT programme.
References [1] B Ahmad, Ahmed Alsaedi, Franco Brezzi, L Donatella Marini, and A Russo. Equivalent pro-jectors for virtual element methods.
Computers & Mathematics with Applications , 66(3):376–391, 2013.[2] Odd Andersen, Halvor Møll Nilsen, and Sarah Gasda. Modelling geomechanical impact of co injection using precomputed response functions. In ECMOR XV – 15 th European Conferenceon the Mathematics of Oil Recovery, Amsterdam, Netherlands, 29 August - 1 September2016 . EAGE, 2016.[3] Louren¸co Beir˜ao da Veiga, F Brezzi, A Cangiani, G Manzini, LD Marini, and A Russo.Basic principles of virtual element methods.
Mathematical Models and Methods in AppliedSciences , 23(01):199–214, 2013.[4] Louren¸co Beir˜ao da Veiga, F Brezzi, LD Marini, and A Russo. The hitchhiker’s guide to thevirtual element method.
Mathematical models and methods in applied sciences , 24(08):1541–1573, 2014.[5] Louren¸co Beir˜ao da Veiga, Franco Brezzi, and L Donatella Marini. Virtual elements for linearelasticity problems.
SIAM Journal on Numerical Analysis , 51(2):794–812, 2013.[6] Louren¸co Beir˜ao da Veiga, Konstantin Lipnikov, and Gianmarco Manzini.
Mimetic FiniteDifference Method for Elliptic Problems , volume 11. Springer, 2014.[7] P. G. Bergan and M. K. Nyg˚ard. Finite elements with increased freedom in choosing shapefunctions.
Int. J. Numer. Meth. Engng. , 20(4):643–663, Apr 1984.[8] Zhonghai Ding. A proof of the trace theorem of sobolev spaces on lipschitz domains.
Pro-ceedings of the American Mathematical Society , 124(2):591–600, 1996.[9] Arun L Gain, Cameron Talischi, and Glaucio H Paulino. On the virtual element methodfor three-dimensional linear elasticity problems on arbitrary polyhedral meshes.
ComputerMethods in Applied Mechanics and Engineering , 282:132–160, 2014. [10] Emmanuel J. Gringarten, Guven Burc Arpat, Mohamed Aymen Haouesse, Anne Dutranois,Laurent Deny, Stanislas Jayr, Anne-Laure Tertois, Jean-Laurent Mallet, Andrea Bernal, andLong X. Nghiem. New grids for robust reservoir modeling.
SPE Annual Technical Conferenceand Exhibition , 2008.[11] KnutAndreas Lie, Stein Krogstad, Ingeborg Skjelkv˚ale Ligaarden, Jostein Roald Natvig,Halvor Nilsen, and B˚ard Skaflestad. Open-source MATLAB implementation of consistentdiscretisations on complex grids.
Comput. Geosci. , 16:297–322, 2012.[12] Bradley Mallison, Charles Sword, Thomas Viard, William Milliken, and Amy Cheng. Un-structured cut-cell grids for modeling complex reservoirs.
SPE Journal , 19(02):340–352, Apr2014.[13] The MATLAB Reservoir Simulation Toolbox, version 2016a, 7 2016.[14] Open Porous Media initiative. Open datasets, 2015. .[15] David K Ponting. Corner point geometry in reservoir simulation. In