Modified flux-vector-based Green element method for problems in steady-state anisotropic media

ARTICLE IN PRESS Engineering Analysis with Boundary Elements 33 (2009) 368–387

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Modified flux-vector-based Green element method for problems in steady-state anisotropic media P. Lorinczi a,, S.D. Harris b, L. Elliott c a b c

School of Earth and Environment, University of Leeds, Leeds LS2 9JT, UK Rock Deformation Research Ltd, School of Earth and Environment, University of Leeds, Leeds LS2 9JT, UK Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, UK

a r t i c l e in f o

a b s t r a c t

Article history: Received 15 August 2007 Accepted 16 June 2008 Available online 13 August 2008

In this study we present a new numerical technique for solving problems in steady-state heterogeneous anisotropic media, namely the ‘flux-vector-based’ Green element method (‘q-based’ GEM) for anisotropic media. This method, which is appropriate for problems where the permeability has either constant or continuous components over the whole domain, is based on the boundary element method (BEM) formulation for direct, steady-state flow problems in anisotropic porous media, which is applied to finite element method (FEM) meshes. For situations involving media discontinuities, an extension of this ‘q-based’ GEM formulation is proposed, namely the modified ‘q-based’ GEM for anisotropic media. Numerical results are presented for various physical problems that simulate flow in an anisotropic medium with diagonal layers of different permeabilities or around faults and wells, and they show that the new method, with the extensions proposed, is very suitable for steady-state problems in such media. & 2008 Elsevier Ltd. All rights reserved.

Keywords: Green element method Flux-vector-based GEM Anisotropy Heterogeneous porous medium Nodal flux condition Permeability tensor

1. Introduction The study of physical problems involving anisotropic media has received much interest in recent decades. However, analytically, the differential equations of such problems are difficult to solve. In this article we investigate the solution of steady-state problems in heterogeneous and anisotropic media in the presence of faults and wells using the flux-vector-based Green element method (‘q-based’ GEM). The concept of the ‘q-based’ GEM was introduced by Pecher et al. [1], in order to maintain the high-order accuracy of the GEM, diminished by the approximations used in the standard applications of the GEM, see Taigbenu [2]. This approach is based on the idea that the number of linearly independent vectors is the same as the number of dimensions of their vector space and retains the normal fluxes as unknowns. Taigbenu [3] used the flux-based GEM formulation, based on calculating directly the normal fluxes at the internal nodes, for solving heat conduction problems. The high level of accuracy of this formulation is attained with coarse grids. The numerical results obtained by the flux-based GEM formulation are compared with those obtained using the standard GEM of Taigbenu [2] and it

 Corresponding author. Tel.: +44 113 3436620.

E-mail address: [email protected] (P. Lorinczi). 0955-7997/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2008.06.004

is concluded that the solution provided by the flux-based GEM formulation is far more comprehensive. Popov and Power [4,5] presented a multi-domain decomposition integral equation method for domain dominant problems, for the linear and non-linear convective–diffusion equations. As in the GEM, in their approach the solution domain is divided into subdomains. However, in this approach the domain integrals are transformed into surface integrals at the contour of each subdomain via the dual reciprocity method (DRM), using interpolation functions. The system is solved for the internal field variables as well as for their derivatives, without any additional interpolation. The ‘q-based’ GEM was used by Lorinczi et al. [6], who extended it to the modified ‘q-based’ GEM. This latter method is applicable to heterogeneous isotropic media with rapid and highorder changes in the permeability of the media and therefore enables problems involving flow through layered sequences, across partially sealing faults and around wells to be treated. This paper is concerned with introducing the ‘q-based’ GEM for anisotropic media. Several examples, governed by the diffusion equation, with the components of the permeability tensor being continuous (either constant or variable) are solved using this method. For situations involving media discontinuities, an extension of this formulation is proposed, namely the modified ‘q-based’ GEM for anisotropic media, based, as for the isotropic case, on satisfying a nodal flux condition at each node and the

ARTICLE IN PRESS P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

continuity of the tangential pressure gradient across the elements sharing a node. This novel approach is applied to physical problems within a geological setting, namely simulating flow in an anisotropic medium with inclined layers of different permeabilities, and around faults and wells.

369

example, Chang et al. [7] or Garroni and Menaldi [8]) G2 ðx; x0 Þ ¼ 

jK ij j1=2 lnðRg Þ 2p

(5)

where jK ij j is the determinant of the inverse of the matrix K ij and the geodesic distance Rg is defined by d X

2. Flux-vector-based Green element method for anisotropic media

R2g ¼

2.1. Mathematical formulation

The points where Gd ðx; x0 Þ become singular, namely xi ¼ x0i, are called poles.

We consider an anisotropic medium in a bounded domain flux-vector-based (‘q-based’) GEM for anisotropic media is derived by applying the classical boundary element method (BEM) for steady-state anisotropic media to elements of a finite element method (FEM) mesh. In anisotropic media, provided that laminar flow conditions prevail, that is at low Reynolds numbers, Darcy’s law can be formulated as K

m

rp

d X K ij qp m qxj j¼1

2.3. Integral equations The governing partial differential equation (4) is transformed into an equivalent integral equation by means of the fundamental Green’s function and Green’s second identity, namely  Z  q q lðxÞpðxÞ ¼ Gd ðx; x0 Þ pðx0 Þ  pðx0 Þ Gd ðx; x0 Þ dGx0 þ þ qn qn G Z 0 0 þ Gd ðx; x ÞFðx Þ dLx0 (7) L

(1)

where q is the flux or velocity vector, K is the heterogeneous and anisotropic permeability tensor, m is the fluid viscosity and p is the fluid pressure. The components of the flux q are then written as follows: qi ¼ 

for i ¼ 1; d

(2)

where K ij ðxÞ and qi represent the components of the tensor KðxÞ and the vector q, respectively, and xj denote the Cartesian coordinates. The components of the tensor can be continuous functions of the space variables. The governing equation of steady-state, single-phase flow in a heterogeneous and anisotropic porous medium can be derived from Darcy’s law and the conservation of mass to have the following form:

where lðxÞ is the corner coefficient at x, namely lðxÞ ¼ 1 if x 2 L, lðxÞ ¼ 12 if x 2 G and G is straight at x, lðxÞ depends on the angle at the corner points of L, see Symm and Pitfield [10] or Brebbia et al. P [9], and q=qnþ ¼ di;j¼1 K ij cosðn; xi Þq=qxj , with cosðn; xi Þ the direction cosines of the normal n to the surface G. 2.4. Numerical discretisation of the integral equations Some form of numerical approximation is necessary in order to obtain a suitable reduction of the boundary integral equation (7) to an algebraic form that can be solved by a numerical approach, as this integral equation may rarely be solved analytically. After the discretisation of the domain imposed by the GEM, Eq. (7) is discretised and it becomes

lðeÞ ðxÞpðeÞ ðxÞ ¼

Ne Z X e¼1

r  ðKðxÞrpÞ ¼ FðxÞ in L

i;j¼1

K ij

q2 p ðxÞ ¼ FðxÞ qxi qxj

(4)

The tensor KðxÞ is assumed to be symmetric and positive-definite so that Eq. (4) is of the elliptic type. Different boundary conditions may be imposed on the boundary qL ¼ G, namely Dirichlet, Neumann or Robin conditions.

The 2D fundamental Green’s function for the steady-state diffusion equation in anisotropic media is given by (see, for

 Gd ðx; x0 Þðq  nÞðx0 Þ

(8)

LðeÞ

in which the superscript ðeÞ denotes a typical element from the domain discretisation, pðeÞ ðxÞ represents the value of p at the source point x in the element ðeÞ; Ne is the total number of polygonal elements that are used in discretising the domain, and GðeÞ and LðeÞ are, respectively, the boundary and domain of the typical element. The same idea as that from the ‘q-based’ Green element method presented in Pecher et al. [1] is used, namely that the maximum number of linearly independent vectors is the same as the dimension of their vector space. The normal component of the flux qn is replaced by the components qx ; qy (and qz in three dimensions), where qn¼

d X

qi ni

i¼1

¼ 2.2. Fundamental Green’s function

GðeÞ



 q Gd ðx; x0 Þ dGðeÞ pðx0 Þ x0 þ qn  Z þ Gd ðx; x0 ÞFðxÞ dLðeÞ 0 x

(3)

where FðxÞ is an internal/external source forcing function that incorporates the fluid viscosity. The GEM requires that the computational domain is discretised by suitable polygonal elements that collectively, represent the shape of the domain. We will assume that the physical properties of the medium (that is, the components K ij ; i; j ¼ 1; d) are constant within each element of the discretisation, but that they can change from element to element. In this case, for each element ðeÞ we have ðeÞ ðeÞ K ij  K¯ ij , where K¯ ij is some mean value of K ij within that element, i; j ¼ 1; d. This approximation essentially treats the entire domain as a piecewise homogeneous medium. Then the governing equations reduce to several partial differential equations of the following form: d X

(6)

i;j¼1

L  Rd , where d is the dimensionality of the space domain. The

q¼

K ij ðxi  x0i Þðxj  x0j Þ

8 < qx nx þ qy ny

in two dimensions

: qx nx þ qy ny þ qz nz

in three dimensions

(9)

with n the unit outward vector normal to the boundary. We will further refer to a rectangular, non-uniform grid; thus in this case each element has four nodes (corners).

ARTICLE IN PRESS 370

P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

At this stage, as in the isotropic case, the pressure p, as well as the components qx , qy (and qz in three dimensions) of the normal flux and F, is approximated by the linear interpolation functions Oj ; j ¼ 1; 4, of the Lagrange family. For example, qx ¼

4 X

Oj qx;j

(10)

Lij ¼ ½Ajþ1 ðxi Þ  C jþ1 ðxi Þ  C j ðxi Þ

for i; j ¼ 1; 4

(21)

where Aj ðxÞ ¼ Aðx; xj1 ; xj Þ

(22)

Bj ðxÞ ¼ Bðx; xj1 ; xj Þ

(23)

C j ðxÞ ¼ Cðx; xj1 ; xj Þ

(24)

Dj ðxÞ ¼ Dðx; xj1 ; xj Þ

(25)

j¼1

where we denote by qx;j the value of qx at node j. In Eq. (8) we consider x ¼ xi (where i ¼ 1; 4) and x0 ¼ x to denote a source and a field node, respectively, related to the 2D case, namely d ¼ 2. Then, by using the approximations mentioned above, Eq. (8) is applied to each node from each element and the system obtained can be written more compactly as Ne X

ðRðeÞ p  LðeÞ q  T ðeÞ FjÞ ¼ 0 ij j ij j ij

(11)

e¼1

where the element matrices Rij ; Lij and T ij can be expressed as Z q ðeÞ G2 ðx; xi ÞOj dGðeÞ þ dij l (12) RðeÞ ij ¼ ðeÞ qnþ G LðeÞ ¼ ij T ðeÞ ¼ ij

Z GðeÞ

Z LðeÞ

G2 ðx; xi ÞOj dGðeÞ

(13)

G2 ðx; xi ÞOj dLðeÞ

(14)

The expressions for these element matrices are presented more explicitly in the next section. 2.5. Evaluation of integrals for a rectangular grid 2.5.1. Boundary integrals The evaluation of the element matrices Rij and Lij for rectangular elements involves the integration of combinations of powers of the variable z (the natural coordinate along a boundary element line, which increases from zero to unity) and either G2 ðx; xi Þ or qG2 ðx; xi Þ=qnþ. Thus their calculation reduces to the evaluation of integrals of the following type: Z Aðx; x01 ; x02 Þ ¼ G2 ðx; x0 Þ dGx0 for x ¼ ðx; yÞ 2 L (15) G

Bðx; x01 ; x02 Þ ¼ Cðx; x01 ; x02 Þ ¼ Dðx; x01 ; x02 Þ ¼

Z G

Z

qG2 ðx; x0 Þ dGx0 qnþ G2 ðx; x0 Þz dGx0

for x ¼ ðx; yÞ 2 L

for x ¼ ðx; yÞ 2 L

(16)

(17)

G

Z G

qG2 ðx; x0 Þz dGx0 qnþ

for x ¼ ðx; yÞ 2 L

(18)

where G2 ðx; x0 Þ ¼

jK ij j1=2 ln 2p

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K 11 ðx  x0 Þ2 þ 2K 12 ðx  x0 Þðy  y0 Þ þ K 22 ðy  y0 Þ2

(19)

(thus each of the integrals Aj ðxÞ; Bj ðxÞ; C j ðxÞ and Dj ðxÞ calculates the given integral on the line having the end-points xj1 ; xj , for the source node x, with the integration being performed in an anticlockwise direction). The value of the coefficients on the diagonal can easily be calculated once the off-diagonal coefficients are all known, i.e. (see Brebbia et al. [9]) X Rii ¼  Rij for i ¼ 1; 4 (26) jai

The element matrix Lij is divided, as in the isotropic case, into two parts: Lbij and Lfij to account, respectively, for the flux before and following the node, in an anticlockwise direction of traversing the sides of the element and these are defined as follows: Lbij ¼ C j ðxi Þ

for i; j ¼ 1; 4

Lfij ¼ ½Ajþ1 ðxi Þ  C jþ1 ðxi Þ

(27) for i; j ¼ 1; 4

2.5.2. Domain integrals The domain integrals that are evaluated depend on the interpolation functions, the type of elements used in the discretisation of the domain and the position of the source node relative to the other nodes of the element. When rectangular elements and polynomial functions of the Lagrange family for the four-node rectangular element are used, the domain integrals that have to be evaluated are of the form   Z Z jK ij j Dy Dx  x m y n IL;mn ¼  2p 0 Dx Dy 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 11 2 12  ln K x þ 2K xy þ K 22 y2 dx dy (29) where Dx and Dy are the dimensions of an element from the domain discretisation. These integrals have been analytically evaluated for m; n ¼ 0; 1 and their expressions are presented in Lorinczi [12]. There are three unknowns at each internal node, namely the pressure and the two components of the normal flux. Although the original ‘q-based’ GEM theory does not require it (see Pecher et al. [1]), for every internal node we use the maximum number of equations, namely four, for generating the system to be solved. Therefore, we will obtain an over-specified system of equations, which can be solved by the conjugate gradient method (CGM), see Shewchuk [13]. The CGM is effective for systems of the form AX ¼ B

can be found, see Shewchuk [13], from

Rij ¼ Bjþ1 ðxi Þ  Djþ1 ðxi Þ þ Dj ðxi Þ

A> AX ¼ A> B

(20)

(30)

where X is an unknown vector, B is a known vector and A is a known, square, symmetric, positive-definite (or positive-indefinite) matrix. However, the CGM can be used to solve systems where A is not symmetric, not positive-definite and even not square. A solution to the least squares problem

K ij are the components of the inverse of matrix K and G is the straight line having the end-points x01 ¼ ðx01 ; y01 Þ; x02 ¼ ðx02 ; y02 Þ. The integrals A; B; C and D have been evaluated analytically for every x 2 R2 and any straight line G in the BEM using linear boundary elements, see Mera [11]. If we denote by xj ¼ ðxj ; yj Þ; j ¼ 1; 4, the four corners of a rectangular element, with the convention that x0 ¼ x4 and x5 ¼ x1 , then matrices Rij and Lij can be expressed as for i; j ¼ 1; 4; iaj

(28)

min kAX  Bk2 X

(31)

(32)

ARTICLE IN PRESS P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

A> A is symmetric and positive (for any X; X > A> AX ¼ kAXk2 X0Þ. If A is square and non-singular, the solution to Eq. (32) is the solution to Eq. (30). If A is not square, and the system is overconstrained (has more linearly independent equations than variables) then there may be or may not be a solution to (30), but it is always possible to find a value of X that minimises expression (31), the sum of the squares of the errors of each linear equation.

371

8

7

Kij,3

Kij,4 (3)

(4)

4

3.1. Internal source nodes To illustrate this approach, we refer to a rectangular grid and consider first the internal source node 5 in Fig. 1. The permeability tensors corresponding to each of the four elements in this figure are denoted by K ij;m ; i; j ¼ 1; 2, where m ¼ 1; 4 labels the element number. We consider first separate flux vectors within each of the elements (1)–(4). As in the isotropic case, we denote the x components of the flux LB at the source node 5 by qLT x and qx for the left-top and the leftbottom components of qx in the positive x direction, and by qRT x and qRB x for the right-top and the right-bottom components of qx in the negative x direction, respectively. Similarly, for the y RB components of the flux, we denote by qLB y and qy the left-bottom and the right-bottom components of qy in the positive y direction, RT and by qLT y and qy the left-top and the right-top components of qy in the negative y direction, respectively. The pressure gradient must be continuous across the four interfaces around node 5 as follows. The tangential pressure gradient qp=qx is continuous across the interfaces between elements (1) and (3), and (2) and (4), while qp=qy is continuous across the interfaces between elements (1) and (2), and (3) and (4). The x components of the flux can be defined with the notations introduced as  ð1;3Þ  ð1;2Þ qp qp qLB  K 12;1 x ¼  K 11;1 qx qy

6

5

3. Modified ‘q-based’ GEM for anisotropic media—development for a rectangular grid The ‘q-based’ GEM for heterogeneous and anisotropic media described in Section 2 is suitable for anisotropic media in which the permeability is either the same over the whole domain, with the tensor K ij having constant components, or where the permeability is such that the tensor K ij has continuous components. In the latter case an average will be used for every element, leading to only small differences in the tensor components between the elements. However, if there is a discontinuity in the permeability of the domain, then the simple ‘q-based’ GEM for anisotropic media is no longer suitable. This is particularly evident when large differences in the permeability are prescribed across the permeability discontinuities. The different values of the permeability that can apply for the elements sharing a node lead to an inconsistent definition of the tangential flux across the associated element boundaries, as the definition of the flux incorporates the value of the permeability. In this section we present an extension of the ‘q-based’ GEM for anisotropic media that is aimed at overcoming the difficulties encountered when a discontinuity in the permeability of the media occurs, namely the modified ‘q-based’ GEM for anisotropic media. This approach is based on the same idea as the modified ‘q-based’ GEM for isotropic media (see Lorinczi et al. [6]), namely satisfying a nodal flux condition at each node and the continuity of the tangential pressure gradient across the elements sharing a node.

9

(1)

(2)

Kij,1

Kij,2

1

2

3

Fig. 1. The permeabilities of the elements sharing an internal node.

qRB x ¼  K 11;2 qLT x ¼  K 11;3 qRT x ¼  K 11;4



ð2;4Þ



ð1;3Þ



ð2;4Þ

qp qx

 K 12;2

qp qx

 K 12;3

qp qx

 K 12;4



ð1;2Þ



ð3;4Þ



ð3;4Þ

qp qy qp qy qp qy

(33)

where the superscript of the tangential pressure gradient denotes the two elements that share the interface on which this derivative applies. For example, the notation ðqp=qxÞð1;3Þ indicates that this tangential pressure gradient applies at the interface between elements (1) and (3), but the notation also ensures the continuity of this tangential pressure gradient. Furthermore, this value applies for both elements (1) and (3) in the immediate vicinity of node 5. Similarly, the y components of the flux can be defined as qLB y ¼  K 12;1 qRB y ¼  K 12;2 qLT y ¼  K 12;3 qRT y ¼  K 12;4



ð1;3Þ



ð2;4Þ



ð1;3Þ



ð2;4Þ

qp qx

 K 22;1

qp qx

 K 22;2

qp qx

 K 22;3

qp qx

 K 22;4



ð1;2Þ



ð1;2Þ



ð3;4Þ



ð3;4Þ

qp qy qp qy qp qy qp qy

(34)

However, in the ‘q-based’ GEM for anisotropic media, the components of the flux within an element and in the vicinity of a node depend not only on the tangential pressure derivative coming from just one edge of the element but also on the tangential pressure derivatives coming from both edges of the element, namely on both qp=qx and qp=qy. Therefore, the continuity of the tangential pressure gradient across the elements sharing a node cannot lead in this case to a simple relation between two components of the flux coming from neighbouring LT LB RB elements, as we had between qLB x and qx or between qy and qy , for example, in the isotropic case, see Lorinczi et al. [6]. Thus, in this case, rather than expressing our equations in terms of the pressure and two flux components, we will instead consider two components of the pressure gradient from one element, namely ðqp=qxÞð1;3Þ and ðqp=qyÞð1;2Þ , as the principal unknowns at node 5, and all the components of the flux at node 5 will be expressed in terms of these unknowns. So in this approach the unknowns at each node are the pressure and two pressure gradient terms.

ARTICLE IN PRESS 372

P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

In order to express all the flux components at the source node 5 in terms of ðqp=qxÞð1;3Þ and ðqp=qyÞð1;2Þ , we use the nodal flux condition at this node. To illustrate this condition, we consider a small square of side 2 centred on node 5 and with the interfaces between the elements sharing node 5 being perpendicular to the sides of the square. Thus the nodal flux conditions can be written as follows: LT RT RB qLB x þ qx þ qx þ qx ¼ 0 LB LT RT qy þ qy þ qy þ qRB y ¼ 0

(35)

in the x and y directions, respectively, and we can take the limit as  ! 0. Introducing expressions (33) and (34) into Eqs. (35), after the limit as  ! 0 has been taken, together with the continuity of the pressure gradient conditions, results in a system that can be solved for ðqp=qxÞð2;4Þ and ðqp=qyÞð3;4Þ , thus enabling these tangential pressure components to be expressed in terms of ðqp=qxÞð1;3Þ and ðqp=qyÞð1;2Þ as follows:  ð2;4Þ  ð1;3Þ  ð1;2Þ qp qp qp ¼ c1;x þ c1;y qx qx qy  ð3;4Þ  ð1;3Þ  ð1;2Þ qp qp qp ¼ c2;x þ c2;y (36) qy qx qy where c1;x ¼ ððK 12;4  K 12;3 ÞðK 12;3  K 12;1 Þ þ ðK 22;4 þ K 22;3 ÞðK 11;1 þ K 11;3 ÞÞ=D

(37)

c1;y ¼ ððK 12;3  K 12;4 ÞðK 22;1 þ K 22;2 Þ þ ðK 22;4 þ K 22;3 ÞðK 12;1  K 12;2 ÞÞ=D

(38)

(39)

c2;y ¼ ððK 12;2  K 12;4 ÞðK 12;1  K 12;2 Þ þ ðK 11;2 þ K 11;4 ÞðK 22;1 þ K 22;2 ÞÞ=D

(40)

D ¼ ðK 12;2  K 12;4 ÞðK 12;4  K 12;3 Þ þ ðK 11;2 þ K 11;4 ÞðK 22;4 þ K 22;3 Þ

(41)

In this way, for an internal node, the definitions of the components of the pressure gradient (and also of the flux components) will incorporate the components of the permeability tensors from all four elements sharing that node. Expressions (36) are introduced into the equations obtained by integrating in turn over each element. In this way, the unknowns at each node for which the global system is solved are the pressure and the two components of the pressure gradient denoted by ðqp=qxÞð1;3Þ and ðqp=qyÞð1;2Þ , which are determined directly from the system of equations obtained. Once the two components of the pressure gradient have been determined, we use the nodal flux conditions, given by (35), to express the values of the two components of flux. These conditions can be rewritten as LT RT RB qLB x þ qx ¼ ðqx þ qx Þ LT RT RB qLB y þ qy ¼ ðqy þ qy Þ

(42)

We introduce the following notations: LT qLB x þ qx ; 2 LB RB qy þ qy qBy ¼ ; 2

qLx ¼

qRx ¼ qTy ¼

RT qRB x þ qx 2 RT qLT y þ qy

2

  ð1;3Þ 1 K 12;3 c2;x qp qLx ¼  qRx ¼  ðK 11;1 þ K 11;3 Þ  2 2 qx  ð1;2Þ 1 qp  ðK 12;1 þ K 12;3 c2;y Þ 2 qy  ð1;3Þ 1 qp ðK 12;1 þ K 12;2 c1;x Þ 2 qx  ð1;2Þ 1 qp  ðK 22;1 þ K 22;2 þ K 12;2 c1;y Þ 2 qy

(43)

with qLx and qRx oriented in the positive and in the negative x directions, respectively, and qBy and qTy oriented in the positive and

(44)

qBy ¼ qTy ¼ 

(45)

where c1;x ; c1;y ; c2;x and c2;y are given by Eqs. (37)–(40). 3.2. Side source nodes We refer first to the side source nodes on the boundaries parallel to the Ox-axis, namely to node 2 in Fig. 1. With the same notation as for an internal node, for side source node 2 there will RT be four components of the flux, namely qLT x and qx (oriented in the positive and the negative x directions, respectively), and qLT y and qRT y (both oriented in the negative y direction), defined as follows: qLT x ¼  K 11;1 qRT x ¼  K 11;2

c2;x ¼ ððK 12;2  K 12;4 ÞðK 11;1 þ K 11;3 Þ þ ðK 11;2 þ K 11;4 ÞðK 12;1  K 12;3 ÞÞ=D

in the negative y directions, respectively, leading to qLx ¼ qRx and to qBy ¼ qTy . To derive expressions for these flux components, all the quantities are written in terms of ðqp=qxÞð1;3Þ and ðqp=qyÞð1;2Þ . The newly introduced components of qx and qy will be given by

qLT y ¼  K 12;1 qRT y ¼  K 12;2



ð1Þ



ð2Þ



ð1Þ



ð2Þ

qp qx qp qx qp qx qp qx

 K 12;1  K 12;2  K 22;1  K 22;2



ð1;2Þ



ð1;2Þ



ð1;2Þ



ð1;2Þ

qp qy qp qy qp qy qp qy

(46)

where the superscripts (1) and (2) refer to the tangential pressure derivatives at the boundary interfaces of the corresponding elements. The outward normal flux to the boundary, denoted by qn , is assumed to apply at the side source node 2 for both elements (1) and (2) and is either an unknown (when pressure is known at this node) or it is given (when pressure is unknown at node 2). In expressing the flux components in relations (46) there are three components of the tangential pressure gradient that are used. These, together with the pressure and the outward normal flux to the boundary, result in five unknowns at the side source node 2. The equations necessary to find these are the nodal flux conditions, the two integral equations and the boundary condition at the side source node 2. 3.2.1. Prescribed pressure boundary condition at the side source node We first refer to the situation when the pressure is known on the boundary at the side source node 2. We consider ðqp=qxÞð1Þ and ðqp=qyÞð1;2Þ as the principal tangential pressure gradient unknowns at the side source node 2. Therefore we express ðqp=qxÞð2Þ in terms of these using the nodal flux condition. To illustrate this condition, we consider a small square of side 2 centred on node 2. The nodal flux condition in the x direction can be written as follows (after taking the limit as  ! 0): RT qLT x ¼ qx

(47)

ARTICLE IN PRESS P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387 RT Using the expressions for qLT x and qx from Eqs. (46) in condition (47) allows us to find an expression for ðqp=qxÞð2Þ as follows:  ð2Þ  ð1Þ  ð1;2Þ qp K 11;1 qp K 12;1 þ K 12;2 qp ¼  (48) qx K 11;2 qx K 11;2 qy

which is used in the equations arising from elements (1) and (2). In this way, ðqp=qxÞð1Þ and ðqp=qyÞð1;2Þ are determined when the system of equations is solved. Their values are then used to find ðqp=qxÞð2Þ at the side source node 2 using (48). In the y direction, the nodal flux condition is (after taking the limit as  ! 0) RT 2qn ¼ qLT y þ qy

(49)

RT Using the definitions of qLT y and qy given by (46) together with Eq. (48) in condition (49), we obtain the expression for qn as follows:   ð1Þ 1 K 11;1 qp K 12;1 þ K 12;2 qn ¼ 2 K 11;2 qx   ð1;2Þ 1 K 12;1 þ K 12;2 qp K 22;1  K 22;2 þ K 12;2 þ (50) 2 K 11;2 qy

Thus, when using the boundary condition at the node 2, all five unknowns at this node are completely determined. 3.2.2. Prescribed flux boundary condition at the side source node If the flux is prescribed on the boundary at the side source node 2, then, as in the previous case, the nodal flux conditions in the x and y directions lead to Eqs. (48) and (50). However, now condition (50) is used to express ðqp=qyÞð1;2Þ in terms of ðqp=qxÞð1Þ as follows:   ð1Þ 1 K 11;1 qp  ð1;2Þ K 12;1 þ K 12;2 qn  qp 2 K 11;2 qx  (51) ¼  1 K 12;1 þ K 12;2 qy K 22;1  K 22;2 þ K 12;2 2 K 11;2 In this way, at the side source node 2, only one of the unknown pressure tangential derivatives, namely ðqp=qxÞð1Þ, will enter into the integral equations and at the same time the equations will include the value of the flux. Once the pressure and ðqp=qxÞð1Þ have been determined by solving the system of equations, we can find ðqp=qyÞð1;2Þ using relation (51). Likewise, ðqp=qxÞð2Þ can be determined from Eq. (48) and thus all five unknowns at this node are completely determined. For a node situated on the boundaries parallel to the Oy-axis, namely to the side source node 4 in Fig. 1, a similar treatment can be applied. In all the above situations, once that the unknowns are determined at each side source node, the values of the unknown components of the flux can be directly calculated. 3.3. Corner source nodes At a corner source node, such as node 1 in Fig. 1, there are two flux components, corresponding to each of the two interfaces RT adjacent to node 1, namely (using the earlier notation) qRT x and qy defined as  ð1Þ  ð2Þ qp qp qRT  K 12;1 x ¼  K 11;1 qx qy  ð1Þ  ð2Þ q p q p qRT  K 22;1 (52) y ¼  K 12;1 qx qy where the superscripts (1) and (2) of the pressure derivative denote the edges of the element between nodes 1 and 4, and 1 and 2, respectively. The outward normal flux components that are considered at node 1 are denoted by qn1 and qn2 (normal to the

373

edges joining nodes 1 and 4, and 1 and 2, respectively, and oriented in the negative x and y directions, respectively). There are thus five unknowns at the corner source node 1, namely the pressure, ðqp=qxÞð1Þ ; ðqp=qyÞð2Þ ; qn1 and qn2 . The conditions used to determine these unknowns are the integral equation arising from element (1), the nodal flux condition and the boundary conditions at the corner node 1. For a corner source node in a rectangular grid we consider three different boundary condition formulations. 3.3.1. Both flux components at the corner source node If both flux components are specified at the corner source node, then the nodal flux condition is applied by considering a small square of side 2, with  ! 0, centred on node 1. In the x and y directions, the nodal flux conditions take the simple form qRT x ¼ qn1 ;

qRT y ¼ qn2

(53)

and when combined with Eqs. (52) these provide two relationships between qn1 ; qn2 and the two unknown pressure derivatives ðqp=qxÞð1Þ and ðqp=qyÞð2Þ . By using these conditions together with the two boundary conditions and the integral equation generated on element (1), we can determine all five unknowns at the corner source node 1. 3.3.2. Pressure and one flux component at the corner source node If the pressure and one flux component are specified at the source node 1, then the nodal flux conditions in the x and y directions are again given by (53). However, in this case, after we combine these equations with (52), which yields qn1 ¼ K 11;1

qn2 ¼ K 12;1





qp qx qp qx

ð1Þ

ð1Þ

 K 12;1

 K 22;1





qp qy qp qy

ð2Þ

ð2Þ

(54a)

(54b)

we must express one pressure derivative component in terms of the other. If the flux is known on the boundary parallel to the Oy-axis, then we use (54a) to write "  ð1Þ  ð2Þ # qp 1 qp (55) ¼ qn1  K 12;1 K 11;1 qx qy and, if the flux is known on the boundary parallel to the Ox-axis, we will use (54b) to write "  ð2Þ  ð1Þ # qp 1 qp ¼ qn2  K 12;1 (56) K 22;1 qy qx Therefore, if the flux is known on the boundary parallel to the Oy-axis, for example, then ðqp=qyÞð2Þ is determined directly by solving the system. We then find ðqp=qxÞð1Þ using Eq. (54a) and qn2 using (54b). Finally, by the use of the boundary conditions at the node 1, all the unknowns are completely determined. 3.3.3. Pressure only at the corner source node If, however, only the pressure is specified at the corner source node 1, then we make the assumption that the pressure derivative components ðqp=qxÞð1Þ and ðqp=qyÞð2Þ are equal in magnitude (leading to flux components equal in magnitude), as insufficient conditions are available to determine the two pressure derivative components. Thus the number of unknowns at the node 1 reduces to four in this case. These are determined using conditions (53), the integral equation on element (1) and the pressure boundary condition.

ARTICLE IN PRESS 374

P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

4.1.1. Example 1: homogeneous and anisotropic permeability In this example we consider the differential operator

3.4. Principal directions of the permeability tensor In order to have a better physical representation of some of the problems, we provide the principal values of the permeability tensor, namely K 1 and K 2 (also called the maximum and minimum principal values, respectively), and the angle of rotation y from the horizontal x-axis to the direction of the maximum principal value K 1 (measured in an anticlockwise direction). The coefficients of the permeability tensor with reference to the ðx; yÞ coordinate system are given by 9 > K 11 ¼ K 1 cos2 y þ K 2 sin2 y > = K 12 ¼ K 21 ¼ ðK 1  K 2 Þ cos y sin y (57) > > ; K ¼ K cos2 y þ K sin2 y 22

2

1

In all the problems where the permeability is given by the principal directions of the permeability tensor, the principal directions are illustrated in the figures of the numerical results.

3.5. Well test simulation Some problems simulating two wells in an anisotropic medium will be investigated using the modified ‘q-based’ GEM for anisotropic media. For these examples the internal/external source term from Eq. (3) is given in a similar way as to the isotropic case, namely by FðxÞ ¼ 

Nw X

qj dðx  xj Þdðy  yj Þ

on L

(58)

j¼1

where qj is the given strength of the jth well, Nw represents the number of wells in the domain, d is the Dirac delta function and ðxj ; yj Þ; j ¼ 1; Nw , are the coordinates of the points where the wells are located. Domain integrals for the examples involving wells have been avoided by specifying the two components of the flux at each corner of the cells containing the wells, namely jqx;j j ¼ jqy;j j ¼

jqj j ; 4ðDxÞj

j ¼ 1; N w

L½pðx; yÞ ¼

q2 p q2 p q2 p þ2 2  qx2 qxqy qy

(61)

and define the diffusion equation in an anisotropic medium with a non-zero external/internal forcing term given by

L½pðx; yÞ ¼ 2

(62)

when the components of the permeability tensor K are K 11 ¼ 1; K 12 ¼ K 21 ¼ 0:5 and K 22 ¼ 2 on the whole of the domain. The analytical pressure distribution is given by p ¼ x2 þ 2y

(63)

when the following boundary conditions are imposed: pðx; 0Þ ¼ x2 ;

pðx; 1Þ ¼ x2 þ 2;

qx ð0; yÞ ¼ 1;

qx ð1; yÞ ¼ 1;

0pxp1 0pyp1

(64)

The numerical results have been investigated for different mesh sizes. A comparison between the numerical and the analytical solutions for the pressure is presented in Fig. 2 for 128  128 elements. As expected from the discussion on the classical GEM inaccuracy, presented in Lorinczi et al. [6], the use of linear interpolation functions in problems with quadratic solutions does not provide an exact numerical solution. Improved numerical solutions could have been obtained by using higher-order elements, with quadratic functional representation for the unknowns. However, good agreement between the two sets of solutions for this example can still be observed with linear functions. Tables 1 and 2 present the numerical and the analytical solutions for the pressure and the two components of the flux at four points, namely ð0:25; 0:25Þ; ð0:75; 0:25Þ; ð0:125; 0:75Þ and ð0:875; 0:75Þ, for different numbers of elements. A convergence of the numerical solution towards the analytical solution for both the pressure and the flux when increasing the number of elements is observed, with the largest relative error when using 128  128 elements of 0:46%.

(59)

where j denotes a well, qx;j ; qy;j are the two components of the flux from a corner of the cell representing the well and ðDxÞj is the dimension of this cell. These conditions are specified using the normal pressure derivatives, which are the unknowns in the modified ‘q-based’ GEM for anisotropic media, as follows:

qp qp  K 12;j ¼ qx;j ; qx qy qp qp  K 22;j ¼ qy;j ; K 12;j qx qy

K 11;j

j ¼ 1; N w j ¼ 1; N w

(60)

4. Numerical results and discussion 4.1. Numerical results when known analytical solutions are available In order to investigate the performance of the numerical method proposed in Section 2, namely the ‘q-based’ GEM for anisotropic media, we consider some examples for which the analytical solution is known in the plane domain L ¼ fðx; yÞ : 0pxp1; 0pyp1g. In the first example, the permeability tensor is constant everywhere in the domain, while in the second example it has continuous heterogeneous components. Uniform grids are used in discretising the domain for these problems.

Fig. 2. The numerical solution obtained by the ‘q-based’ GEM for anisotropic media for 128  128 elements (—-) and the analytical solution (– – – –) for the pressure for Example 1.

ARTICLE IN PRESS P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

375

Table 1 Numerical results obtained for different numbers of elements and the analytical solution at the points ð0:25; 0:25Þ and ð0:75; 0:25Þ for Example 1 Mesh

88 16  16 32  32 64  64 128  128 Analytical

(0.25, 0.25)

(0.75, 0. 25)

p

qx

qy

p

qx

qy

0.566447 0.565604 0.565059 0.564677 0.564394 0.562500

0.505849 0.503925 0.503131 0.502664 0.502320 0.500000

3.770668 3.766934 3.764027 3.761938 3.760389 3.750000

1.066579 1.065727 1.065163 1.064765 1.064471 1.062500

0.494277 0.496283 0.497070 0.497516 0.497840 0.500000

3.272431 3.266921 3.263892 3.261816 3.260280 3.250000

Table 2 Numerical results obtained for different numbers of elements and the analytical solution at points ð0:125; 0:75Þ, and ð0:875; 0:75Þ for Example 1 Mesh

88 16  16 32  32 64  64 128  128 Analytical

(0.125, 0.75)

(0.875, 0.75)

p

qx

qy

p

qx

qy

1.519907 1.519013 1.518424 1.518007 1.517698 1.515625

0.746117 0.747633 0.748172 0.748454 0.748657 0.750000

3.856685 3.858378 3.861731 3.863753 3.865216 3.875000

2.269414 2.268600 2.268078 2.267712 2.267441 2.265625

0.754507 0.752821 0.752229 0.751899 0.751655 0.750000

3.101176 3.108620 3.111100 3.113139 3.114676 3.125000

4.1.2. Example 2: heterogeneous and anisotropic permeability In this example, the governing equation is given by (3), where the components of the heterogeneous permeability tensor are continuous functions defined now as

Fig. 3. The numerical solution obtained by the ‘q-based’ GEM for anisotropic media for 128  128 elements (—-) and the analytical solution (– – – –) for the pressure for Example 2.

Table 3 Numerical results obtained for different numbers of elements and the analytical solution at points (0.25, 0.25) and (0.625, 0.25) for Example 2 Mesh

K 11 ðx; yÞ ¼ 1 þ x þ y K 12 ðx; yÞ ¼ K 21 ðx; yÞ ¼ 0:5 þ y K 22 ðx; yÞ ¼ 1 þ 2x þ y

(65)

The external/internal forcing term for this example has the expression F ¼ ð3 þ x þ yÞex  cos y þ ð1 þ 2x þ yÞ sin y

(66)

88 16  16 32  32 64  64 128  128 Analytical

(0.25, 0.25)

(0.625, 0.25)

p

qx

qy

p

qx

qy

1.539175 1.537918 1.536856 1.536062 1.535462 1.531429

2.676129 2.661183 2.658772 2.658002 2.657403 2.652722

2.698647 2.681892 2.677823 2.675251 2.673190 2.658615

2.123522 2.122093 2.121048 2.120279 2.119697 2.115649

4.265559 4.244319 4.240541 4.239157 4.238067 4.229645

3.870906 3.861309 3.856220 3.851813 3.848335 3.823465

and the exact solution is given by p ¼ ex þ sin y

(67)

when the following boundary conditions are prescribed:

Table 4 Numerical results obtained for different numbers of elements and the analytical solution at points ð0:25; 0:75Þ and ð0:625; 0:75Þ for Example 2

pðx; 0Þ ¼ ex ;

Mesh

pðx; 1Þ ¼ ex þ sin 1;

qx ð0; yÞ ¼ ½1 þ y þ ð0:5 þ yÞ cos y;

0pxp1 0pyp1

qx ð1; yÞ ¼ ½ð2 þ yÞe þ ð0:5 þ yÞ cos y;

0pyp1

(68)

The components of the tensor K are averaged on each element using the formula Z 1 K ij;e ¼ K ðx; yÞ dx dy (69) Se LðeÞ ij where Se is the area of the element ðeÞ. By this averaging, the components of the tensor K are constant on each element, and they differ slightly from one element to another. The numerical results have been investigated for different numbers of elements. A comparison between the numerical and the analytical solutions for pressure for 128  128 elements is presented in Fig. 3. A very good agreement between the two solutions is clearly in evidence. The numerical results for the pressure and the two components of the flux for the points (0.25, 0.25), (0.625, 0.25), (0.25, 0.75) and (0.625, 0.75) in comparison to the analytical solution are illustrated in Tables 3 and 4, and a convergence of the numerical solution

88 16  16 32  32 64  64 128  128 Analytical

(0.25, 0.75)

(0.625, 0.75)

p

qx

qy

p

qx

qy

1.971001 1.970113 1.969391 1.968848 1.968435 1.965664

3.461814 3.469315 3.471858 3.473466 3.474692 3.482661

3.218356 3.218766 3.224016 3.228148 3.231199 3.251331

2.556442 2.555326 2.554462 2.553813 2.553319 2.549884

5.326674 5.336090 5.339527 5.341421 5.342808 5.351695

4.472126 4.481353 4.490034 4.495988 4.500390 4.530374

towards the analytical solution when increasing the number of elements is observed, with a relative error smaller than 0:66% when a maximum number of elements are used.

4.2. Numerical results for physical problems In this section physical problems in anisotropic media are investigated. Examples involving geological features, such as layers of material with different permeabilities, faults or wells, are

ARTICLE IN PRESS 376

P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

simulated with different boundary conditions imposed. In all these situations, no analytical solutions are available and the problems have been solved using the modified ‘q-based’ GEM for anisotropic media. To investigate the performance of this approach, we first tested it on problems with an available analytical solution. 4.2.1. Domains without wells The examples where there are no wells are governed by Eq. (3) in the absence of any external/internal forcing term. Pressure conditions are imposed on the top and bottom boundaries and noflow on the remainder boundaries. 4.2.1.1. Example 3: layers of different anisotropic permeabilities. This example simulates fluid flow in an anisotropic medium with three diagonal layers of different anisotropic permeabilities. The domain is represented by the rectangle L ¼ fðx; yÞ 2 R2 : 0pxp1; 0pyp2g, with the boundary conditions pðx; 0Þ ¼ 1; qx ð0; yÞ ¼ 0;

pðx; 2Þ ¼ 0;

0pxp1

(70) pffiffiffi The are delimited by the straight lines 3x  y  0:25 ¼ 0 players ffiffiffi and 3x pyffiffiffiþ 0:75 ¼ 0. For simplicity, we will refer to the layer satisfying 3x  y  0:25X0 as layer 1, the p middle one as layer 2 ffiffiffi and the top one, namely that satisfying 3x  y þ 0:75p0, as layer 3 (see Fig. 4). For all the elements used in the discretisation of the domain, the coordinates of their centre have been used to determine the layer to which the element belongs. The permeability tensor is defined everywhere in the domain by means of its principal values and by the angle of rotation. We denote by K 1;i and K 2;i and by yi the principal values of the permeability tensor and the angle of rotation corresponding to layer i, respectively, i ¼ 1; 3. We consider the following three situations: (i) K 1;1 (ii) K 1;1 K 1;3 (iii) K 1;1

¼ 6, ¼ 6, ¼ 2, ¼ 7,

qx ð1; yÞ ¼ 0;

0pyp2

K 1;2 ¼ 4, K 1;3 ¼ 2, K 2;i ¼ 1, yi ¼ p=3, i ¼ 1; 3, K 2;1 ¼ 1, y1 ¼ 5p=6, K 1;2 ¼ 4, K 2;2 ¼ 1, y2 ¼ p=3, K 2;3 ¼ 1, y ¼ 5p=6, K 1;2 ¼ 5, K 1;3 ¼ 2, K 2;i ¼ 1 and yi ¼ 5p=6; i ¼ 1; 3.

In all three cases, as expected, the fluid tends to locally flow in the direction of the largest component of the principal value of the permeability tensor. Specifically, in case (i) the fluid tends to flow towards the right in all layers. In case (ii) the flow tends to go towards the left in layer 1, then in the second layer it suddenly changes the direction, flowing virtually parallel to the vertical side walls of the medium in their vicinity, before changing its direction again in layer 3 back towards the left. In case (iii) the fluid maintains its initial direction towards the left. These tendencies for the fluid flow in each case are clearly evident in Fig. 4, representing the flux vector field. The pressure contours are illustrated in Fig. 5 and they show the changes in the pressure gradient at the interfaces between the layers. The results have been investigated using different numbers of elements and they are presented here for 64  128 elements, since for finer meshes they are expected to provide similar results. The flow is primarily influenced by the magnitude of the largest principal value and by the orientation of this. The arrows are almost proportional to this value and it is clearly in evidence in Fig. 4 how, for the same situation, the flux arrows change their magnitude according to the permeability of the layers. 4.2.1.2. Example 4: anisotropic medium containing a fault with isotropic permeability. The purpose of this example is to show the effect of the anisotropy of the medium on the fluid flow across the domain. The same domain and boundary conditions as in the previous example are imposed. One fault is considered in the

domain, by including cells of lower permeability, of thickness 2t f ¼ 103 and of length 0.75, located over ½0; 0:75  ½1  t f ; 1 þ t f . The medium, except for the fault cells, is considered to have the same permeability everywhere, which is given by ! 5 K 12 K¼ (71) K 12 5 where K 12 takes different values, varying between 4 and 4, namely 4; 3; 2; 1; 0; 1; 2; 3; 4. In each of these cases, the permeability of the fault is the same, namely ! 0 102 Kf ¼ (72) 0 102 The anisotropy of the medium has an important influence on the fluid flow. Fig. 6 illustrates the flux vectors for different situations, namely when K 12 ¼ 4; K 12 ¼ 2; K 12 ¼ 2 and K 12 ¼ 4. There is an obvious tendency of the fluid flow to assume a different direction depending on the permeability of the medium. It can be seen that as the component K 12 approaches extreme values (4 and 4), the flux field deviates from the direct vertical flow to a larger extent. This is most apparent from the arrows close to the top and the bottom boundaries. The fluid avoids the fault with a different intensity in each case. The pressure plots, given in Fig. 7, reflect the changes in the directions of the fluid flow. The results presented here have been obtained for 64  129 elements. To have an even better image of how the anisotropy of the medium influences the fluid flow pattern, two problems based on Example 4 (with the fault position and the boundary conditions maintained the same) have been investigated. The first problem maintains a constant thickness for the fault, namely 2t f ¼ 103, and varies its permeability, so that the diagonal components K 11;f ¼ K 22;f have the following values: 101 ; 102 ; 103 and 104 , with K 12;f ¼ 0 in each case. In each of these situations, the permeability of the medium is given all the values above, plus the two extra values of K 12 ¼ 4:9 and 4.9, chosen so that their value is close to the limiting cases of 5 and 5, respectively. In this way, the effect of the anisotropy on the fluid flow for different permeabilities of the fault can be distinguished. The second problem consists in maintaining the permeability of the fault as constant, namely K 11;f ¼ K 22;f ¼ 102, K 12;f ¼ 0, and changing the thickness of the fault to the values 101 ; 102 ; 103 . As before, K 12 is given all the above-mentioned values from 4:9 to 4:9. For both problems, the bulk permeability is calculated using the formula given by qav ¼ K bulk

Dp ly

(73)

where qav is the flux (fluid velocity) averaged over the total of cross section (the top boundary in this case) and Dp is the pressure drop over the domain length ly ¼ 2 in the present situation. The bulk permeability is plotted as a function of K 12 =K 11 , so that a value of zero on the horizontal axis refers to the isotropic situation. The results have been obtained for 16  33 elements. Fig. 8 illustrates the situation when the permeability of the fault is changed and the thickness remains constant. It can be seen that the bulk permeability tends to increase as the fault permeability increases for any medium permeability. The differences in the bulk permeability values are the largest in the vicinity of the isotropic medium permeability and they become smaller towards the extremities where K 12 reaches the values 4:9 and 4:9. Fig. 9 represents the plot of the bulk permeability when the permeability of the fault is constant and the thickness varies. There is a tendency of the bulk permeability to increase as the fault becomes thinner.

ARTICLE IN PRESS P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

377

Fig. 4. Representation of the flux vectors obtained by the modified ‘q-based’ GEM for anisotropic media for 64  128 elements for (a) case (i), for (b) case (ii) and for (c) case (iii) for Example 3.

ARTICLE IN PRESS 378

P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

Fig. 5. The numerical results obtained by the modified ‘q-based’ GEM for anisotropic media for 64  128 elements (—-) for the pressure for (a) case (i), (b) case (ii) and (c) case (iii) for Example 3.

ARTICLE IN PRESS P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

379

Fig. 6. Representation of the flux vectors obtained by the modified ‘q-based’ GEM for anisotropic media for 64  129 elements when the permeability of the media is given by K 11 ¼ K 22 ¼ 5 and (a) K 12 ¼ 4, (b) K 12 ¼ 2, (c) K 12 ¼ 2 and (d) K 12 ¼ 4 for Example 4.

ARTICLE IN PRESS 380

P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

Fig. 7. The numerical results obtained by the anisotropic modified ‘q-based’ GEM for 64  129 elements (—-) for the pressure when the permeability of the media is given by K 11 ¼ K 22 ¼ 5 and (a) K 12 ¼ 4, (b) K 12 ¼ 2, (c) K 12 ¼ 4 and (d) K 12 ¼ 2 for Example 4.

4.2.1.3. Example 5: layers of different anisotropic permeabilities across a fault. This example is a simulation of fluid flow in an anisotropic porous medium when there is a fault (parallel to the Ox-axis) across the whole of the middle of the domain. The domain is taken to be the rectangle L ¼ fðx; yÞ 2 R2 : 0pxp1:3;

0pyp2g, with the boundary conditions pðx; 0Þ ¼ 1; qx ð0; yÞ ¼ 0;

pðx; 2Þ ¼ 0;

0pxp1:3

qx ð1:3; yÞ ¼ 0;

0pyp2

(74)

ARTICLE IN PRESS P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

381

Fig. 8. Representation of the bulk permeability K bulk as a function of K 12 =K 11 for different values of the fault permeability: K 12;f ¼ 0 and K 11;f ¼ 101 (-  - -  - ), K 11;f ¼ 102 (—-), K 11;f ¼ 103 ð    Þ and K 11;f ¼ 104 (– – – –), using 16  33 elements for 2t f ¼ 103 and K 12 varying between 4:9 and 4:9.

Fig. 9. Representation of the bulk permeability K bulk as a function of K 12 =K 11 for different fault thicknesses: 103 (—-), 102 (-  -  -  - ), and 101 ð    Þ using 16  33 elements for constant permeability of the fault given by K 11;f ¼ 102 ; K 12;f ¼ 0, and K 12 varying between 4:9 and 4:9.

The fault has a thickness of 2t f ¼ 103 and is located throughout ½0; 1:3  ½1  t f ; 1 þ t f . The medium consists of layers of three different permeabilities. Below the fault, p three layers are conffiffiffi sidered, delimited by the straight lines 3 x  y p0:9 pffiffiffi ffiffiffi ¼ 0 and 3x  y ¼ 0. These layers are referred to as follows: 3x  yp0 is

layer 1, the middle one is layer 2 and the remaining, pffiffiffi 3x  y  0:9X0, is layer 3 (see Fig. 10). Above p the ffiffiffi fault, there are also pffiffiffithree layers, delimited by the lines p 3ffiffiffi x  y  0:25 ¼ 0 and 3x  y þ 0:65 ¼ 0. The layer satisfying 3x  y þ 0:65p0 is referred to as layer 4, the

ARTICLE IN PRESS 382

P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

Fig. 10. Representation of (a) the pressure and (b) flux vectors obtained by the modified ‘q-based’ GEM for anisotropic media with a fault of permeability K 11;f ¼ 103 ; K 22;f ¼ 1, K 12;f ¼ 0 and 80  129 elements for Example 5.

pffiffiffi middle one is layer 5 and the remaining, 3x  y  0:25X0, is layer 6. The permeability of the layers is given by the principal values and the rotation angle, which is the same for each layer, namely y ¼ p=3. The permeabilities of layers 1 and 4 are the same as those of layers 2 and 5 and of layers 3 and 6 and they are all given by the principal values, as follows: K 1;1 ¼ K 1;4 ¼ 3;

K 2;1 ¼ K 2;4 ¼ 1

K 1;2 ¼ K 1;5 ¼ 5;

K 2;2 ¼ K 2;5 ¼ 1

K 1;3 ¼ K 1;6 ¼ 2;

K 2;3 ¼ K 2;6 ¼ 1

(75)

The layers above and below the fault are thus considered as being displaced by the fault, see Fig. 10. Two different permeabilities of the fault are considered. First, the permeability is given by K 11;f ¼ 103, K 22;f ¼ 1, K 12;f ¼ 0. In this case, the fluid is allowed to pass through the fault easily, as the fault is much more permeable in the y direction than in the x direction. This is evident in Fig. 10 where the pressure and the flux vector fields are represented for 80  129 elements. The permeability of the layers influences the magnitude of the flux vectors. The flow in the middle layers, from both above and below the fault, is more intense than in the other layers, as these have a higher permeability, reflected by their principal values. The second investigated case for this example is when the permeability of the fault has the components K 11;f ¼ 1, K 22;f ¼ 103 , K 12;f ¼ 0, so the fault is much less permeable in the y direction than in the x direction. All the other conditions remain the same. In this situation, much less fluid is allowed to pass across the fault. The magnitude of the overall flux vectors is smaller than that in the previous case. Fig. 11(b) illustrates the difficulty the fluid has in passing through the fault imposed by the low permeability in the y direction. The fluid also has a tendency to flow mainly in the middle layers where the permeability is higher. The pressure for this example is represented in Fig. 11(a).

4.2.1.4. Example 6: layers of anisotropic permeability with an overlap area across a fault. Example 6 again simulates a fluid flow in an anisotropic medium containing a fault, with the medium consisting of layers of different permeabilities. The same domain and boundary conditions are considered as in the previous example. The fault is located across the domain, over ½0; 1:3  ½1  t f ; 1 þ t f , where the thickness is 2t f ¼ 102 . However, in this example the fault has a different permeability along its length, namely in the middle of the fault, over ½0:26; 1:04  ½1  t f ; 1 þ t f , the permeability components are K 11;f ¼ 1; K 22;f ¼ 0:05; K 12;f ¼ 0, while in the rest of the fault K 22;f ¼ 103 , with the other components the same as in the middle section of the fault. There are two layers of the same permeability in the domain. pffiffiffi One is below the fault, delimited by the straight lines 3 xy pffiffiffi 0:6 ¼ 0 and 3x  y  0:1 ¼ 0 (see Fig. 12). p The ffiffiffi other one is above the pffiffiffi fault, delimited by the straight lines 3x  y  0:2 ¼ 0 and 3x  y þ 0:3 ¼ 0. Hence, if the fault was ignored, the two central layers would have a small overlap area. For these two layers the permeability is the same, given by the principal values K 1 ¼ 5, K 2 ¼ 3 and the angle of rotation y ¼ p=3. In the rest of the domain, the permeability is defined by K 1 ¼ K 2 ¼ 104 and the same angle of rotation. In these conditions, the fluid is all forced to flow in the middle layers, as is clearly in evidence in Fig. 12(b), where the flux vectors are represented when using 80  129 elements. At the centre, the fluid is forced to go through the part of the fault which shares a common part with both the middle layers. In the rest of the domain the flow has an extremely low intensity due to the low permeability imposed. The corresponding pressure contours are plotted in Fig. 12(a).

4.2.1.5. Example 7: layers of anisotropic permeability with no overlap area. In this example all the conditions remain as for Example 6, except for the location of the central layer above the fault. This is

ARTICLE IN PRESS P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

383

Fig. 11. Representation of (a) the pressure and (b) flux vectors obtained by the modified ‘q-based’ GEM for anisotropic media with a fault of permeability K 11;f ¼ 1; K 22;f ¼ 103 , K 12;f ¼ 0 and 80  129 elements for Example 5.

Fig. 12. Representation of (a) the pressure and (b) flux vectors obtained by the modified ‘q-based’ GEM for anisotropic media with 80  129 elements for Example 6.

now moved towards the left p and pffiffiffi ffiffiffi is now delimited by the straight lines 3x  y þ 0:4 ¼ 0 and 3x  y þ 0:9 ¼ 0 (see Fig. 13). Thus there is no overlap between the central layers across the fault, a factor that has a huge influence on the fluid flow. The fluid is still forced to flow through the two layers, but with a much lower intensity. When the fluid reaches the fault virtually all the fluid is forced to travel along it until it reaches the more permeable layer above the fault, due to the low permeability of the rest of the

medium. The pressure and the flux vector field are illustrated in Fig. 13, when using 80  129 elements.

4.2.2. Domains with wells The problems in this section investigate the fluid flow in an anisotropic medium in a domain with sealing boundary conditions. The domain is represented in each case by the rectangle

ARTICLE IN PRESS 384

P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

Fig. 13. Representation of (a) the pressure and (b) flux vectors obtained by the modified ‘q-based’ GEM for anisotropic media with 80  129 elements for Example 7.

Fig. 14. Representation of (a) the pressure and (b) flux vectors obtained by the modified ‘q-based’ GEM for anisotropic media using 82  162 elements for Example 8.

L ¼ fðx; yÞ 2 R2 : 0pxp1; 0pyp2g with the boundary conditions qy ðx; 0Þ ¼ 0;

qy ðx; 2Þ ¼ 0;

0p xp1

qx ð0; yÞ ¼ 0;

qx ð1; yÞ ¼ 0;

0pyp2

(76)

Two wells are considered in the domain in each case. One well is located at ð0:25; 1:75Þ and is considered as a source (injection) with strength q1 ¼ 1, and the other well is considered as a sink (production) with strength q2 ¼ 1 and is located at ð0:75; 0:25Þ.

ARTICLE IN PRESS P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

385

Fig. 15. Representation of the flux vectors obtained by the modified ‘q-based’ GEM for anisotropic media with 82  148 elements for the permeability of the media K 11 ¼ K 22 ¼ 5 and (a) K 12 ¼ 4, (b) K 12 ¼ 2, (c) K 12 ¼ 2 and (d) K 12 ¼ 4 for Example 9.

ARTICLE IN PRESS 386

P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

4.2.2.1. Example 8: wells within diagonal layers of anisotropic permeability. In this example the medium is represented by three layers of different layers are delimited by the pffiffiffi permeabilities. Thep ffiffiffi straight lines 3x  y  0:25 ¼p0ffiffiffiand 3x  y þ 0:75 ¼ 0. We will refer to the layer satisfying 3x  y  0:25X0 as layer 1, the middle one as layer 2 and the top one, namely that satisfying pffiffiffi 3x  y þ 0:75p0, as layer 3 (see Fig. 14). We will denote by K 1;i and K 2;i and by yi the principal values of the permeability tensors

and the angle of rotation corresponding to layer i, respectively, i ¼ 1; 3, and we consider the following values for these: K 1;1 ¼ 6;

K 1;2 ¼ 4;

K 1;3 ¼ 2;

K 2;i ¼ 1;

yi ¼

5p , 6

i ¼ 1; 3. It can be seen that, as the angle of rotation is the same in the three layers, the flow maintains similar orientations in all the layers,

Fig. 16. Representation of the pressure obtained by the modified ‘q-based’ GEM for anisotropic media with 82  148 elements for the permeability of the media K 11 ¼ K 22 ¼ 5 and (a) K 12 ¼ 4, (b) K 12 ¼ 2, (c) K 12 ¼ 2 and (d) K 12 ¼ 4 for Example 9.

ARTICLE IN PRESS P. Lorinczi et al. / Engineering Analysis with Boundary Elements 33 (2009) 368–387

with only slight variations in flow direction from the source to the sink. The pressure and the flux vectors are represented in Fig. 14, when using 82  162 elements. 4.2.2.2. Example 9: wells in anisotropic media containing faults. The last example is a simulation of wells in a domain containing two faults. It presents a comparison between several situations when the permeability of the medium changes and all the other conditions are the same. Each fault is represented by a thin layer of thickness 2t f ¼ 102 and of length 0.8 by including cells of lower permeability in the domain. One fault is located over ½0; 0:8  ½1:5  t f ; 1:5 þ t f  and the other one over ½0:2; 1  ½0:5  t f ; 0:5 þ t f . The faults are of the same permeability, given by K 11;f ¼ K 22;f ¼ 103 ; K 12;f ¼ 0. Over the rest of the domain a constant permeability is imposed. Four situations are investigated, namely ! 5 K 12 K¼ (77) K 12 5 where K 12 takes the values 4, 2, 2 and 4. The flux vector field for each of these situations is presented in Fig. 15. The difference in flow for the permeabilities considered can be noticed around the wells, near the edges of the faults and between the faults. It can be seen that when K 12 ¼ 4 the arrows have a similar magnitude in the region between the faults, while for K 12 ¼ 4 the arrows from the centre of the domain have a larger magnitude than those in the vicinity of the boundaries. The pressure contours change their orientation with the permeability of the medium, as shown in Fig. 16.

5. Conclusions In this article we have presented a numerical technique for solving steady-state flow problems in anisotropic media. The ‘q-based’ GEM for anisotropic media is obtained by applying the theory from the BEM for such media over FEM meshes. It uses the same concept as the ‘q-based’ GEM for isotropic media, with the number of unknowns at each internal node remaining at three. The ‘q-based’ GEM for anisotropic media can be used for problems in which the permeability is homogeneous or has continuous components. In the latter case every component of the permeability tensor is averaged on each element. However, if the permeability tensor is greatly different within the domain, namely when a discontinuity in the permeability of the media occurs, then the ‘q-based’ GEM for anisotropic media is no longer adequate. A modified ‘q-based’ GEM has been developed for such situations, which is based on satisfying a nodal flux condition at each node and on the continuity of the tangential pressure gradient across the element boundaries. The ‘q-based’ GEM for anisotropic media has been successfully used for solving problems for which an analytical solution was

387

available, namely problems involving constant and continuous permeability tensor components. The numerical results were found to be in very good agreement with the exact solutions. Physical problems with discontinuous permeabilities have been investigated using the modified ‘q-based’ GEM for anisotropic media. Geological features such as layers of different permeabilities, faults and wells have been modelled. The impact of the different permeabilities of the media has been illustrated in each case for both the pressure and the flux vector fields. In summary, from the numerical results obtained it can be concluded that the ‘q-based’ GEM for anisotropic media with the extension presented here is a suitable and reliable technique for solving physical steady-state problems in such media. Discontinuities in the media permeabilities can be naturally treated, which makes the present numerical method a useful tool for geological anisotropic problems.

Acknowledgement Piroska Lorinczi would like to acknowledge the financial support received from the ORS and the Rock Deformation Research Group, University of Leeds, UK. References [1] Pecher R, Harris SD, Knipe RJ, Elliott L, Ingham DB. New formulation of the Green element method to maintain its second-order accuracy in 2D/3D. Eng Anal Boundary Elem 2001;25:211–9. [2] Taigbenu AE. The Green element method. Boston: Kluwer Academic Publishers; 1999. [3] Taigbenu AE. Improvements in heat conduction calculations with flux-based green element method. In: Advances in boundary integral methods— proceedings of the fifth UK conference on boundary integral methods, The University of Liverpool, Liverpool, 2005, p. 190–9. [4] Popov V, Power H. The DRM-MD integral equation method for the numerical solution of convection–diffusion equation. In: Proceedings of the European boundary element method symposium, 1998, p. 67–80. [5] Popov V, Power H. The DRM-MD integral equation method: an efficient approach for the numerical solution of domain dominant problems. Int J Numer Methods Eng 1999;44(3):327–53. [6] Lorinczi P, Harris SD, Elliott L. Modelling of highly-heterogeneous media using a flux-vector-based green element method. Eng Anal Boundary Elem 2006;30:818–33. [7] Chang YP, Kang CS, Chen DJ. The use of fundamental Green’s functions for the solution of heat conduction in anisotropic media. Int J Heat Mass Transfer 1973;16:1905–18. [8] Garroni MG, Menaldi JL. Green functions for second order parabolic integrodifferential problems. Harlow: Longman Scientific and Technical; 1992. [9] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques: theory and applications in engineering. Berlin: Springer; 1984. [10] Symm GT, Pitfield RA. Solutions of the Laplace’s equation in two dimensions. NPL Report NAC 44; 1974. [11] Mera NS. A comparison of boundary element method formulations for steady state anisotropic heat conduction problems. Eng Anal Boundary Elem 2001; 25:115–28. [12] Lorinczi P. Applications of the Green element method to flow in heterogeneous porous media. PhD thesis, University of Leeds, 2006. [13] Shewchuk JR. An introduction to the conjugate gradient method without the agonizing pain. Technical Report CMU-CS-94-125, School of Computer Science, Carnegie Mellon University, Pittsburgh, PA; 1994.