Fourier Series Expansion in a Non-Orthogonal System of Coordinates for the Simulation of 3D DC Borehole Resistivity Measurements
D. Pardo[a], V. M. Calo[b], C. Torres-Verd´ın[a], and M. J. Nam[a] a Department b Institute
of Petroleum and Geosystems Engineering, The University of Texas at Austin
for Computational Engineering and Sciences (ICES), The University of Texas at Austin
Abstract We describe a new method to simulate 3D borehole resistivity measurements at zero frequency (DC). The method combines the use of a Fourier series expansion in a non-orthogonal system of coordinates with an existing 2D goal-oriented higher-order self-adaptive hp-Finite Element algorithm. The new method is suitable for simulating measurements acquired with borehole logging instruments in deviated wells. It delivers high-accuracy simulations and it enables a considerable reduction of the computational complexity with respect to available 3D simulators, since the number of Fourier modes (basis functions) needed to solve practical applications is limited (typically, below 10). Furthermore, numerical results indicate that the optimal 2D grid based on the 0-th Fourier mode (also called central Fourier mode) can be employed to efficiently solve the final 3D problem, thereby, avoiding the expensive construction of optimal 3D grids. Specifically, for a challenging through-casing resistivity application, we reduce the computational time from several days (using a 3D simulator) to just two hours (with the new method), while gaining accuracy. The new simulation method can be easily extended to different physical phenomena with similar geometries, as those arising in the simulation of 3D borehole electrodynamics and sonic (acoustics coupled with elasticity) measurements. In addition, the method is especially suited for inversion, since we demonstrate that the number of Fourier modes needed for the exact representation of the materials is limited to only one (the central mode) for the case of borehole measurements acquired in deviated wells. Key words: Fourier Series Expansion, Non-Orthogonal System of Coordinates, hp-FEM, Goal-Oriented Adaptivity, Borehole Measurements
Preprint submitted to Elsevier
5 December 2007
1
2 3 4 5 6 7 8 9 10 11 12 13
14 15 16 17 18 19 20 21 22 23 24 25
26 27 28 29 30
1
Introduction
Since the Schlumberger brothers acquired the first borehole resistivity measurement in 1927, borehole logging measurements are widely used by the oilindustry for hydrocarbon reservoir characterization and surveillance. A variety of new logging instruments has been developed and employed during the last decades in virtually all existing hydrocarbon reservoirs worldwide. Despite the success of well-logging measurements, the planning and drilling of a single well may cost several millions of dollars, and the corresponding results (logs) may sometimes be difficult to interpret. To improve the interpretation of results, and thus, to better quantify and determine the existing subsurface materials and increase hydrocarbon recovery, diverse numerical methods have been developed to perform borehole simulations as well as to invert well-log measurements. Most numerical methods used by the oil-industry are based on 1D and 2D algorithms. 1D results provide fast interpretation of subsurface material properties, enabling real-time modifications on the well trajectory in the case of logging-while-drilling (LWD) instruments. Despite the fact that 1D results are typically “corrected” (modified) using semi-analytical formulas to account for modeling of 2D and 3D geometries, their accuracy is compromised in the presence of logging instruments, mud-filtrate invasion effects, anisotropy, casing, etc. [1]. On the contrary, 2D axial-symmetric simulations (see [2–4]) enable accurate modeling of logging measurements, invasion effects, anisotropy, and casing, provided that both the geometry and sources are axial-symmetric. This is an important restriction that implies that the well trajectory is vertical, that is, perpendicular to the subsurface material layers. If the source is not axial-symmetric, it is possible to employ a Fourier series expansion for the source, and to solve the resulting sequence of problems (one problem for each Fourier mode) independently using a 2D axial-symmetric simulator. This method involves independent solutions of various 2D problems, and thus, it is referred as 2.5D method (see [5]).
38
For general 3D geometries, such as those resulting from simulations of deviated wells, a number of simulators have been developed (see, for example, [6–13]). Despite the sophistication of some of those methods, they all have failed to provide accurate results in a limited amount of time, due to the high complexity associated with 3D simulations in arbitrary geometries. Nevertheless, it is becoming increasingly important to accurately and efficiently simulate various logs in deviated wells, since highly deviated wells span longer distances within hydrocarbon layers, thereby providing a higher level of hydrocarbon recovery.
39
The main technical contribution of this paper is that we take advantage of the
31 32 33 34 35 36 37
2
47
fact that typical geometries arising in the simulation of measurements acquired in deviated wells (see Fig. B.1) are almost two-dimensional if we consider a particular non-orthogonal system of coordinates (defined in Section 2.2). A first approach toward using non-orthogonal systems of coordinates was developed by [14], where they describe an oblique coordinate system suitable for deviated wells. However, the applicability of their work is limited because they assume that the borehole is devoid of materials, that is, no logging instrument is present.
48
[Fig. 1 about here.]
40 41 42 43 44 45 46
62
In this paper, we describe a new general method for the numerical simulation of resistivity logging measurements acquired in deviated wells. We first divide the domain into three-different subdomains: (1) the logging instrument, where we employ a cylindrical system of coordinates, (2) the formation material, where we employ an oblique system of coordinates, and (3) the borehole segment located between the logging instrument and the formation material (containing fluid and possibly steel casing), where we construct a system of coordinates intended to match the two previously defined systems of coordinates and such that the resulting change of coordinates is globally continuous, bijective, and invertible (see Fig. B.2 for additional geometrical details). We note that a deviated well penetrating horizontal layers (Fig.B.1) is simply a rotation of a vertical well penetrating dipping layers (Fig. B.2). Thus, the important feature in these type of problems is the dip angle between the well and the layers in the formation.
63
[Fig. 2 about here.]
49 50 51 52 53 54 55 56 57 58 59 60 61
64 65 66 67 68 69 70 71 72 73 74 75 76 77
78 79
After defining the above non-orthogonal system of coordinates, we notice that material properties are constant along one direction (the “quasi-azimuthal” direction). In addition, the metric [15] associated with the change of coordinates from a reference 2D grid to the physical 3D geometry can be decomposed exactly (without approximations) in terms of only five Fourier modes in the quasi-azimuthal direction. Thus, a total of five Fourier modes are necessary to account for all material and geometrical properties. This surprising fact implies that the resulting stiffness matrix is penta-diagonal in terms of the Fourier modal coefficients. Furthermore, numerical results indicate that, in most applications, only a few modes (fewer than ten) are necessary to obtain an adequate approximation to the exact solution. This new method provides a dramatic reduction on the computational complexity with respect to conventional 3D simulators. The method is suitable for forward and inverse problems, as well as for multi-physic applications. In the remainder of this paper, we analyze the following topics: In Section 2, we formally derive the new method outlined above, and present the final varia3
80 81 82 83 84 85
tional formulation. The resulting formulation requires several modifications on an existing 2D self-adaptive goal-oriented hp-Finite Element Method (FEM), which are described in detail in Section 3. Numerical results included in Section 4 are intended to validate the method and to study its applicability to everyday logging-operations. Further applications of this method are discussed in Section 5. The main conclusions of this paper are summarized in Section 6.
89
This paper also incorporates two appendices: The first one describes the Fourier series modal coefficients of the metric associated with the change of coordinates for deviated wells. The second one describes a change of coordinates suitable for borehole-eccentered measurements.
90
2
86 87 88
91 92 93 94 95 96 97 98 99
100
101 102
In this Section, we first derive a variational formulation for electromagnetics at zero frequency. Second, we introduce a non-orthogonal system of coordinates suitable for deviated wells in a borehole environment, and we discuss its main properties. Third, we describe a variational formulation for zerofrequency Maxwell’s equations in the new system of coordinates. Then, we employ a Fourier series expansion in the quasi-azimuthal direction to derive the corresponding variational formulation in terms of the Fourier modal coefficients. Finally, we briefly describe the method employed for solving the resulting formulation via a self-adaptive goal-oriented hp-FEM.
2.1
105 106 107 108 109 110
Variational Formulation
At DC (i.e., zero frequency), the electromagnetic phenomena (governed by Maxwell’s equations) reduces to the so called conductive media equation, i.e., ∇ · (σ∇u) = −∇ · Jimp ,
103
104
Method
(1)
where σ 6= 0 ∈ L∞ (Ω) 1 is the conductivity tensor, Jimp represents the prescribed impressed electric current sources, and u is the scalar electric potential. The above equation should be understood in the distributional sense, i.e. it is satisfied in the classical sense in subdomains of regular material data, but it also implies appropriate interface conditions across material interfaces. We note that the electric field is given by E = −∇u in the case of simply connected domains. 1
If σ = 0 in part of the domain, see [16].
4
111 112 113
To derive the variational formulation for the conductive media equation, we first define the L2 -inner product of two (possibly complex- and vector-valued) functions g1 and g2 as hg1 , g2 iL2 (Ω) =
114
Z
g1∗ g2 dV ,
(2)
Ω
115
116 117 118 119
where g1∗ denotes the adjoint (conjugate transpose) of function g1 . By multiplying test function v ∈ H 1 (Ω) = {u ∈ L2 (Ω) : ∇u ∈ L2 (Ω)} by equation (1), and by integrating by parts over the domain Ω ⊂ R 3 , we obtain the following variational formulation after incorporating the essential (Dirichlet) boundary condition (BC): Find
1 u ∈ uD + HD (Ω) such that:
h∇v
, σ∇uiL2 (Ω) = v , ∇ · Jimp
120
D
(3)
E L2 (Ω)
+ hv , hiL2 (ΓN )
1 ∀v ∈ HD (Ω) ,
124
where uD is a lift (typically uD = 0) of the essential BC data uD (denoted with the same symbol), h = σ∇u·n is a prescribed flux on ΓN , n is the unit normal 1 outward (with respect to Ω) vector, and HD (Ω) = {u ∈ H 1 (Ω) : u|ΓD = 0} is the space of admissible test functions associated with problem (3).
125
2.2
121 122 123
126 127 128
Non-Orthogonal Coordinate System for Deviated Wells
For deviated wells, as the one described in Fig. B.2, we introduce the following non-orthogonal coordinate system ζ = (ζ1 , ζ2 , ζ3 ) in terms of the Cartesian coordinate system x = (x1 , x2 , x3 ): x1
129
x2
= ζ1 cos ζ2 = ζ1 sin ζ2
,
(4)
x3 = ζ3 + θ0 f1 (ζ1 ) cos ζ2
5
130 131
where θ0 = tan θ, θ is the dip angle 2 , and f1 is defined for given values ρ1 and ρ2 as: 0 ζ
1 − ρ1 ρ2 f1 (ζ1 ) = f1 = ρ2 − ρ1 ζ 1
132
133
ζ1 < ρ1
137 138 139 140 141
142 143
ζ1 > ρ2
ζ1 < ρ1
ρ2 f10 (ζ1 ) = f10 = ρ2 − ρ1 1
134
136
(5)
with the corresponding derivative given by 0
135
ρ1 ≤ ζ1 ≤ ρ2 ,
ρ1 < ζ1 < ρ2 .
(6)
ζ1 > ρ2
Intuitively, ρ1 is defined so that ζ1 < ρ1 covers “subdomain I” of Fig. B.2. In this subdomain, we have defined a cylindrical system of coordinates. Additionally, ρ2 is defined so that ζ1 > ρ2 covers “subdomain III” of Fig. B.2. We employ an oblique system of coordinates over this subdomain. Finally, “subdomain II” of Fig. B.2 is intended to “glue” subdomain I with subdomain III so the resulting system of coordinates is globally continuous, bijective, and exhibits a positive Jacobian. 1 ,x2 ,x3 ) and its inverse associated with the above The Jacobian matrix J = ∂(x ∂(ζ1 ,ζ2 ,ζ3 ) change of coordinates are given by:
(
J =
144
∂xi ∂ζj
)
= i,j=1,2,3
−ζ1 sin ζ2
cos ζ2 sin ζ2
ζ1 cos ζ2
θ0 f10 cos ζ2
−θ0 f1 sin ζ2
0 0
(7)
, and
1
J
145
−1
=
cos ζ2 sin ζ2 − ζ1 f1 sin2 ζ2
−θ0
146
ζ1
!
+
f10
2
cos ζ2
sin ζ2 cos ζ2 ζ1 θ0 sin ζ2 cos ζ2
0 f1 − f10 ζ1
!
0 ,(8)
1
respectively, where det(J ) = |J | = ζ1 . 2
The dip angle is defined in borehole geophysics as the angle between the well trajectory and a normal vector to the layer boundaries. For example, if formation layers are horizontal, a ninety-degree dip angle corresponds to a horizontal well.
6
147 148
The corresponding metric tensor G = {G nm }n,m=1,2,3 = J T J is given by (see [15] for details about metrics)
G=
149
(θ0 cos ζ2 f10 )2
−θ02 f1 f10
1 + −θ02 f1 f10 sin ζ2 cos ζ2
ζ12 + (θ0 f1 sin ζ2 )2
θ0 f10 cos ζ2
150
0
1/ζ12
154 155
156 157 158
cos ζ2
, −θ0 f1 sin ζ2 1
(9)
0 G −1 =
151
153
−θ0 f1 sin ζ2
with its inverse G −1 = {G nm }n,m=1,2,3 given by: 1
152
sin ζ2 cos ζ2
θ0 f10
−θ0 f10 cos ζ2
θ0 f1
−θ0 f10 cos ζ2 sin ζ2 .(10) θ0 f1 2 ζ1 f 1 1 + θ02 (( sin(ζ2 ))2 + (cos(ζ2 )f10 )2 ) ζ1
sin ζ2 ζ12
Appendix A provides formulas for all the Fourier modal coefficients of tensor metric G and its inverse with respect to variable ζ2 . Here, we emphasize that only five Fourier modal coefficients are necessary to exactly reproduce the tensor matrix (and its inverse) associated with the above change of coordinates. In summary, the described change of coordinates has three essential properties that make it suitable and attractive for simulations of resistivity measurements along deviated wells:
163
• It is globally continuous, bijective, and with positive Jacobian. • Material properties are constant with respect to the quasi-azimuthal direction ζ2 . • Only five Fourier modes in terms of ζ2 are necessary to reproduce the tensor metric and its inverse (see Appendix A).
164
2.3
159 160 161 162
Variational Formulation in an Arbitrary System of Coordinates
171
Given an arbitrary (possibly non-orthogonal) system of coordinates ζ = (ζ1 , ζ2 , ζ3 ), as the one defined above, in this subsection we derive the corresponding variational formulation for the conductive media equation. By selecting the Cartesian coordinate system x = (x1 , x2 , x3 ) as our reference system of coordinates, our change of coordinates is described by the mapping x = ψ(ζ), which is assumed to be bijective, with positive Jacobian determinant, and globally continuous (see [17], Chapter XII).
172
Given arbitrary scalar-valued functions u = u(x), v = v(x), we denote u˜ :=
165 166 167 168 169 170
7
173 174
175
176 177
u ◦ ψ = u˜(ζ), and v˜ := v ◦ ψ = v˜(ζ). Thus, making use of the chain rule, we obtain ∇u =
3 X
∂ u˜ ∂ζn ˜ T ∂u exi = J −1 , ∂ζ i,n=1 ∂ζn ∂xi
∂ u˜ ∂ u˜ is the vector with the n-th component being , and exi is the ∂ζ ∂ζn unit vector in the xi -direction. Therefore, where
*
h∇v , σ∇uiL2 (Ω) = J 178
*
=
179
−1 T
+
˜ ∂˜ v T ∂u ˜ −1 , σJ ∂ζ ∂ζ
∂˜ v ˜ T ∂u ˜ −1 , J −1 σJ ∂ζ ∂ζ
L2 (Ω)
(12)
+
, L2 (Ω)
˜ = σ ◦ ψ. Similarly, we obtain where σ D
180
(11)
E
hv, f iL2 (Ω) = v˜, f˜
, and
L2 (Ω)
D
E
˜ hv, hiL2 (ΓN ) = v˜, h
L2 (ΓN )
,
(13)
181
˜ = h ◦ ψ. where f = ∇ · Jimp , f˜ = f ◦ ψ, and h
182
Extending the ideas of [18] to the electrostatic case, we introduce the tensor T
183
184
185
186 187 188
189 190
˜ N EW := J −1 σJ ˜ −1 |J |, σ
(14)
and functions f˜N EW := f˜|J |
;
˜ N EW := h|J ˜ S |, h
(15)
where |J S | is the determinant of the Jacobian associated with the change of variables corresponding to 2D surface ΓN . Our new space of admissible test v T ∂˜ 1 ˜D functions is given by H (Ω) = {˜ v ∈ L2 (Ω) : v˜|Γ˜ D = 0 , J −1 ∈ L2 (Ω)}. ∂ζ ˜ = Ω ◦ ψ and dropping the ˜ symbol from the notaThen, integrating over Ω tion, we arrive at the original variational formulation (3) in terms of our new 8
191
coordinate system, with new material and load data, namely, Find * ∂v
1 u ∈ uD + HD (Ω) such that:
∂u , σ N EW ∂ζ ∂ζ
192
+
= hv , fN EW iL2 (Ω) + hv , hN EW iL2 (ΓN ) 1 ∀v ∈ HD (Ω) ,
193 194 195 196 197 198
(16)
L2 (Ω)
where our L2 inner-product definition does not include the determinant of the Jacobian |J | corresponding to the change of variables, since information about the determinant of the Jacobian |J | is already included in the new material coefficient σ N EW and load data fN EW and hN EW . Thus, for arbitrary functions g1 and g2 defined on the ζ-coordinate system, our inner-product is described as hg1 , g2 iL2 (Ω) =
199
Z
g1∗ g2 dζ1 dζ2 dζ3 .
(17)
Ω
200
201 202 203 204
205
2.4
Fourier Series Expansion
Let ζ2 (a variable in the new coordinate system) be defined in a bounded domain, for example, [0, 2π). Then, any function w in the new coordinate system is periodic (with period 2π) with respect to ζ2 , and therefore, w can be expressed in terms of its Fourier series expansion, namely,
w=
l=∞ X
wl ejlζ2 ,
(18)
l=−∞
206 207 208 209 210 211 212
213
√ 1 R 2π −jlζ2 where j = −1 is the imaginary unit, ejlζ2 are the modes, and wl = 2π dζ2 0 we are the modal coefficients that are independent of variable ζ2 . For convenience, we define symbol Fl , such that when applied to a scalar-valued function w, it produces the l-th Fourier modal coefficient wl , and when applied to a vector (or matrix) it produces a vector (or matrix) of the same dimensions, with each of the components being equal to the l-th Fourier modal coefficient corresponding to the component of the original vector (or matrix). For example, Fl (u v w) = (Fl u Fl v Fl w) = (ul vl wl ) .
9
(19)
214
215
∂u , uD , σ N EW , fN EW , ∂ζ and hN EW , variational formulation (16) can be expressed as Using the Fourier series expansion representation for u,
Find * ∂v 216
217 218
219 220 221
∂ζ
1 u = Fl (u)ejlζ2 ∈ Fl (uD )ejlζ2 + HD (Ω) such that:
!
, Fp (σ N EW )Fl
227 228
E L2 (ΓN )
230 231 232
(20)
1 ∀v ∈ HD (Ω) .
1 u = Fl (u)ejlζ2 ∈ Fl (uD )ejlζ2 + HD (Ω) such that:
∂v ∂ζ
!
, Fk−l (σ N EW )Fl
∂u ∂ζ
!+
(21)
L2 (Ω2D )
= hFk (v) , Fk (fN EW )iL2 (Ω2D ) + hFk (v), Fk (hN EW )iL2 (ΓN (Ω2D )) 1 ∀Fk (v)ejkζ2 ∈ HD (Ω) .
In the above formula, for each equation k, we are employing Einstein’s summation convention for −∞ ≤ l ≤ ∞. However, if we employ the non-orthogonal coordinate system described in Section 2.2, and under the additional (realistic) assumption that σ is constant as a function of ζ2 , we note that Fk−l (σ N EW ) ≡ 0 if |k − l| > 2. Therefore, the infinite series in terms of l reduces for each k to a finite sum with at most five terms, namely l = k − 2, ..., k + 2. We obtain P P 1 Find u = l Fl (u)ejlζ2 ∈ l Fl (uD )ejlζ2 + HD (Ω) ! !+ * l=k+2 X ∂v ∂u , Fk−l (σ N EW )Fl Fk
229
L2 (Ω)
We select a mono-modal test function v = vk ejkζ2 , where the Fourier modal coefficient vk is a function of ζ1 and ζ3 . Then, variational problem (20) reduces by orthogonality of the Fourier modes in L2 to
226
L2 (Ω)
E
In the above formula, we are employing Einstein’s summation convention, with −∞ ≤ l, p ≤ ∞, and we are assuming that ΓN and ΓD are independent of ζ2 .
k
225
D
= v , Fl (fN EW )ejlζ2
D
222
224
+
+ v, Fl (hN EW )ejlζ2
Find * F
223
∂u j(l+p)ζ2 e ∂ζ
l=k−2 =
∂ζ
∂ζ
such that:
L2 (Ω2D )
(22)
hFk (v) , Fk (fN EW )iL2 (Ω2D ) + hFk (v), Fk (hN EW )iL2 (ΓN (Ω2D )) 1 (Ω) . ∀Fk (v)ejkζ2 ∈ HD
In order to implement variational problem (22) in a Finite Element code, we need to relate the Fourier series modal coefficients of the derivative of a function to the Fourier series modal coefficients of the function itself. In order 10
233
to establish this relation, we first note that ∂w ∂(Fl (w)ejlζ2 ) = Fl ∂ζ1 ∂ζ1
∂(Fl (w)ejlζ2 ) ∂w = Fl ∂ζ3 ∂ζ3
Fl
240 241 242 243 244 245 246 247 248
249 250 251 252 253 254
255 256
!
ejlζ2 .
∂w ∂ζ2
!
P
:= Fl
∂(
k
wk ejkζ2 ) ∂ζ2
!
!
= Fl
X
jkζ2
jkwk e
= jlFl (w) . (24)
k
Finally, by combining Equations (23) and (24), we obtain Fl
238
239
(23)
By invoking the Fourier series expansion of w, one obtains
236
237
ejlζ2 ,
∂(Fl (w)ejlζ2 ) = jlFl (w) ejlζ2 , ∂ζ2
234
235
!
2.5
∂w ∂ζ
!
=
∂(Fl (w)ejlζ2 ) −jlζ2 e . ∂ζ
(25)
A Self-Adaptive Goal-Oriented hp-FEM
Each of the above Fourier modal coefficients represents a 2D function in terms of variables ζ1 and ζ3 . Furthermore, variational problem (21) constitutes a system of linear equations in terms of 2D functions (Fourier modal coefficients). To solve the above system of linear equations, it is necessary to select a software capable of simulating 2D-DC problems. The choice of the 2D software is somehow arbitrary, since the corresponding Fourier series expansion in terms of the quasi-azimuthal component ζ2 is independent of the numerical algorithm employed to solve each 2D problem with respect to variables ζ1 and ζ3 . In this paper, we have selected as our starting point a 2D self-adaptive goaloriented hp-FEM. This goal-oriented hp-FEM delivers exponential convergence rates in terms of the error in the quantity of interest versus the number of unknowns and CPU time [19,20]; the outstanding performance of the hp-FEM for simulating diverse resistivity logging measurements has been documented in [2–4]. A description of the hp-FEM can be found in [17]. We refer to [16] for technical details on the goal-oriented adaptive algorithm.
11
257
258 259 260 261
3
Implementation
In this Section, we first assume that we have a software capable of solving 2D-DC problems for an arbitrary material conductivity tensor σ. Then, we describe the modifications that are necessary to simulate 3D borehole measurements acquired in deviated wells using the method described in Section 2. T
262 263 264 265 266 267 268 269 270
271 272 273 274
First, we compute the new material coefficient σ N EW := J −1 σJ −1 |J |, and all its Fourier modes σ N EW,i , i = −2, −1, 0, 1, 2. Analytic computation of these coefficients may be long and tedious. Therefore, we employ symbolic computations using Maple [21]. Accordingly, to include the final symbolic expressions into our FORTRAN code, we utilize the existing routine “Fortran” within the software Maple to produce Fortran source code. We only need to compute the positive coefficients, since negative Fourier modal coefficients are the complex conjugate of the corresponding positive Fourier modal coefficients, ¯ N EW,−i for every i, because σ N EW is a real-valued tensor. that is, σ N EW,i = σ Second, we define as many equations as number of Fourier modal coefficients we want to solve. This number may be modified during execution. Third, we need to modify the structure of the stiffness matrix to account for various equations (Fourier modes). To that end, we introduce the following notation: *
275
276 277
∂v ∂ζ
(k, k − l, l) := Fk
!
, Fk−l (σ N EW )Fl
∂u ∂ζ
!+
.
Then, according to Formula (21), we obtain the following structure for the stiffness matrix A for the example case of seven Fourier modes:
278
A=
0 0 0 0 (−3,0,−3) (−3,−1,−2) (−3,−2,−1) 0 0 0 (−2,1,−3) (−2,0,−2) (−2,−1,−1) (−2,−2,0) 0 0 (−1,2,−3) (−1,1,−2) (−1,0,−1) (−1,−1,0) (−1,−2,1) . 0 0 (0,2,−2) (0,1,−1) (0,0,0) (0,−1,1) (0,−2,2) 0 0 (1,2,−1) (1,1,0) (1,0,1) (1,−1,2) (1,−2,3) 0 0 0 (2,2,0) (2,1,1) (2,0,2) (2,−1,3) 0
279 280 281 282
(26)
L2 (Ω2D )
0
0
0
(3,2,1)
(3,1,2)
(27)
(3,0,3)
In the above matrix, rows and columns are associated with test and trial Fourier modal basis functions, respectively. We emphasize that the resulting stiffness matrix is, in general, penta-diagonal, since σ N EW,k−l = 0 for every |k − l| > 2. Furthermore, for subdomain III, we have σ N EW,k−l = 0 for every 12
283 284 285
286 287 288 289 290 291 292 293
294 295 296 297 298 299 300
301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317
318 319 320 321 322 323
|k − l| > 1, and thus, the stiffness matrix becomes tri-diagonal. For subdomain I, we have σ N EW,k−l = 0 if k − l 6= 0, and the corresponding stiffness matrix becomes simply diagonal. We note that it is possible to reproduce the above structure for the stiffness matrix at different levels, such as the degree-of-freedom (d.o.f.) level (several equations per d.o.f), the element stiffness matrix level (several element stiffness matrices per element), or the global stiffness matrix level (several global stiffness matrices). A different level selection entails a different ordering of the unknowns in the final global stiffness matrix, which may affect the performance of a direct solver. In order to simplify the implementation, we have used multiple equations structure at the d.o.f. level. The resulting system of linear equations needs to be solved with either a direct solver or an iterative solver. It seems natural to use an iterative solver based on a block-Jacobi preconditioner accelerated with a Krylov-based subspace optimization method, where the blocks of the preconditioner are defined by the 2D problems corresponding to the diagonal entries of matrix (27). However, in this paper we use a direct solver of linear equations to avoid additional numerical errors possibly introduced by the iterative solver of linear equations. To solve the linear system of equations, we use the sequential version of the multifrontal direct solver MUMPS (version 4.6.2) [22–24], with the ordering of the unknowns provided by METIS (version 4.0) [25]. The interface with the direct solver used in this work is based on the assembled stiffness matrix format. Notice that the element-by-element matrix interface assumes that element stiffness matrices are dense. In our case, and due to the sparse coupling structure of the Fourier modes—see matrix (27)—, element stiffness matrices are sparse. To illustrate the importance of this observation, we implemented both the element-by-element and the assembled stiffness matrix interfaces with MUMPS, and we compared results obtained with the same solver for a particular problem containing 6145 d.o.f. and 15 Fourier modes. For this problem, solver MUMPS interfaced with the assembled stiffness matrix format used approximately half the memory and was four times faster than when interfaced with the element-by-element stiffness matrix format. Specifically, in a machine equipped with a 2 Ghz AMD Opteron processor, solver MUMPS spent 88 seconds when using the assembled format vs. 364 seconds when using the element-by-element format. In the remainder of this section, we first discuss advantages and disadvantages of different gridding possibilities, and we describe the gridding techniques that we use in the work. Finally, we describe a sum factorization algorithm that is critical to significantly reducing the computational time during Gaussian integration of the stiffness matrix. The use of this integration technique (or of some other advanced integration method) is crucial for elements with high 13
327
polynomial order of approximation p. Otherwise, the time spent during integration may be larger than the total CPU time used by the remaining parts of the FEM code (including the LU factorization of the system of linear equations).
328
3.1
324 325 326
329 330 331 332 333
Gridding
We note that it is possible to employ different grids for each diagonal term of the stiffness matrix defined in (27). In this case, different spaces for test and trial functions will appear on various non-diagonal entries of the matrix. Nevertheless, to simplify the implementation, we employ a unique grid for all entries in matrix (27).
347
To construct an optimal hp-grid for each problem, we first manually select a coarse grid that is conforming to the geometry of the problem. Then, we employ a self-adaptive goal-oriented hp-FEM (see [16]). The self-adaptive algorithm utilizes a fine-grid solution to guide optimal mesh refinements. This fine-grid solution provides an error estimate function (not just a number) over coarser grids, and is used to perform optimal hp-refinements. The major limitation of this approach is that the computation of the fine-grid solution may be time and memory consuming. To minimize the latter computational cost, it is possible to use adaptivity over a few Fourier modes (perhaps only the central Fourier mode, or one Fourier mode at a time). However, the resulting hp-adapted grid constructed in this manner may provide inaccurate solutions. In Section 4.3, we analyze the total CPU time and the corresponding accuracy of the solution when different numbers of Fourier modes are considered during hp-adaptivity.
348
3.2
334 335 336 337 338 339 340 341 342 343 344 345 346
Sum Factorization
358
Construction of the stiffness matrix requires the integration of Equation (26) over 2D elements. For a higher-order quadrilateral finite element of degree p, it is customary to use a Gaussian integration rule of degree p + 1, which is exact for polynomials of degree 2p + 1. Thus, it provides exact integration if the material coefficients are constant (or linear) within the element when the element mapping is affine. Since the Jacobian matrix associated with the change of coordinates is not constant within each element, new material coefficient σ N EW is not constant within each element. Nevertheless, we shall still employ a Gaussian integration rule of degree p + 1 and simply disregard the integration error.
359
In addition, we accelerate the integration procedure using the sum factoriza-
349 350 351 352 353 354 355 356 357
14
360 361 362
363 364 365 366 367
tion algorithm described in [26]. The main advantage of a sum factorization algorithm for 2D integration is that the number of operations is reduced from p6 to p5 . Implementing a sum factorization algorithm requests the decomposition of each 2D shape function u = u2D as the product of two 1D shape functions: one in the horizontal direction (φζ1 ), and the second in the vertical direction (φζ3 ). A similar decomposition follows for the gradient operator, that is, ∇ζ u = (∇ζ1 φζ1 ) · (∇ζ3 φζ3 ), where we define
φζ3 0
∇ζ1 φζ1 := (φ0ζ1
368
jlφζ1
φζ1 ) ; ∇ζ3 φζ3 := 0
0
0
φζ3 0 0
φ0ζ3
.
(28)
375
After performing the above decomposition, Table B.1 describes the number of operations used by the sum factorization algorithm. In that table, symbol “X” indicates that a loop over the corresponding column is necessary to compute the row entry. Symbol “O” expresses the savings due to the sum factorization algorithm. No loop is necessary with respect to that column in order to compute the corresponding row entry, as opposed to standard integration, where symbol “O” should be replaced with symbol “X”.
376
[Table 1 about here.]
369 370 371 372 373 374
377
4
378
4.1
379 380 381 382 383 384 385 386
387 388
Numerical Results
Sources, Receivers, and Boundary Conditions (BCs)
4.1.0.1 Sources. It is customary in computer-aided simulations to model electrodes as Dirac’s delta functions. However, since the exact solution corresponding to a Dirac’s delta load has infinite energy, this load should not be used in combination with self-adaptive codes. Thus, we employ finite-size electrodes. Our model loop electrodes have a radius of 7 cm and a square cross-section of 2 cm × 2 cm. These dimensions are commensurate with those of logging instruments. We impose a constant f = ∇ · Jimp on the volume occupied by source electrodes.
4.1.0.2 Receivers. In the numerical simulations described in this paper, we assume that our logging instrument is equipped with one current (emitter) 15
389 390 391 392 393 394 395 396
397 398 399
400
401 402 403
electrode and two (or three) voltage (collector) electrodes. We note that realistic logging instruments typically incorporate a large number of electrodes (10-20) in order to perform several simultaneous measurements, and thus, enhance the focusing (vertical and horizontal resolution) of electrical currents while minimizing measurements errors. Despite the reduced number of electrodes used in our simulations, logging instruments assumed in this paper have been designed to reproduce the same basic physical principles of those customarily used in actual field operations. Depending upon the number of receivers, we define two different quantities of interest. The first one considers two measurement electrodes, and is defined as the difference of potential u between electrodes RX1 and RX2 , that is: L1 (u) =
|ΩRX1 |
Z
u dV −
ΩRX1
406 407 408 409 410
411 412 413 414 415
Z
|ΩRX2 |
u dV ,
(29)
ΩRX2
R
1 |ΩRX1 |
Z
u dV −
ΩRX1
+
405
1
where |ΩRXi | = ΩRX 1dV . The second one considers three electrodes, and i is defined as the second difference of potential u between three consecutive electrodes RX1 , RX2 , and RX3 , that is, L2 (u) =
404
1
2
Z
|ΩRX2 | 1
Z
|ΩRX3 |
u dV
ΩRX2
(30)
u dV .
ΩRX3
4.1.0.3 Boundary Conditions (BCs). A variety of BCs can be used to truncate the computational domain. For borehole geophysical applications, it is customary to use homogeneous Dirichlet BC in a sufficiently large spatial domain. This strategy is justified by the rapid decay of the electric potential in lossy media. Here, we follow the same approach and apply Dirichlet BCs in a sufficiently large domain. The remainder of this section is divided into three parts. The first part is intended to validate the simulation code. In the second part, we study the error due to the truncation of the Fourier series expansion in challenging problems. Finally, the third part draws physical conclusions from the simulations of practical 3D borehole measurements.
16
416
417 418 419
4.2
Validation of the Code
In this subsection, we select three model problems for which numerical solutions are already available, and compare them to those obtained with the software described in this paper.
430
We consider the three problems illustrated in Fig. B.3. Voltage electrodes are located 1 m and 1.25 m above the current electrode, respectively. Analytical solutions are available for all three problems (c.f., [27]). We consider the case of a deviated well. Because the formation material is unbounded, homogeneous, and isotropic, we know that the solution should be identical for all possible dip angles. In particular, it should coincide with the axial-symmetric solution (0 angle) that we compute with an existing 2D code that has already been verified [27] against existing analytical solutions. Using this high-accuracy 2D solution as the exact solution, we study the convergence of our method in terms of the relative error in the quantity of interest with respect to the number of Fourier modes used to construct the solution for different dip angles.
431
[Fig. 3 about here.]
420 421 422 423 424 425 426 427 428 429
432 433 434 435 436 437 438 439 440
441 442 443 444 445 446 447 448 449 450 451 452 453 454
Fig. B.4 displays the convergence history for each of the three problems described above versus the number of Fourier modes. The reference solution is computed with the 2D code for the axial-symmetric case (0 degrees). We observe that, for a sixty-degree deviated well, we obtain three (or more) significant digits (below 0.1% error) of the exact (2D) solution using only nine Fourier modes. For a thirty-degree deviated well, we obtain four (or more) significant digits (below 0.01% error) of the exact (2D) solution using only five Fourier modes. As expected, to achieve a given tolerance error an increase in dip angle also requires an increase in the number of Fourier modes. In Fig. B.4 we also observe that all curves are concave in a log-log scale, which indicates the exponential convergence of the method. This (exponentially) fast convergence is peculiar of spectral (higher-order) methods when applied to smooth solutions [28]. We note that the Fourier series expansion is a spectral method and that, in our applications, the solution on the quasiazimuthal direction is smooth because materials are assumed to be constant on the quasi-azimuthal direction and hence no geometrically-induced singularity occurs. Therefore, the exponential convergence displayed in Fig. B.4 is expected to hold for all our applications. In addition, our method has a fixed maximum bandwidth (up to five Fourier modes interact with each other), as opposed to traditional spectral methods, where the bandwidth grows unbounded as one adds more basis functions. In summary, our method combines the advantages of spectral methods (exponential convergence) while maintaining a low computational cost.
17
[Fig. 4 about here.]
455
456
4.3
Performance and Error Analysis
459
To perform an error analysis of our method, we consider two different types of logging instruments that are widely used by the logging industry: Throughcasing resistivity instruments, and Laterolog instruments.
460
4.3.1
457 458
461 462 463 464 465 466 467 468 469 470 471 472 473
474 475 476 477 478 479 480 481
482 483 484 485 486
Through-Casing Borehole Measurements
First, we select a through-casing model problem (see Fig.B.5). In this type of applications, a casing (thin metallic pipe) is placed against the wall of the borehole to avoid the mechanical collapse of the well. Since a metallic (steel) casing is a good electric conductor, electric currents travel long distances within the casing in the vertical direction. At the same time, a small amount of current leaks into the electrically conductive formation. This current leakage, which is several orders of magnitude smaller than the electric field itself, is proportional to the conductivity of the formation. Furthermore, Kaufmann demonstrated in [29] that the second vertical derivative of the electric potential 3 within the borehole is proportional to the current leaked into the formation, and therefore, also proportional to the formation conductivity. Based on this principle, logging measurements provide useful information about the conductivity of the subsurface formations penetrated by the steel-cased well. Computer simulations of through-casing resistivity measurements are very challenging because of three reasons: (1) It is necessary to consider a long computational domain in the vertical direction (sometimes, thousands of meters) or to employ a sophisticated truncation method—such as a Perfectly Matched Layer (PML) [30]—to account for currents within casing; (2) There exists a large contrast in material properties (typically, 9-14 orders of magnitude), and; (3) There exists a large dynamic range (quotient between the maximum value of the solution and the solution in the quantity of interest). These difficulties have discouraged the use of computer-aided simulations to analyze through-casing resistivity measurements in deviated wells. The only existing work in this area can be found in [13], where a 3D self-adaptive hpFEM was used to simulate the problem described in Fig. B.5. In that work, several days of CPU time were needed to perform simulations for 80 different 3
Through-casing resistivity logging instruments measure the second difference of the potential along the trajectory of the well given by Equation (30), which is an approximation of the second derivative of the potential along the trajectory of the well.
18
492
logging positions along the deviated well. With the method described in this paper, for the first time we perform accurate simulations of through-casing resistivity measurements in deviated wells using CPU times within one or two hours for 80 logging points (that is, 1-2 minutes per logging position). A performance and error analysis of our method based on the model problem described in Fig. B.5 follows below.
493
[Fig. 5 about here.]
487 488 489 490 491
494 495 496 497
For the hp self-adaptive goal-oriented refinement strategy, we select a tolerance error of 1% in the quantity of interest. That is, we request that the difference in the quantity of interest corresponding to the solutions in the final coarse and fine (globally hp-refined) grids, respectively, remain below 1%.
508
Executing the adaptive algorithm with many Fourier modes is computationally expensive. Thus, we shall restrict ourselves to calculations performed with only either one or five Fourier modes for the construction of adapted hp-grids. Thus, the final hp-grid may not be optimal due to the fact that we employ a few Fourier modes for its construction. Therefore, we shall consider the case of possibly p-enriching the final hp-grid for increased accuracy. Once the final hpgrid has been constructed, we shall employ a larger number of Fourier modes to compute the final solution in order to guarantee its accuracy. Specifically, we will use either nine or fifteen Fourier modes for the computation of the final solution. Table B.2 describes the eight possible algorithms we use to solve each problem.
509
[Table 2 about here.]
498 499 500 501 502 503 504 505 506 507
510 511 512 513 514 515 516 517 518 519 520
521 522
In Table B.3, we report the total time 4 associated with each of the eight algorithms defined in Table B.2, when considering a 200 m (vertical) × 100 m (horizontal) computational domain in a sixty-degree deviated well. CPU times correspond to the computation of a full log, consisting of 80 different logging positions. We observe a large difference in the CPU time as a function of the employed algorithm (case number): As we consider more Fourier modes and/or we globally p-enrich the final grid, we observe an increase of the CPU time. For a 2000 m (vertical) × 1000 m (horizontal) domain, the grids needed to achieve a similar accuracy contain more d.o.f., and the total CPU time increases by approximately 50% for cases I through IV, and by approximately 15% for cases V through VIII. We emphasize the importance of using fast integration rules and an adequate interface for the solver of linear equations, as described in Section 3. Both 4
CPU timings obtained with the Linux command “time” have been used to compute the total CPU times. The command “gprof” was used to assess the percentage of time spent on each subroutine.
19
533
the integration and the solution of the system of linear equations embody a significant portion (20% - 50%) of the the total CPU time, as indicated in Table B.3. When a standard integration routine is used, as opposed to the fast sum factorization integration technique, the CPU time spent during integration triplicates for the case of algorithm VIII, and thus, the overall CPU time of the entire code almost doubles. If an interface with the same solver (MUMPS) based on dense element-by-element stiffness matrices is used, then the solver time increases by a factor of four or five when considering fifteen Fourier modes, which severely deteriorates the overall performance of the code. The memory used by solver MUMPS also increases significantly as we consider more Fourier modes.
534
[Table 3 about here.]
523 524 525 526 527 528 529 530 531 532
548
Fig. B.6 (top panel) displays the simulated measurements corresponding to the through-casing resistivity problem described in Fig. B.5 in a sixty-degree deviated well, with the resistivity of casing equal to 10−5 Ω· m. Different curves correspond to the eight algorithms (case numbers) considered in this paper. We observe that the measured signal is (almost) proportional to the formation conductivity. Selecting the solution corresponding to ‘case VIII’ as the reference solution, Fig. B.6 (bottom panel) displays the relative error of the quantity of interest corresponding to the remaining seven cases. For most logging positions, the relative error is consistently below 10%. This level of accuracy is enough from an engineering point of view to properly assess the resistivity of the formation. Nonetheless, substantially larger errors (above 100%) are observed in the highly resistive layer when employing algorithms I and III, which indicates the accuracy limitations of those algorithms. Algorithms II and IV provide a good compromise between accuracy and CPU time.
549
[Fig. 6 about here.]
535 536 537 538 539 540 541 542 543 544 545 546 547
550 551 552 553 554
555 556 557 558 559 560 561 562
Fig. B.7 displays similar results to those shown in Fig. B.6, when considering a larger computational domain, specifically, a domain of size 2000 m (vertical) × 1000 m (horizontal). These results further illustrate the inaccurate solutions delivered by algorithms I and III. Again, algorithms II and IV provide the best compromise between accuracy and CPU time. From the results shown in Fig. B.7 (bottom-left panel), we emphasize the discrepancy between solutions obtained with this method and the solution obtained with the 3D hp-FEM utilized in [13]. We note that at the points where 0.3 m < z < 0.7 m (points of discrepancy between 3D and 2D solutions), the 3D hp-FE solution did not converge at the desired level of accuracy (3%), and the simulation was stopped after various days of computations. With the new method presented in this paper, we have reduced the computational time from several days to less than two hours (if we employ, for example, 20
563 564 565 566 567
algorithm IV). Also, with the new method we obtain (almost) identical results as we increase the number of Fourier modes and/or enrich the final hp-grid. This consistency of results among various gridding algorithms and number of Fourier modes indicates that the discretization errors are small, and therefore, that the solutions obtained with the new method are accurate.
571
The above observations confirm that the new method is substantially more accurate than the 3D hp-FEM. In summary, with the new method we dramatically reduce the computational time while we gain accuracy on the final solution.
572
[Fig. 7 about here.]
568 569 570
582
Finally, Fig. B.8 analyzes the numerical error due to the truncation of the computational domain. We compare the solution obtained with a 200 m (vertical) × 100 m (horizontal) domain against the solution obtained with a 2000 m (vertical) × 1000 m (horizontal) domain. A relative difference below 20% is obtained at all logging points. Furthermore, with the exception of a few logging points located across the most resistive layer of the formation, the discrepancy between both solutions remains below 10%. These differences can be neglected from the engineering point of view in the case of through-casing resistivity measurements, since the uncertainty of actual field measurements is often above the observed level of discrepancy.
583
[Fig. 8 about here.]
573 574 575 576 577 578 579 580 581
584
585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600
4.3.2
Laterolog Resistivity Measurements
We now consider Laterolog-type resistivity measurements acquired at zero frequency (DC). This type of measurements are widely used by the logging industry when the borehole mud is more electrically conductive than the surrounding formation. Fig. B.9 describes the corresponding logging instrument, electrodes and materials. The current (emitter) electrode is excited by prescribing a constant flux with total current equal to 1 A, that is, f = ∇ · Jimp is equal to one over the volume of the current (emitter) electrode. For the sake of simplicity, we avoid the use of voltage sources (prescribing the electric potential at the source) and bucking electrodes (used to maintain a zero difference of potential between two electrodes). We note, however, that voltage sources and bucking electrodes may enhance the dependence of the measurements upon the formation resistivity, and therefore, facilitate a fast inversion of the measurements. They can be easily modeled using a non-homogeneous Dirichlet BC (for the voltage source) and a slight modification of the resulting system of linear equations as described, for example, in [31] (for modeling bucking electrodes).
21
[Fig. 9 about here.]
601
612
Table B.4 summarizes the CPU time spent by each of the eight algorithms defined in Table B.2 for the case of a sixty-degree deviated well. CPU times correspond to the computation of a full log, consisting of 80 different logging positions. Each logging position is separated by a true vertical distance of 5 cm. In similar fashion to the results summarized in Table B.3, we observe large differences in the CPU times as a function of the employed algorithm (case number). However, we need only roughly half of the time for simulating Laterolog resistivity measurements compared to that needed to simulate through-casing resistivity measurements. This behavior occurs because fewer unknowns are necessary to simulate Laterolog resistivity measurements to achieve a similar level of accuracy.
613
[Table 4 about here.]
602 603 604 605 606 607 608 609 610 611
614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630
Fig. B.10 (left panel) displays the simulated measurements corresponding to the Laterolog resistivity problem described in Fig. B.9 in a sixty-degree deviated well. Different curves correspond to the eight algorithms (case numbers) considered in this paper. We observe a strong correlation between the simulated signal and the formation conductivity. Numerically, we observe that curves obtained using algorithms II,IV,V,VI,VII, and VIII are (almost) identical. By selecting the solution corresponding to ‘case VIII’ as reference, Fig. B.10 (right panel) displays the relative error of the quantity of interest corresponding to the remaining seven cases. At most logging positions, the relative error remains below 2% for algorithms II,IV,V,VI, and VII. This level of accuracy is exceptionally good from the engineering point of view when assessing the resistivity of the formation. Nonetheless, substantially larger errors (above 60%) are observed at various logging points when employing algorithms I and III, as it was also the case with through-casing resistivity measurements. Again, we observe the accuracy limitations associated with these two algorithms: For Laterolog instruments, algorithms II and IV seem to also provide the best compromise between accuracy and CPU time.
632
Algorithm IV, which only utilizes one Fourier mode for the adaptivity, is used to simulate resistivity measurements in the remainder of this paper.
633
[Fig. 10 about here.]
631
634
635 636 637
4.4
Well-Logging Applications
In this section, we illustrate the applicability of this method of solution by studying, for example, the effect of water invasion in through-casing resistivity measurements for various dip angles. We present the first results, which 22
638 639 640
are of great interest to the logging industry, and provide a clear physical interpretation of the water invasion effect for different layers in the formation and for various dip angles.
654
We consider the through-casing resistivity problem described in Fig. B.5 and a computational domain of 2000 m (vertical) × 1000 m (horizontal). Fig. B.11 compares simulation results for various dip angles corresponding to casing with resistivity equal to 10−5 Ω· m (left panel), and 2.3 × 10−7 Ω· m (right panel), respectively. We observe that measurements simulated in highly deviated wells are proportional to the average of the conductivity of formation materials surrounding the receivers. As indicated in Fig. B.11, the minimum and maximum recorded measurements increase as we decrease the dip angle. As a function of the casing resistivity, we observe a dramatic change of the intensity of the received signal, as physically expected. However, we observe qualitatively similar results for different values of casing resistivity. This behavior is consistent with the result for vertical wells provided by [29], where the author indicates that simulation results as a function of casing resistivity should be qualitatively identical.
655
[Fig. 11 about here.]
641 642 643 644 645 646 647 648 649 650 651 652 653
656 657 658
659 660 661 662 663 664 665 666
667 668 669 670 671 672 673 674 675 676 677
In the remainder of this paper, we assume a casing resistivity equal to 2.3 × 10−7 Ω· m, since this is a typical value of casing conductivity encountered in actual field applications. We now consider the effect of water-base mud-filtrate invasion in the two layers of resistivities equal to 0.01 Ω· m (layer 1) and 10000 Ω· m (layer 2), respectively. In so doing, we assume a piston-like radial water invasion with radial length of invasion equal to 10 cm and 50 cm, respectively. For the invaded part of the conductive layer (layer 1), we assume that the layer resistivity after being invaded with water is equal to 1 Ω· m. For the invaded part of the resistive layer (layer 2), we assume that the resulting resistivity is equal to 10 Ω· m. Fig. B.12 displays simulation results for the case of invading the conductive layer with water at different dip angles. The effect of water invasion on the simulated measurements is qualitatively similar for all dip angles. We observe a strong effect on the results due to water injection. With only 10 cm of radial length of water invasion, the simulated measurements decrease by approximately one order of magnitude. A larger effect ensues when the radius of invasion increases to 50 cm. Further increase of the radius of water invasion hardly affects the measurements, since the radial length of investigation of these logging instruments is relatively short (10-70 cm). We also observe that the effect of water invasion into layer 1 is non-local, as it affects measurements above and below the layer where water invasion takes place. 23
684
Fig. B.13 displays simulation results for the case of invading the resistive layer with water for different dip angles. We observe that the effect of water invasion on the simulated measurements is qualitatively similar for all dip angles. However, in this case we observe that the effect of water invasion into layer 2 is local, and that it only affects the measurements within the resistive layer. As physically expected, a larger measurement variation ensues when the radius of invasion increases.
685
[Fig. 12 about here.]
686
[Fig. 13 about here.]
678 679 680 681 682 683
687
688 689
690 691 692 693 694 695 696
697 698 699 700 701 702 703 704
705 706 707 708 709 710 711 712
5
Discussion About Further Applications
The method presented in this paper is efficient, reliable, and accurate in various engineering applications, for two reasons: • Material properties in the newly defined non-orthogonal system of coordinates are constant in the quasi-azimuthal direction, and thus, one Fourier mode is sufficient to reproduce exactly the material properties (in our case, the conductivity σ). • The Jacobian matrix expressing the change of coordinates from Cartesian to the non-orthogonal system of coordinates can be represented exactly with only five Fourier modes. Thus, the use of this method is limited by the geometrical description of the problem. Arbitrary 3D geometries will, in general, not satisfy the two properties described above. However, the method is physics independent, and can be applied to simulate diverse measurements in deviated wells, such as those based on electrodynamics, sonic propagation (acoustics coupled with elasticity), flow in porous media, geomechanics, etc. An application of this method to electrodynamic measurements shall be presented in a forthcoming paper. Despite the geometrical limitations inherent to our method, there exists another particular geometry of interest for logging measurements for which our formulation is also suitable: measurements acquired with the logging instrument eccentered from the axis of the well. Appendix B describes a suitable change of variables for borehole-eccentered measurements. Thus, it is possible to simulate borehole-eccentered measurements acquired in deviated wells by composing the change of coordinates for deviated wells with the change of coordinates for borehole-eccentered measurements. 24
713 714 715 716 717 718
The presented method may also be used in combination with Perfectly Matched Layers (PML) to truncate the computational domain. Indeed, the interpretation of PML in terms of a change of variables in the complex plane (described in [32]), makes the implementation of a PML trivial by simply composing the change of variables used in our method with the change of variables pertaining to the PML implementation.
731
In addition, the method described in this paper is ideal for the solution of inverse problems. Since the dip angle of the well is often measured by the logging instrument, and therefore, known a priori, only one Fourier mode is needed to reproduce exactly the material properties. In other words, material properties are constant with respect to the quasi-azimuthal variable, whereupon the inverse problem becomes a 2D problem. Using different grids for the forward and inverse problems (as proposed in [33]), we realize that only a 2D grid with just one Fourier mode is needed for reproducing the inverse solution (material coefficients). This observation about the dimensionality of the inverse problem greatly reduces the computational effort needed to solve it. We firmly believe that the method proposed in this paper will have a great practical impact in the logging industry, as it allows accurate and inexpensive simulations of forward and inverse borehole problems.
732
6
719 720 721 722 723 724 725 726 727 728 729 730
733 734 735 736 737 738 739 740
741 742 743 744 745 746 747
748 749 750
Conclusions
We have introduced and successfully tested a new simulation method based on the use of a non-orthogonal system of coordinates with a Fourier series expansion in one direction. The method is suitable for the simulation of borehole resistivity logging measurements acquired in deviated wells. For these geometries, material coefficients are constant in the new system of coordinates, and only five Fourier modes are necessary to reproduce exactly the new materials constructed by incorporating the change of coordinates. The new method is suitable for solving forward and inverse problems. The implementation of the new method of solution is based on the superposition of 2D algorithms. Its implementation requires a fraction of the time needed to develop a conventional 3D simulator. In order to achieve efficient computer performance, special care was taken during integration, where a sum factorization algorithm is employed. Also, it is essential to use an adequate solver (or solver interface) that takes advantage of the sparsity of the ensuing finite-element matrices. The new solution method delivers exponential convergence rates in terms of the error in the quantity of interest versus the number of Fourier modes. In addition, the bandwidth of the corresponding system of linear equations 25
751 752
remains bounded (each Fourier mode only interacts with no more than five Fourier modes).
761
We have validated the method, and illustrated its efficiency by solving various forward problems based on borehole electrostatic measurements. Results indicate that accurate solutions are obtained using only a limited number of Fourier modes for the solution (typically, below 10), thereby enabling a significant complexity reduction. Specifically, for through-casing borehole resistivity measurements, the computational time was dramatically reduced from several days (when using a 3D-adaptive hp-FEM code) to less than two hours (when using the new method). In addition, the consistency and reliability of the results indicates that we also gain accuracy.
762
ACKNOWLEDGMENTS
753 754 755 756 757 758 759 760
768
This work was financially supported by The University of Texas at Austin’s Joint Industry Research Consortium on Formation Evaluation sponsored by Anadarko, Aramco, Baker Atlas, British Gas, BHP-Billiton, BP, Chevron, ConocoPhillips, ENI E&P, ExxonMobil, Halliburton, Marathon, Mexican Institute for Petroleum, Hydro, Occidental Petroleum, Petrobras, Schlumberger, Shell E&P, Statoil, TOTAL, and Weatherford International, Ltd.
769
A
763 764 765 766 767
770
771
772
Fourier Series Expansion of the Metric Associated with the NonOrthogonal System of Coordinates for Deviated Wells
2π 1 Z G nm e−jkζ2 dζ2 of tensor All non-zero Fourier modal coefficients G k = 2π 0 matrix G with respect to variable ζ2 are given by
G0 =
1 + 0
0.5(θ0 f10 )2
0
0
0
ζ12 + 0.5(θ0 f1 )2
, 0
0
0
G1 =
0
0.5θ0 f10
1
0
0.5θ0 f10
0
0.5jθ0 f1
0.5jθ0 f1
0
26
(A.1)
,
(A.2)
0 2 0.25(θ0 f1 )
0.25jθ02 f1 f10
G 2 = 0.25jθ02 f1 f10
0 0
2
,
(A.3)
G −1 = G 1 , and G −2 = G 2 .
(A.4) (A.5)
−0.25(θ0 f1 )
0
0
0
−1
2π 1 Z nm −jkζ2 )k = G e dζ2 of the 2π
773
All non-zero Fourier modal coefficients (G
774
inverse tensor matrix G −1 with respect to variable ζ2 are given by
0
(G
−1
1 0
)0 =
0
0
0
1/ζ12
0
0
1 + 0.5θ02 ((f1 /ζ1 )2 + (f10 )2 )
,
(A.6)
,
(A.7)
,
(A.8)
(G −1 )−1 = (G −1 )1 , and
(A.9)
0
(G
−1
)1 =
0
−0.5θ0 f10
−0.5θ0 f10
0
0
−0.5jθ0 f1 /ζ12
−0.5jθ0 f1 /ζ12
0
(G
−1
)2 =
0 0
0
0
0
0
0
0
0.25θ02 (−(f1 /ζ1 )2 + (f10 )2 ) (G
−1
)−2 = (G
−1
)2 .
(A.10)
779
REMARK: We note that G 0 6= diag(1, ζ12 , 1) when θ0 6= 0. This fact implies that the axial-symmetric formulation is not the optimal 2D formulation for approximating results in deviated wells. Furthermore, the optimal 2D formulation (in the Fourier sense) for approximating results in deviated wells stems from the approximation G = G 0 .
780
B
775 776 777 778
781
782 783
Non-Orthogonal Coordinate System for Borehole-Eccentered Logging Instruments
For borehole-eccentered logging instruments, as the one described in Fig. B.14 (left panel), we introduce the following non-orthogonal coordinate system ζ = 27
784
(ζ1 , ζ2 , ζ3 ) in terms of the Cartesian coordinate system x = (x1 , x2 , x3 ): x1
785
x2
= f1 (ζ1 ) + ζ1 cos ζ2 ,
= ζ1 sin ζ2
(B.1)
x3 = ζ3
786
787
where f1 is defined by the formula
f1 (ζ1 ) = f1 =
ρ0 ζ
1
ρ1
ζ1 < ρ1 − ρ2 ρ0 − ρ2
0
788
(B.2)
ζ1 > ρ2
The corresponding derivative is given by 0
789
ρ1 ≤ ζ1 ≤ ρ2 .
ζ1 < ρ1
ρ0 f10 (ζ1 ) = f10 = ρ1 − ρ2 0
ρ1 < ζ1 < ρ2 .
(B.3)
ζ1 > ρ2
797
Intuitively, ρ1 is defined such that ζ1 < ρ1 covers the area occupied by the eccentered logging instrument, which corresponds to the interior part of the black circle shown in Fig. B.14 (left panel). Outside the borehole, which is identified by the dotted circle shown in Fig. B.14 (left panel), we employ a cylindrical coordinate system. Finally, the area between the logging instrument and the borehole wall is intended to “glue” all subdomains so that the resulting system of coordinates is globally continuous, bijective, and with positive Jacobian.
798
[Fig. 14 about here.]
790 791 792 793 794 795 796
799 800
The Jacobian matrix associated with the above change of coordinates is given by: (
801
J =
∂xi ∂ζj
0 f1
)
= i,j=1,2,3
+ cos ζ2
sin ζ2
0
802 803
−ζ1 sin ζ2
0
ζ1 cos ζ2
. 0
0
1
(B.4)
Accordingly, the determinant of the Jacobian |J | is equal to |J | = ζ1 [1 + f10 cos(ζ2 )] > 0 if ρ0 + ρ1 < ρ2 . 28
804 805 806
For the case of eccentered measurements, the described change of coordinates for borehole-eccentered measurements has three essential properties that make it suitable and attractive for numerical simulations:
811
• It is globally continuous, bijective, and with positive Jacobian. • Material properties are constant with respect to the quasi-azimuthal direction ζ2 . • Only a few Fourier modes in terms of ζ2 are necessary to approximate the tensor metric and its inverse.
812
References
807 808 809 810
813 814 815
816 817 818 819
820 821 822 823
824 825 826
827 828 829
830 831
832 833 834
835 836 837
838 839 840
[1] X. Lu, D. L. Alumbaugh, One-dimensional inversion of three-component induction logging in anisotropic media, SEG Expanded Abstract 20 (2001) 376– 380. [2] D. Pardo, L. Demkowicz, C. Torres-Verdin, M. Paszynski, Simulation of resistivity logging-while-drilling (LWD) measurements using a self-adaptive goal-oriented hp-finite element method, SIAM Journal on Applied Mathematics 66 (2006) 2085–2106. [3] D. Pardo, C. Torres-Verdin, L. Demkowicz, Simulation of multi-frequency borehole resistivity measurements through metal casing using a goal-oriented hp-finite element method, IEEE Transactions on Geosciences and Remote Sensing 44 (2006) 2125–2135. [4] D. Pardo, C. Torres-Verdin, L. Demkowicz, Feasibility study for twodimensional frequency dependent electromagnetic sensing through casing, Geophysics 72 (2007) F111–F118. [5] L. Tabarovsky, M. Goldman, M. Rabinovich, K. Strack, 2.5D modeling in electromagnetic methods of geophysics, Journal of Applied Geophysics 35 (1996) 261–284. [6] J. Zhang, R. L. Mackie, T. R. Madden, 3-D resistivity forward modeling and inversion using conjugate gradients, Geophysics 60 (1995) 1312–1325. [7] V. L. Druskin, L. A. Knizhnerman, P. Lee, New spectral Lanczos decomposition method for induction modeling in arbitrary 3-D geometry, Geophysics 64 (3) (1999) 701–706. [8] G. A. Newman, D. L. Alumbaugh, Three-dimensional induction logging problems, part 2: A finite-difference solution, Geophysics 67 (2) (2002) 484– 491. [9] S. Davydycheva, V. Druskin, T. Habashy, An efficient finite-difference scheme for electromagnetic logging in 3D anisotropic inhomogeneous media, Geophysics 68 (5) (2003) 1525–1536.
29
841 842
843 844
845 846 847
848 849 850 851
852 853 854
855 856
857 858 859 860
861 862
863 864 865
866 867 868
869 870 871
872 873
874 875 876
877 878 879
[10] T. Wang, S. Fang, 3-D electromagnetic anisotropy modeling using finite differences, Geophysics 66 (5) (2001) 1386–1398. [11] T. Wang, J. Signorelli, Finite-difference modeling of electromagnetic tool response for logging while drilling, Geophysics 69 (1) (2004) 152–160. [12] D. B. Avdeev, A. V. Kuvshinov, O. V. Pankratov, G. A. Newman, Threedimensional induction logging problems, part 1: An integral equation solution and model comparisons, Geophysics 67 (2002) 413–426. [13] D. Pardo, C. Torres-Verdin, M. Paszynski, Simulation of 3D DC borehole resistivity measurements with a goal-oriented hp finite element method. Part II: Through casing resistivity instruments, Computational Geosciences, in press. Preprint available at: www.ices.utexas.edu/%7Epardo. [14] A. Abubakar, P. M. van den Berg, Nonlinear inversion in electrode logging in a highly deviated formation with invasion ising an oblique coordinate system, IEEE Transactions on Geoscience and Remote Sensing 38 (2000) 25–38. [15] R. Aris, Vectors, Tensors, and the Basic Equations of Fluid Mechanics, PrenticeHall, 1962. [16] D. Pardo, L. Demkowicz, C. Torres-Verdin, L. Tabarovsky, A goal-oriented hp-adaptive finite element method with electromagnetic applications. Part I: electrostatics, International Journal for Numerical Methods in Engineering 65 (2006) 1269–1309. [17] L. Demkowicz, Computing with hp-Adaptive Finite Elements. Volume I: One and Two Dimensional Elliptic and Maxwell Problems, Chapman and Hall, 2006. [18] A. Ward, J. B. Pendry, Calculating photonic Green’s functions using a nonorthogonal finite-difference time-domain method, Phys. Rev. B 58 (1998) 7252–9. [19] I. Babuska, B. Guo, Approximation properties of the h-p version of the finite element method, Computer Methods in Applied Mechanics and Engineering 133 (3-4) (1996) 319–346. [20] C. Schwab, p- and hp-finite element methods, Numerical Mathematics and Scientific Computation, The Clarendon Press Oxford University Press, New York, 1998, theory and Applications in Solid and Fluid Mechanics. [21] www.maplesoft.com, MAPLE: Math Software for Engineers, Educators and Students (2007). [22] P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, Multifrontal parallel distributed symmetric and unsymmetric solvers, Computer Methods in Applied Mechanics and Engineering 184 (2000) 501–520. [23] P. R. Amestoy, I. S. Duff, J. Koster, J.-Y. L’Excellent, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM Journal of Matrix Analysis and Applications 23 (1) (2001) 15–41.
30
880 881 882
883 884
885 886 887
888 889 890
891 892
893 894
895 896 897 898
899 900 901
902 903 904
905 906 907 908
[24] P. R. Amestoy, A. Guermouche, J.-Y. L’Excellent, S. Pralet, Hybrid scheduling for the parallel solution of linear systems, Parallel Computing 32 (2006) 136– 156. [25] www-users.cs.umn.edu/ karypis/metis, Partitioning Algorithms (2007).
METIS
-
Family
of
Multilevel
[26] J. Kurtz, Fully automatic hp-adaptivity for acoustic and electromagnetic scattering in three dimensions, Ph.D. thesis, The University of Texas at Austin (May 2007). [27] M. Paszynski, L. Demkowicz, D. Pardo, Verification of goal-oriented hpadaptivity, Computers and Mathematics with Applications 50 (2005) 1395– 1404. [28] J. P. Boyd, Chebyshev and Fourier Spectral Methods, Second Edition, DOVER Publications, Inc., 2000. [29] A. A. Kaufman, The electrical field in a borehole with casing, Geophysics 55 (1) (1990) 29–38. [30] D. Pardo, L. Demkowicz, C. Torres-Verdin, C. Michler, PML enhanced with a self-adaptive goal-oriented hp finite-element method and applications to through-casing borehole resistivity measurements, Submitted to: SIAM Journal on Scientific Computing. Preprint available at: www.ices.utexas.edu/%7Epardo. [31] B. I. Anderson, Modeling and Inversion Methods for the Interpretation of Resistivity Logging Tool Response, Ph.D. thesis, Delft University of Technology (2001). [32] W. C. Chew, W. H. Weedon, A 3D perfectly matched medium from modified Maxwell’s equations with streched coordinates, Microwave Opt. Tech. Lett. 7 (1994) 599–604. [33] W. Bangerth, A framework for the adaptive finite element solution of large inverse problems. I. Basic techniques, Tech. Rep. 2004-39, Institute for Computational Engineering and Sciences, The University of Texas at Austin (2004).
31
909
910 911 912 913
914 915 916 917 918 919 920 921 922 923 924 925 926
927 928 929 930
931 932 933 934 935
936 937 938 939 940 941 942
List of Figures
B.1 3D description of the geometry of a resistivity logging instrument operating in a deviated well. We display several materials in the formation as well as the mud-filtrate invasion effect occurring in the proximity of the well.
35
B.2 Cross section showing the 3D geometry of a resistivity logging instrument in a vertical well penetrating dipping layers. Oblique circles in the left panel indicate the “quasi-azimuthal” direction ζ2 in a non-orthogonal system of coordinates. (x1 , x2 , x3 ) represents the Cartesian system of coordinates, and (ζ1 , ζ2 , ζ3 ) represents the new non-orthogonal system of coordinates. The new system of coordinates is different in each of the three sub-domains (right panel). Subdomain I corresponds to the logging instrument, subdomain II to the borehole, and subdomain III to the formation. The new system of coordinates is globally continuous, as indicated by the parameterization. Symbol ‘O’ indicates the origin of both systems of coordinates.
36
B.3 Description of three different simulation problems in a 100 m × 200 m computational domain, composed of (possibly) a mandrel—Material I—, a borehole—Material II—, and a uniform material in the formation—Material III—.
37
B.4 Convergence history as a function of the number of Fourier modes included in the numerical simulation. Different curves correspond to the three problems described in Fig. B.3 for two dip angles: 30 degrees—solid curves—, and 60 degrees—dashed curves—, respectively.
38
B.5 2D cross-section of an axial-symmetric through-casing resistivity problem in a borehole environment. Measurements consist of one current electrode (emitter) and three voltage electrodes (collectors). The subsurface earth formation is assumed to be composed of four different horizontal layers with varying resistivities in the range from 0.01 Ω· m to 10000 Ω· m.
39
32
943 944 945 946 947 948 949
950 951 952 953 954 955 956 957
958 959 960 961 962 963
964 965 966 967 968 969
970 971 972 973 974 975
976 977 978 979 980 981
B.6 Simulated through-casing resistivity measurements in a sixty-degree deviated well. Size of computational domain: 200 m (vertical) × 100 m (horizontal). Different curves correspond to different algorithms summarized in Table B.2. Top panel: solution in the quantity of interest as a function of logging position. Bottom panel: Relative error with respect to the reference solution corresponding to algorithm (case) VIII.
40
B.7 Simulated through-casing resistivity measurements in a sixty-degree deviated well. Size of computational domain: 2000 m (vertical) × 1000 m (horizontal). Different curves correspond to different algorithms summarized in Table B.2. Top-left, top-right, and bottom-left panels: solution in the quantity of interest as a function of logging position. Bottom-right panel: Relative error with respect to the reference solution corresponding to algorithm (case) VIII.
41
B.8 Simulated through-casing resistivity measurements in a sixty-degree deviated well. Different curves correspond to difference sizes of computational domain. Left panel: simulated measurements. Right panel: Relative error using the solution on the largest spatial domain (2000 m × 1000 m) as the reference solution.
42
B.9 2D cross-section of an axial-symmetric Laterolog resistivity problem in a borehole environment. Measurements are based on one current (emitter) and two voltage (collector) electrodes. The subsurface earth formation is assumed to be composed of four different horizontal layers with varying resistivities, from 0.5 Ω· m to 20 Ω· m.
43
B.10 Simulated Laterolog resistivity measurements acquired in a sixty-degree deviated well. Different curves correspond to different versions of the simulation algorithm. Left panel: solution in the quantity of interest as a function of logging position. Right panel: Relative error with respect to the reference solution corresponding to algorithm (case) VIII.
44
B.11 Simulated through-casing resistivity measurements. Size of computational domain: 2000 m (vertical) × 1000 m (horizontal). Different curves correspond to different dip angles: 0, 30, 45, and 60 degrees. Left panel: casing resistivity equal to 10−5 Ω· m. Right panel: casing resistivity equal to 2.3 × 10−7 Ω· m.
45
33
982 983 984 985 986 987 988 989
990 991 992 993 994 995 996 997
998 999 1000 1001 1002 1003 1004 1005
B.12 Simulated through-casing resistivity measurements. Casing resistivity equal to 2.3 × 10−7 Ω· m. Size of computational domain: 2000 m (vertical) × 1000 m (horizontal). Different panels correspond to different dip angles: 0 (top-left), 30 (top-right), 45 (bottom-left), and 60 degrees (bottom-right). Different curves correspond to different radial lengths of water invasion into the conductive layer (layer 1): no invasion, 10 cm invasion, and 50 cm invasion.
46
B.13 Simulated through-casing resistivity measurements. Casing resistivity equal to 2.3 × 10−7 Ω· m. Size of computational domain: 2000 m (vertical) × 1000 m (horizontal). Different panels correspond to different dip angles: 0 (top-left), 30 (top-right), 45 (bottom-left), and 60 degrees (bottom-right). Different curves correspond to different radial lengths of water invasion into the resistive layer (layer 2): no invasion, 10 cm invasion, and 50 cm invasion.
47
B.14 Left panel: Top view of the geometry describing the location of a resistivity tool eccentered in the borehole. The radius of the logging instrument is equal to ρ1 , and the radius of the borehole is equal to ρ2 , with the distance between the center of the logging instrument and the center of the borehole equal to ρ0 . Right panel: Iso-lines corresponding to the change of coordinates for borehole-eccentered measurements described in Formula B.1, with ρ0 = 0.3, ρ1 = 0.5, and ρ2 = 1.
48
34
Borehole Invasion Effect
-
*
Electrodes Logging Instrument
-
-
Fig. B.1. 3D description of the geometry of a resistivity logging instrument operating in a deviated well. We display several materials in the formation as well as the mud-filtrate invasion effect occurring in the proximity of the well.
35
ζ3
x3 x2 6 c- x1
ζ2
x3 x2 6 c- x1
6ζ1
O SU BD OM
ζ3
ζ3
6 ζ1 QQ s
6 ζ1 Q s Q
AI
N
III
SUBDOMAIN II SUBDOMAIN I SUBDOMAIN II
PA
ζ2
RA M
ET
ER
SU
BD
IZ
OM
AT IO
N
AI
N
III
Fig. B.2. Cross section showing the 3D geometry of a resistivity logging instrument in a vertical well penetrating dipping layers. Oblique circles in the left panel indicate the “quasi-azimuthal” direction ζ2 in a non-orthogonal system of coordinates. (x1 , x2 , x3 ) represents the Cartesian system of coordinates, and (ζ1 , ζ2 , ζ3 ) represents the new non-orthogonal system of coordinates. The new system of coordinates is different in each of the three sub-domains (right panel). Subdomain I corresponds to the logging instrument, subdomain II to the borehole, and subdomain III to the formation. The new system of coordinates is globally continuous, as indicated by the parameterization. Symbol ‘O’ indicates the origin of both systems of coordinates.
36
t
Rxt2
Problem 1: Uniform material Material I: 1 Ω· m Material II: 1 Ω· m Material III: 1 Ω· m
0.25 m
Rx1
200 m
1m t
Problem 2: LWD type Material I: 10−5 Ω· m Material II: 10 Ω· m Material III: 1 Ω· m
Material III
Material II 0.05 m
0.05 m
Material I
Tx
Problem 3: Laterolog type Material I: 104 Ω· m Material II: 0.2 Ω· m Material III: 1 Ω· m
99.9 m
Fig. B.3. Description of three different simulation problems in a 100 m × 200 m computational domain, composed of (possibly) a mandrel—Material I—, a borehole—Material II—, and a uniform material in the formation—Material III—.
37
Fig. B.4. Convergence history as a function of the number of Fourier modes included in the numerical simulation. Different curves correspond to the three problems described in Fig. B.3 for two dip angles: 30 degrees—solid curves—, and 60 degrees—dashed curves—, respectively.
38
5 Ω· m t
100-1000 m
Casing resistivity: 10−5 − 10−7 Ω· m Casing thickness: 0.0127 m
0.25 m Rx3 t
0.25 m Rx2 Rx1
0.1 m
0.25 m
5 Ω· m
Tx
0.01 Ω· m
5 Ω· m
100-1000 m
t
Metal casing
1m
10000 Ω· m
0.5 m
t
100-1000 m
Fig. B.5. 2D cross-section of an axial-symmetric through-casing resistivity problem in a borehole environment. Measurements consist of one current electrode (emitter) and three voltage electrodes (collectors). The subsurface earth formation is assumed to be composed of four different horizontal layers with varying resistivities in the range from 0.01 Ω· m to 10000 Ω· m.
39
5 Ω· m
0.5 0.25
10000 Ω· m
0 0.01 Ω· m
−0.25 −0.5 −0.75
5 Ω· m
−1
−8
−6
1
0.5
5 Ω· m
10000 Ω· m 0.01 Ω· m
−0.75 −1
5 Ω· m −2
10
0
10 Relative error (in %)
10000 Ω· m
0 0.01 Ω· m
−0.25 −0.5 −0.75
5 Ω· m
−8
−6
1
0
−0.5
0.25
−4
10 10 Second difference of potential (V)
0.25
−0.25
Case V Case VI Case VII Case VIII
5 Ω· m
0.5
10
Case I Case II Case III Case IV
0.75
0.75
−1
−4
10 10 Second difference of potential (V)
Vertical position of receivers (m)
Vertical position of receivers (m)
0.75
1
Case I Case II Case III Case IV
Vertical position of receivers (m)
Vertical position of receivers (m)
1
5 Ω· m
0.25 0 −0.25 −0.5
10000 Ω· m 0.01 Ω· m
−0.75 −1
2
10
Case V Case VI Case VII
0.75 0.5
10
5 Ω· −2 m
10
0
10 Relative error (in %)
2
10
Fig. B.6. Simulated through-casing resistivity measurements in a sixty-degree deviated well. Size of computational domain: 200 m (vertical) × 100 m (horizontal). Different curves correspond to different algorithms summarized in Table B.2. Top panel: solution in the quantity of interest as a function of logging position. Bottom panel: Relative error with respect to the reference solution corresponding to algorithm (case) VIII.
40
5 Ω· m
0.5 0.25
10000 Ω· m
0 0.01 Ω· m
−0.25 −0.5 −0.75 −1
5 Ω· m
−8
−6
10 10 Second difference of potential (V) 1 Vertical position of receivers (m)
Vertical position of receivers (m)
0.75
1
Case I Case II Case III Case IV
0.5
5 Ω· m
10000 Ω· m 0.01 Ω· m
−0.75 −1
5 Ω· m−8
−6
10 10 Second difference of potential (V)
0.25
10000 Ω· m
0 0.01 Ω· m
−0.25 −0.5 −0.75
5 Ω· m
−8
−6
1
0
−0.5
0.5
−4
10 10 Second difference of potential (V)
0.25
−0.25
Case V Case VI Case VII Case VIII
5 Ω· m
−1
−4
3D Case II Case IV Case VIII
0.75
0.75
10
Vertical position of receivers (m)
Vertical position of receivers (m)
1
5 Ω· m
0.25 0 −0.25 −0.5
10000 Ω· m 0.01 Ω· m
−0.75 −1
−4
10
Case V Case VI Case VII
0.75 0.5
10
5 Ω· −2 m
10
0
10 Relative error (in %)
2
10
Fig. B.7. Simulated through-casing resistivity measurements in a sixty-degree deviated well. Size of computational domain: 2000 m (vertical) × 1000 m (horizontal). Different curves correspond to different algorithms summarized in Table B.2. Top-left, top-right, and bottom-left panels: solution in the quantity of interest as a function of logging position. Bottom-right panel: Relative error with respect to the reference solution corresponding to algorithm (case) VIII.
41
0.75
1
200x100 2000x1000
5 Ω· m
Vertical position of receivers (m)
Vertical position of receivers (m)
1
0.5 0.25
10000 Ω· m
0 0.01 Ω· m
−0.25 −0.5 −0.75 −1
5 Ω· m
−8
−6
10 10 Second difference of potential (V)
5 Ω· m
0.5 0.25
10000 Ω· m
0 0.01 Ω· m
−0.25 −0.5
5 Ω· m
−0.75 −1
−4
10
200x100 m
0.75
5
10 15 Relative error (in %)
20
Fig. B.8. Simulated through-casing resistivity measurements in a sixty-degree deviated well. Different curves correspond to difference sizes of computational domain. Left panel: simulated measurements. Right panel: Relative error using the solution on the largest spatial domain (2000 m × 1000 m) as the reference solution.
42
0.25 m
10 Ω· m
100 m
Mandrel
1.75 m
Mandrel resistivity: 10000 Ω· m Mandrel thickness: 0.07 m
t
Rx2 t
1m
0.5 Ω· m
0.5 m
10 Ω· m
1m
20 Ω· m
100 m
Rx1
t
0.2 Ω· m
0.75 m
Tx
0.1 m
100 m
Fig. B.9. 2D cross-section of an axial-symmetric Laterolog resistivity problem in a borehole environment. Measurements are based on one current (emitter) and two voltage (collector) electrodes. The subsurface earth formation is assumed to be composed of four different horizontal layers with varying resistivities, from 0.5 Ω· m to 20 Ω· m.
43
2.5 2
2.5 Case I Case II Case III Case IV
10 Ω· m
1.5 Vertical position of receivers (m)
Vertical position of receivers (m)
1.5 1 0.5 0
Case V Case VI Case VII Case VIII
2
20 Ω· m 0.5 Ω· m
−0.5 −1
10 Ω· m 1
20 Ω· m
0.5 0
0.5 Ω· m −0.5 −1
10 Ω· m −1.5
−1.5
10 Ω· m −2 −2 10
−1
10 Second difference of potential (V)
−2 −2 10
0
10
2.5 2
10 Ω· m
Case V Case VI Case VII
2 1.5 Vertical position of receivers (m)
Vertical position of receivers (m)
0
10
2.5 Case I Case II Case III Case IV
1.5 1 0.5 0
−1
10 Second difference of potential (V)
20 Ω· m 0.5 Ω· m
−0.5 −1
10 Ω· m
1 0.5
20 Ω· m
0
0.5 Ω· m −0.5 −1
10 Ω· m
10 Ω· m
−1.5 −2 −1 10
−1.5
0
1
10 10 Relative error (in %)
−2 −1 10
2
10
0
1
10 10 Relative error (in %)
2
10
Fig. B.10. Simulated Laterolog resistivity measurements acquired in a sixty-degree deviated well. Different curves correspond to different versions of the simulation algorithm. Left panel: solution in the quantity of interest as a function of logging position. Right panel: Relative error with respect to the reference solution corresponding to algorithm (case) VIII.
44
1
1.5
0 degrees 30 degrees 45 degrees 60 degrees
5 Ω· m
Vertical position of receivers (m)
Vertical position of receivers (m)
1.5
0.5 10000 Ω· m
0 0.01 Ω· m
−0.5
1
5 Ω· m
0.5 10000 Ω· m
0 0.01 Ω· m
−0.5
5 Ω· m
−1 −8 10
0 degrees 30 degrees 45 degrees 60 degrees
5 Ω· m −7
−6
−5
10 10 10 Second difference of potential (V)
−1 −10 10
−4
10
−9
−8
−7
10 10 10 Second difference of potential (V)
−6
10
Fig. B.11. Simulated through-casing resistivity measurements. Size of computational domain: 2000 m (vertical) × 1000 m (horizontal). Different curves correspond to different dip angles: 0, 30, 45, and 60 degrees. Left panel: casing resistivity equal to 10−5 Ω· m. Right panel: casing resistivity equal to 2.3 × 10−7 Ω· m.
45
0 degrees
1
30 degrees 1.5
NO INV 10 cm INV 50 cm INV
5 Ω· m
Vertical position of receivers (m)
Vertical position of receivers (m)
1.5
0.5
0
10000 Ω· m 1 Ω· m −→ 0.01 Ω· m
−0.5 5 Ω· m
−1 −10 10
−9
−8
−7
10 10 10 Second difference of potential (V)
1
5 Ω· m
0.5
0
10000 Ω· m 1 Ω· m −→ 0.01 Ω· m
−0.5 5 Ω· m
−1 −10 10
−6
10
NO INV 10 cm INV 50 cm INV
−9
45 degrees NO INV 10 cm INV 50 cm INV
1 5 Ω· m
0.5
0
−0.5
−1 −10 10
10000 Ω· m 1 Ω· m −→ 0.01 Ω· m 5 Ω· m
−9
−8
−7
−6
10
60 degrees 1.5 Vertical position of receivers (m)
Vertical position of receivers (m)
1.5
−8
10 10 10 Second difference of potential (V)
−7
10 10 10 Second difference of potential (V)
1 5 Ω· m
0.5
0
−0.5
−1 −10 10
−6
10
NO INV 10 cm INV 50 cm INV
10000 Ω· m 1 Ω· m −→ 0.01 Ω· m 5 Ω· m
−9
−8
−7
10 10 10 Second difference of potential (V)
−6
10
Fig. B.12. Simulated through-casing resistivity measurements. Casing resistivity equal to 2.3 × 10−7 Ω· m. Size of computational domain: 2000 m (vertical) × 1000 m (horizontal). Different panels correspond to different dip angles: 0 (top-left), 30 (top-right), 45 (bottom-left), and 60 degrees (bottom-right). Different curves correspond to different radial lengths of water invasion into the conductive layer (layer 1): no invasion, 10 cm invasion, and 50 cm invasion.
46
0 degrees
1
30 degrees 1.5
NO INV 10 cm INV 50 cm INV
5 Ω· m
Vertical position of receivers (m)
Vertical position of receivers (m)
1.5
0.5
0
10 Ω· m −→ 10000 Ω· m 0.01 Ω· m
−0.5 5 Ω· m
−1 −10 10
−9
−8
−7
10 10 10 Second difference of potential (V)
1
5 Ω· m
0.5
0
10 Ω· m −→ 10000 Ω· m 0.01 Ω· m
−0.5 5 Ω· m
−1 −10 10
−6
10
NO INV 10 cm INV 50 cm INV
−9
45 degrees NO INV 10 cm INV 50 cm INV
1 5 Ω· m
0.5
0
10 Ω· m −→ 10000 Ω· m 0.01 Ω· m
−0.5
−1 −10 10
5 Ω· m
−9
−8
−7
−6
10
60 degrees 1.5 Vertical position of receivers (m)
Vertical position of receivers (m)
1.5
−8
10 10 10 Second difference of potential (V)
−7
10 10 10 Second difference of potential (V)
1 5 Ω· m
0.5
0
10 Ω· m −→ 10000 Ω· m 0.01 Ω· m
−0.5
−1 −10 10
−6
10
NO INV 10 cm INV 50 cm INV
5 Ω· m
−9
−8
−7
10 10 10 Second difference of potential (V)
−6
10
Fig. B.13. Simulated through-casing resistivity measurements. Casing resistivity equal to 2.3 × 10−7 Ω· m. Size of computational domain: 2000 m (vertical) × 1000 m (horizontal). Different panels correspond to different dip angles: 0 (top-left), 30 (top-right), 45 (bottom-left), and 60 degrees (bottom-right). Different curves correspond to different radial lengths of water invasion into the resistive layer (layer 2): no invasion, 10 cm invasion, and 50 cm invasion.
47
1
Aρ2 A A A
0.6
A
0.4 0.2 y
0.8
A
0.6
0
0.4
ρ1
0.2
ρ0
y
0.8
1
0
−0.2
−0.2
−0.4
−0.4
−0.6
−0.6
−0.8
−0.8
−1 −1
−0.5
0 x
0.5
−1 −1
1
−0.5
0 x
0.5
1
Fig. B.14. Left panel: Top view of the geometry describing the location of a resistivity tool eccentered in the borehole. The radius of the logging instrument is equal to ρ1 , and the radius of the borehole is equal to ρ2 , with the distance between the center of the logging instrument and the center of the borehole equal to ρ0 . Right panel: Iso-lines corresponding to the change of coordinates for borehole-eccentered measurements described in Formula B.1, with ρ0 = 0.3, ρ1 = 0.5, and ρ2 = 1.
48
1006
1007 1008 1009 1010 1011
1012 1013
1014 1015 1016 1017 1018 1019 1020 1021 1022 1023
1024 1025 1026 1027 1028 1029 1030
List of Tables
B.1 Dependence of each component of the integration upon different parameters. In the description below, “l” indicates the integration points, “u,v” are the trial and test functions, respectively, and subscripts “ζ1 ,ζ3 ” denote the horizontal and vertical components, respectively.
50
B.2 Definition of eight different algorithms used for simulations of borehole resistivity measurements.
51
B.3 CPU simulation time for a total of 80 different logging positions for different versions of the numerical algorithm (case number) for the through-casing resistivity problem described in Fig. B.5 in a sixty-degree deviated well, with the resistivity of casing equal to 10−5 Ω· m. Simulations were performed on a computer equipped with a 2 Ghz AMD Opteron Dual Core processor, although only one of the cores was activated during each execution (no parallelism). The last row describes the maximum amount of memory (in Mb) used by the solver of linear equations.
52
B.4 CPU simulation time for a total of 80 different vertical logging positions for different versions of the numerical algorithm (case number) for the Laterolog measurements described in Fig. B.9 in a sixty-degree deviated well. Simulations performed on a computer equipped with a 2 Ghz AMD Opteron Dual Core processor, although only one of the cores was activated during each execution (no parallelism).
53
49
Table B.1 Dependence of each component of the integration upon different parameters. In the description below, “l” indicates the integration points, “u,v” are the trial and test functions, respectively, and subscripts “ζ1 ,ζ3 ” denote the horizontal and vertical components, respectively. l ζ1
l ζ3
σ
X
X
∇ ζ1 u ζ1
X
∇ ζ3 u ζ3 ∇ζ1 vζ1
u ζ1
u ζ3
vζ3
p2
X
p2
X
X
p2
X X
X
σ∇ζ3 uζ3
X
X
X
∇ζ3 vζ3 σ∇ζ3 uζ3
X
X
X
∇ζ3 vζ3 σ∇ζ3 uζ3 ∇ζ1 uζ1
X
O
X
∇ζ1 vζ1 ∇ζ3 vζ3 σ∇ζ3 uζ3 ∇ζ1 uζ1
X
O
X
50
Nr. of Operations p2
X
∇ζ3 vζ3
vζ1
X
p2 p3
X
p4
X
X
p4
X
X
p5
Case Number
I
II
III
IV
1 Fourier mode used for adaptivity
X
X
X
X
5 Fourier modes used for adaptivity Final hp-grid NOT p-enriched
X
Final hp-grid globally p-enriched
X X
9 Fourier modes used for the final solution
X
X
V
VI
VII
VIII
X
X
X
X
X X
X X
X
X
15 Fourier modes used for the final solution X X X Table B.2 Definition of eight different algorithms used for simulations of borehole resistivity measurements.
51
X
X
Case Number
I
II
III
IV
V
VI
VII
VIII
Total Time (minutes)
21’
40’
39’
109’
244’
290’
286’
432’
Solver Time (%)
32%
38%
38%
46%
40%
47%
47%
52%
Integration Time (%)
17%
23%
22%
28%
28%
28%
30%
31%
Memory (Mb) 139 428 360 1052 1712 1712 1712 2704 Table B.3 CPU simulation time for a total of 80 different logging positions for different versions of the numerical algorithm (case number) for the through-casing resistivity problem described in Fig. B.5 in a sixty-degree deviated well, with the resistivity of casing equal to 10−5 Ω· m. Simulations were performed on a computer equipped with a 2 Ghz AMD Opteron Dual Core processor, although only one of the cores was activated during each execution (no parallelism). The last row describes the maximum amount of memory (in Mb) used by the solver of linear equations.
52
Case Number
I
II
III
IV
V
VI
VII
VIII
Total Time (minutes) 11’ 25’ 18’ 83’ 126’ 153’ 158’ 279’ Table B.4 CPU simulation time for a total of 80 different vertical logging positions for different versions of the numerical algorithm (case number) for the Laterolog measurements described in Fig. B.9 in a sixty-degree deviated well. Simulations performed on a computer equipped with a 2 Ghz AMD Opteron Dual Core processor, although only one of the cores was activated during each execution (no parallelism).
53
![[PDF Download] Fourier Series Full pages](https://p.pdfkul.com/img/300x300/pdf-download-fourier-series-full-pages_5a1cd1421723dd7c27483bef.jpg)
