Electromagnetic Well Logging: Models for MWD / LWD Interpretation and Tool Design

By Wilson C Chin

Almost all publications on borehole electromagnetics deal with idealizations that are not acceptable physically, and unfortunately, even these models are company proprietary.  On the other hand, “exact models” are only available through detailed finite element or finite difference analysis, and more often than not, simply describe case studies for special applications.  In either case, the models are not available for general use and the value of the publications is questionable.

This new approach provides a rigorous, fully three-dimensional solution to the general problem, developed over almost two decades by a researcher familiar with practical applications and mathematical modeling.  Completely validated against exact solutions and physics-based checks through over a hundred documented examples, the self-contained model (with special built-in matrix solvers and iteration algorithms) with a “plain English graphical user interface” has been optimized to run extremely fast – seconds per run as opposed to minutes and hours – and then automatically presents all electric and magnetic field results through integrated three-dimensional color graphics.

In addition to state-of-the-art algorithms, basic “utility programs” are also developed, such as simple dipole methods, Biot-Savart large diameter models, nonlinear phase and amplitude interpolation algorithms, and so on.  Incredibly useful to oilfield practitioners, this volume is a must-have for serious professionals in the field, and all the algorithms have undergone a  laborious validation process with real use in the field.

• Language: English
• Publisher: Wiley
• Release date: Mar 19, 2014
• ISBN: 9781118835203

Wilson C Chin

Wilson C. Chin, PhD MIT, MSc Caltech, fluid mechanics, physics, applied math and numerical methods, has published twenty-five research books with John Wiley & Sons and Elsevier; more than 100 papers and 50 patents; and won 5 awards with the US Dept of Energy. He founded Stratamagnetic Software, LLC in 1997, an international company engaged in multiple scientific disciplines.

Book preview

Chapter 1

Motivating Ideas – General Formulation and Results

1.1 Overview

The general, three-dimensional, electromagnetic problem in layered anisotropic media with dip is solved using a full finite difference, frequency domain solution to Maxwell’s equations that does not bear the inherent limitations behind Born, geometric factor, hybrid and linearized integral equation approaches. Several important physical capabilities are introduced. First, transmitter coils, no longer represented by point dipoles, are modeled using eight azimuthally equidistant nodes where complex currents are prescribed. The coil may reside across multiple beds, a feature useful in modeling responses from thinly laminated zones; the transmitter operates in wireline coil alone or Measurement-While-Drilling (MWD) steel collar modes, with or without conductive mud or anisotropic invasion, and with or without borehole eccentricity. Because coil size and near-field details are explicitly considered, accurate simulation of charge radiation from bed interfaces (responsible for polarization horns) and Nuclear Magnetic Resonance (NMR) sensitive volume size and orientation in layered media are both assured.

Second, dipping interfaces are importantly oriented along coordinate planes, eliminating well known numerical noise effects associated with staircase grids. Transmitter and layer-conforming variable mesh systems, which expand in the farfield to reduce computational overhead, are automatically generated by the simulator. Third, costly performance penalties incurred by anisotropic staggered grid formulations are avoided in the vector and scalar potential method, where all complex Helmholtz equations are solved by modern matrix inversion algorithms that intelligently seek high gradient fields, relaxing and suppressing their numerical residuals. Fourth, rapid computing speeds, e.g., seconds to a minute on typical personal computers, make the approach invaluable for array deconvolution, NMR applications, and rigsite log and geosteering analysis. The availability of a single, self-consistent, open-source model eliminates the uncertainties associated with different proprietary formulations solved by different methodologies at different organizations.

Benchmark studies show excellent agreement with analytical dipole solutions in uniform and layered media and with classical Biot-Savart responses for finite loop coils. Suites of results are described, for responses in complicated media, with and without steel mandrels, invasion and borehole eccentricity, for a range of dip angles. Depth-of-penetration simulations, for electric and magnetic fields, are offered, with a view towards integrated resistivity and NMR formation evaluation. The new algorithm, which is extremely stable, fast and robust, is highly automated and does not require user mathematical expertise or intervention. It is hosted by user-friendly Windows interfaces that support approximately thirty complete simulations every hour. Fully integrated three-dimensional, color graphics algorithms display electromagnetic field solutions on convergence. Receiver voltage responses are given along tool axes, together with circumferential contributions in separate plots; detailed tabulated field results are reported in both geologically-focused rectangular and tool-oriented cylindrical coordinates. Features useful to modern logging instrument design and interpretation are available. For example, users may reconfigure transmitter coils to noncircular oblique geometries on the fly (results for elliptical cross-sections and linear geometries used in existing resistivity designs are given later). In addition, users may dynamically rewire nodal outputs in order to experiment with novel transmitter, receiver and formation evaluation concepts or to interrogate problem geologies for additional formation properties.

1.2 Introduction

The interpretation of borehole resistivity logs in layered media with dip is complicated by anisotropy, low-to-high mandrel conductivities, nonzero transmitter coil diameter, borehole eccentricity, multiple wave scattering, and polarized interfaces and charge radiation, interacting effects which cannot be studied using simplifying dipole, geometric factor, hybrid or linearized integral equation models. These approaches restrict the physics for mathematical expediency. As such, they address only specific and narrow aspects of the complete problem, e.g., purely planar layering, axisymmetric analysis, dipping bed effects modeled by vertical and horizontal dipoles, and so on.

A comprehensive model encompassing all of the above effects has been elusive. For example, real formations are not isotropic, but an anisotropic formulation covering the complete frequency spectrum plus real layering effects is not available. Moran and Gianzero (1979), for instance, deal with the induction limit only, and do not address dipping beds and the problems associated with interfacial surface charge. Howard and Chew (1992) tackle these issues, but the numerically intensive isotropic model invokes geometric factor and Born-type assumptions.

Most computational algorithms do not converge for wide ranges of frequency or resistivity and anisotropy contrast. While models are available for induction sondes with non-conducting mandrels, specialized codes are required to handle the high collar conductivities typical in MWD applications. It is usually not possible to simulate induction and MWD runs in the same formation with the same model, thus complicating interpretation and tool design.

And dipole models, often used in induction logging, are inappropriate to MWD because high conductivity collars and large coils preclude simple description. Furthermore, such tools typically log horizontal wells, where cross-sections often reside across multiple thin beds with thicknesses comparable to coil diameters. Point-wise models, moreover, cannot simulate near-field sensitive magnetic volumes accurately, a requirement that bears increasing importance with the acceptance of NMR logging and the need for improved tool design. As noted, vertical and horizontal dipole superpositions are used to model dipping beds, but this breaks down for the large coils found when tool diameters and layer thicknesses are comparable. Improved source models are necessary to good interpretation and, of course, to next-generation array and azimuthal resistivity tool design and pinpoint rf NMR excitation.

Additional difficulties abound. For example, the continuity law σv,1E1Z = σv,2 E2Z for steady-state vertical current, is often applied incorrectly. Only in Howard and Chew (1992) is the "continuity of complex current," as derived in the classic electrodynamics book of Stratton (1941), properly invoked in transient applications. The effects of multiple wave scattering, for instance, are implicitly ignored in Born and linearized integral equation approaches, thus restricting their usage to small conductivity contrasts and dip angles. Recent finite difference methods, despite the apparent generality, are likewise prone to uncertainty. Druskin et al (1999), a case in point, devise a sequential approach that solves the static problem to leading order; it therefore represents a small frequency perturbation expansion that is not necessarily convergent. Thus, speed is achieved at the expense of accuracy. Moreover, this requires that σμωR² << 1 dimensionlessly, where R is the transmitter radius, but the implications of higher conductivities, large coil radius and length scales associated with layered media are not discussed.

Finally, practical problems associated with large equation systems persist. In the general case, eight coupled unknowns describe each spatial point, and their efficient inversion is needed in formation evaluation and geosteering applications. Issues related to fast problem setup, user-friendliness, robust solvers, low hardware cost, software licensing and efficient color graphics displays, also arise in any attempt to address the overall problem. Fortunately, all of the theoretical and practical problems discussed have been solved, and we are pleased to report the elements of a new computational approach and its implications. Our Stratamagnetic Software emXplorer™ modeling system is available for complimentary download from our cloud servers and further information may be obtained directly from the publisher or author.

1.3 Physical Model and Numerical Formulation

Our Stratamagnetic Software emXplorer™ modeling system will be explained in its entirety, but in order to facilitate its description, physical ideas are developed first, and later reinforced by mathematical analysis. This section is followed by a comprehensive one describing validation procedures and practical real-world examples. The application of the model to new ideas in hardware design and formation evaluation is undertaken, and details related to overall software specification and implementation are given.

1.3.1 Design philosophy.

The fully three-dimensional model is designed to honor most of the geometrical and geological details of the logging tool and formation. As such, it may run slower than zero or one-dimensional models, but faster than many three-dimensional models, requiring at most one minute or so per simulation on typical computers. Its objective is easily stated: allow users to test assumptions about the formation, perhaps simplifying ones that justify the use of more elementary electromagnetic models. Much of this book is devoted to results that support our physical validations which are of interest to instrument designers. However, the equally important needs of petrophysical analysts are also addressed; for example, receiver responses both axially and azimuthally can be plotted, and log generation capabilities are available.

1.3.2 New discretization approach.

A comprehensive model in which all of the salient physical features are retained was developed within a finite difference, frequency domain framework. Here, the numerical noise associated with the grid sizes and aspect ratios of conventional staircase approximations to dipping interfaces does not appear, because coordinate surfaces aligned with local bedding planes define special oriented grids (these are not to be confused with the staggered grids discussed later).

Figure 1.1. Staircase versus boundary conforming grids for layered media.

Staircase grids, originally used in computational aerodynamics in the 1970s, have fallen in favor relative to boundary conforming meshes, which resolve derivatives parallel and perpendicular to boundaries more precisely. In other words, classical boundary conditions cannot be accurately implemented despite grid refinement although, on a case-by-case basis, acceptable agreement with analytical solutions can be achieved by trial and error. At the same time, crude point source models of airfoils were replaced by distributed line singularities; these were more stable than point singularities, since they are less infinite. Taken together, these two methods streamlined analysis and reduced computing, vastly improving performance and predictive capabilities. This grid plus source point technology is adapted to the present modeling work.

In our approach, local vertical, z planes are aligned with bed interfaces. Consistent with aerodynamic practice, the transmitter coil is modeled by eight azimuthally equidistant points where current is specified. Although we assume constant frequency and real currents, the underlying iterative model actually allows worst case complex excitation, so that the numerically stable model applies to Fourier components of pulsed transient systems as well.

To enhance resolution, the coil is always discretized taking six constant meshes across in the near-field, as in Figure 1.2; this mesh is geometrically expanded in the farfield to conserve memory and reduce unnecessary computation as in Figure 1.3. Distributed sources are less singular than point sources; they promote stability and increase convergence speed. The six across (or, eight around) source model resolves transmitter coil geometry well and provides the framework for approximate steel drill collar modeling. As shown in Figures 1.4 and 1.5, the wireline coil alone or MWD drill collar nature of any tool is modeled by twenty-one internal points, which may or may not reside across multiple formation layers. This provides the resolution needed for resistivity modeling in thinly laminated zones and NMR magnetic volume simulation in layered media. These models apply to tools of all diameters.

Figure 1.2. Nondipolar source model in layer-oriented coordinates.

Figure 1.3. Variable expanding meshes conserve memory.

Figure 1.4. Detailed transmitter coil models.

Figure 1.5. Eccentered coil with invasion.

1.3.3 Analytical formulation.

The theory behind Maxwell’s equations appears in Stratton (1941). His work can be specialized to transversely isotropic media, leading to the anisotropic model of Moran and Gianzero (1979). Different authors solve different forms of these equations, and in order to understand the differences, we derive all underlying relationships from first principles. Because much of the derivation is generally available, we will sketch the overall approach, and reserve detailed commentary only to original research results.

We begin with Maxwell’s equations for E and B in the usual form, i.e., ∇ × E + ∂B/∂t = 0, ∇ × H − ∂D/∂t = J, ∇ • B = 0, and ∇ • D = Θ. If we resolve J into a source Js plus a conduction current E, where is a diagonal conductivity tensor [σh, σh, σv] and h and v are horizontal and vertical directions in Figure 1.2, the second equation becomes ∇ × H − ∂D/∂t- E = Js. We assume D = ε E and B = μ H where ε and μ are (for simplicity) isotropic inductive capacities. When σh and σv are equal, we recover the model of Druskin et al (1999). In that work, the authors deal with E and B directly. Thus, the Yee (1966) staggered grid algorithm must be used, in which the dependent variables are evaluated alternatively at grid centers and edges, so that fictitious currents do not arise.

Davydycheva and Druskin (1995) extend this formalism to anisotropic media, but the computational overhead is substantial, with workstation solutions requiring hours per run. Together, these staggered plus staircase mesh systems pose formidable obstacles to obtaining useful solutions. But the most severe limitation appears in the iterations. Druskin et al (1999) note that a static solution is solved first, that is in turn improved by frequency-dependent corrections. As such, the method is, whatever the formalism, a de facto expansion in small frequency that is not necessarily convergent. Thus, speed is achieved at the expense of accuracy. Consider, for example, homogeneous media, where the only length scale is the coil radius R. Because σμωR² is the sole dimensionless quantity physically possible, the scheme implicitly requires σμωR² << 1. Convergence is not addressed, and neither are the implications of high conductivity, larger coil radius and additional length scales arising in layered media. The only apparent benefits are solutions for multiple frequencies obtained with minimal effort, but this is possible with any perturbation scheme taken in powers of σμωR².

1.3.4 An alternative approach.

In problems excited by external currents, Feynman et al (1964) note that it is more natural to solve the potential form of Maxwell’s equations. This is especially true in practice, because the formulation involves classical differential operators (Courant and Hilbert, 1989) whose solution algorithms are widely available. These do not require staggered meshes, nor do they produce spurious solutions associated with direct field approaches. But complications related to boundary conditions do arise, which we have identified and addressed.

A completely equivalent alternative to "E and B is A and V," where the latter represent, respectively, well known vector and scalar potentials. As before, the basic derivation is standard, and can be extended to anisotropic conductivity by introducing straightforward changes. The idea is to represent, without loss of generality, the magnetic and electric fields in the form B = ∇ × A and E = − ∂A/∂t − ∇V. When the Lorentz gauge ∇ • A = − με ∂V/∂t − μσhV is used, the equations governing A and V are found as ∇²A − μ ∂A/∂t − με ∂²A/∂t² − μ ( − σh)∇V = − μ Js and ∇²V − μσh ∂V/∂t − με ∂²V/∂t² = − Θ/ε.

These are satisfying and self-consistent because both are diffusive at low frequencies and wave-like at high frequencies. We remark that the propagation resistivity alluded to in MWD/LWD logging actually satisfies the diffusive limit of the governing equations – true wave propagation effects are not encountered except for dielectric applications.

Real world heterogeneities assume two forms. They can be slowly varying, and if so, they manifest themselves as variable coefficients, e.g., this may be the case with formations undergoing mudcake dominated invasion. Or they can be rapidly varying, in which case the differential equations themselves cease to be valid, since locally infinite derivatives do not exist. This is true of the abrupt changes that occur from layer to layer. For such problems, external interfacial matching conditions must be invoked which are also generally formulated in Stratton (1941).

This matching takes the general form n • (B2 − B1) = 0, n • (ε2E2 − μ1E1) = θ, n × (E2 − E1) = 0, and n × (μ1B2 − μ2B1) = 0 where 1 and 2 represent any two adjacent layers, and n is the unit normal to the interface (any number of layers is allowed). Because we choose to deal with A and V, these conditions must be expressed in the new variables. When a proportionality to e+iωt is assumed, and all vector results are expanded Cartesian form, we have the sequence of coupled boundary value problems for the italicized modes Ax, Ay, Az and V in Equations 1.1 to 1.4. These form the basis for our comprehensive model, with our objective being fast, robust and stable solutions. In order to achieve these goals, their properties must be thoroughly understood.

Ax formulation –

(1.1a)

equation

(1.1b)

equation

(1.1c) equation

(1.1d)

equation

(1.1e) equation

Ay formulation –

(1.2a)

equation

(1.2b)

equation

(1.2c) equation

(1.2d)

equation

(1.2e) equation

Az formulation –

(1.3a)

equation

(1.3b)

equation

(1.3c) equation

(1.3d) equation

(1.3e) equation

V formulation –

(1.4a)

equation

(1.4b)

equation

(1.4c) equation

(1.4d)

equation

(1.4e) equation

In the low-frequency, homogeneous medium limit, our equations reduce to the anisotropic model of Moran and Gianzero (1979). Beyond this, all similarities cease: our source is nondipolar, our interfaces are strongly polarized, and our near-field may contain collars, mud, nonuniform invasion and tool eccentricity. Any practical solution must address inter-related issues. For example, the physical problem and grid must be simple to set up and its solution must be straightforward. Computational times, memory resources and user intervention, must be kept to a minimum. Solutions must be physically correct and displayed meaningfully as they become available. These objectives, our insights and experiences, are discussed next.

1.3.5 Solution philosophy.

Although the equations for A and V are well known, and matching conditions have been previously expressed in these variables, the packaging of the boundary value problems in Equations 1.1 to 1.4 is new. Each consists of a complex Helmholtz equation for a complex dependent variable, and each satisfies standard regularity conditions at infinity. But close examination shows that we have expressed interfacial conditions in terms of left side unknown variables and their normal derivatives, while transferring all coupled dependent variables and tangential derivative terms to the right. By doing so, we have recast a complicated set of equations in near-Neumann form to facilitate a near-classical solution, without neglecting the required modal coupling between different rectangular components of A.

In our iterations, all right sides are updated continuously until convergence is achieved. This normal derivative approach is natural and is successfully used in many disciplines. For example, in heat transfer, temperature is continuous at non-insulating boundaries; continuity in heat flux means specifying jumps in normal temperature derivative (Carslaw and Jaeger, 1959). Similarly, in Darcy flow, pressure is continuous at changes in formation permeability, while velocity continuity implies jumps in normal pressure derivative; more general models do exist, for instance, at shale interfaces which support pressure discontinuities (Muskat, 1937; Chin, 2002).

Of course, while these boundary condition analogies apply to heat transfer and Darcy flow, their differential equations are real as opposed to complex. But a more precise analogy is available in aerodynamics, where a scalar potential satisfying a complex Helmholtz equation describes phase lags between aileron motion and induced pressure (Landahl, 1961). In this application, ϕ is discontinuous through wakes, while its normal derivative is continuous; at lifting surfaces, ϕ is continuous, but its normal derivative is not. Thus, aerodynamic methods are applied here, provided the algorithm for ϕ is vectorized to allow multiple A and V unknowns.

This algorithm research was undertaken by the aerospace industry and military in the 1970s since the methods are used in aerodynamic control. The solvers, forerunners to domain decomposition methods, resolve full problems into smaller domains with contrasting computing requirements. They are efficiently coded, using in place, memory conserving methods, while intelligently updating only those spatial regions characterized by rapid physical gradients. These highly specialized matrix solvers have been adopted for use in the present application.

1.3.6 Governing equations.

In Equations 1.1 and 1.2 for Ax and Ay, the dependent variables are continuous at interfaces; their normal derivatives are not, by an amount that depends on Az. By contrast, in Equation 1.3 for the vertical field Az, the normal derivative is continuous while Az itself is not. As discussed in Howard and Chew (1992) and Stratton (1941), it satisfies "(σv,1 + i ωε1)E1z = (σv,2 + i ωε2)E2z based on complex current. The more familiar σv,1 E1z = σv,2 E2z" appears in the DC limit and does not apply. Correct enforcement is essential to proper solutions for interfacial charge, which are iteratively used to excite the scalar potential V in Equation 1.4. A complicating feature of Equation 1.3 is its dependence on cross-derivatives of Ax and Ay, which arise from formation anisotropy (this effect is destabilizing in many E-B approaches, but stabilizing in A-V solutions). Once A and V are both available, magnetic and electric fields are obtained from B = ∇ × A and E = − ∂A/∂t − ∇V.

1.3.7 Finite difference methodology.

We have yet to discuss how Equations 1.1 to 1.4 are to be discretized. Because accuracy is essential and memory resources are scarce, at least in personal computer applications, second-order accurate central difference approximations are always used. Importantly, all four equation systems possess almost identical differential operators, thus simplifying algorithm development, programming and code maintenance. This similarity arises because our formulation appears in rectangular, as opposed to, say, cylindrical or spherical coordinates. These are consistent with boundary conforming grid requirements as expressed in Figure 1.2.

We do not need to list all of our finite difference equations. It suffices to consider the most complicated partial differential equation, namely the anisotropic one underlying Az, since this equation contains all of the model terms. Equation 1.5 applies to the indexed point (xl,ym,zn) and is written for all points in the computational domain. Variations of this are used in different domains with different coefficients, as the medium changes from rock to collar, to mud, to invaded zone. This is also done for Ax, Ay and V. The resulting set of equations, very large in size, is solved simultaneously for the coupled unknowns. Note how the entire right-side of Equation 1.5 vanishes in the isotropic limit when σv = σh.

(1.5)

equation

While Equation 1.5 is easily derived, involving handbook formulas only, the approach used to discretize jump conditions at interfaces requires more care. In elementary texts, standard difference formulas for smooth functions are given, assuming that the function and its derivatives are continuous. As noted, this is often not the case.

For example, when A is double-valued as in Figure 1.6b, dA/dz must be evaluated using (An+1 − An-1 − Δ)/(zn+1 − zn-1). In addition, provision must be made for extra storage, and optimal re-ordering of unknowns for efficient matrix inversion is required. When slopes are discontinuous as in Figure 1.6c, the amount of the jump must appear explicitly in the formula used. For instance, in Equation 1.2d, the left side would be differenced as usual, but the right side would contain latest nonzero values obtained for Az. These changes require detailed book-keeping and considerable care.

Figure 1.6. Functions with jumps and slope discontinuities.

1.4 Validation Methodology

The dearth of nontrivial analytical solutions for three-dimensional benchmarking is well known. While agreement with well logs represents an ultimate goal, self-consistency with other solutions to Maxwell’s equations provides the best evidence that the formulation is rigorous and its results are correctly computed. Comprehensive validation requirements were developed to identify gridding issues, formulation subtleties, potential numerical problems and programming mistakes, and are summarized below. Shortage of space precludes detailed reporting, but all validation files are available upon request.

In the following validations, a 21 × 21 × 21 variable grid was used, noting that second-order accurate differencing plus high resolution near-field meshes are employed (the transmitter coil center is always centered at 11, 11, 11 to facilitate numerical checks for symmetry and antisymmetry). Nonetheless, because eight unknowns for complex A and V obtain per point, 74,088 coupled equations must be solved. In context, a 100 × 100 axisymmetric, fine mesh calculation involves only 20,000 unknowns. We also emphasize that, while we focus on constant properties in this book, the algorithms apply to nonlinear properties as well. For instance, if the function ε = ε(A) is given, the argument A is approximated by its latest numerical value in the derived iterative scheme and the algorithm is run to convergence.

Real benchmarking addresses more than agreement with known results. Much more demanding is the requirement, Can good results be obtained quickly with this mesh, the first time around and without benefit of analytical results, in a field setting? We answer, definitively, "Yes!" Fully converged results can be stably obtained on personal computers in less than one minute, making the algorithm indispensable for general application (high values of dimensionless diffusion in most field applications guarantee fast convergence). Because total time from problem setup to solution display requires approximately a minute, a suite of a half-dozen simulations, say, can be easily defined, computed and fully documented within thirty minutes of desk-time. For example, the detailed color displays shown later are automatically generated upon convergence. We now summarize the validation types we have pursued.

1.4.1 Fundamental physics.

Axisymmetric coil-alone runs in homogeneous isotropic media were performed at several dip angles to ensure that all solutions were physically identical. This check, which ensures that logic loops behave similarly in all directions, is important because the algorithm does not presuppose any symmetries, thus opening the way to more general heterogeneities.

For vertical wells, computed antisymmetries in Ax and Ay about coil center, as well as zero center values, were always achieved. Similar remarks apply to Ax and Az in horizontal well logic loops.

In layered media, a transmitter coil oriented at nonzero dip angle induces charge distributions that form along bed interfaces. Because source current approaches an interface from one direction and departs in the opposite, induced charges must appear antisymmetrically with respect to a line having zero charge. This antisymmetry is also captured well.

We do not assume that A vanishes at the coil axis, a zero that strictly applies axisymmetrically. That A = 0 is automatically captured in such limits lends credibility to problems without geological symmetries, where centerline values are typically nonzero. For general geologies, this zero appears elsewhere in the domain and appears as part of the iterative solution.

Once A and V are available, B and E are obtained using B = ∇ × A and E = − iωA − ∇V in the frequency domain. All expected symmetries, antisymmetries and zeros in each of the B and E quantities are also independently validated in a wide range of test examples.

Higher attenuation anticipated with increases in frequency and conductivity is always achieved.

Drill collar modeling is always invokes six gridblocks across zeroing of internal fields. Setting E = 0 at collar nodes is justified since the collar is grounded in mud and no attempt is made to simulate at skin depth scales.

Numerical phase and amplitude agreement with special exact Bessel function solutions for realistic coil-alone versus steel-mandrel tools derived by Shen (1991) is crucial. This analytical solution applies to homogeneous, isotropic media only.

Receiver responses show proper axial behavior in a uniform media. For horizontal wells in three-layer media, point-wise azimuthal button responses demonstrate the correct symmetries.

1.4.2 Biot-Savart finite coil validations.

The classical result for finite radius (R) coils in infinite, homogeneous, nonconductive media states that the center magnetic field strength Baxial varies inversely with R when the real current I is fixed. This behavior is captured in all low-frequency, high resistivity simulations.

Analytical results show that Baxial decays inversely with distance cubed along the coil axis. This behavior is consistently captured.

Magnitudes for coil center Baxial obtained analytically and numerically agree for large classes of variable grids with high mesh amplification rates. Truncation errors are low because dipole fields decay rapidly in space. This also supports rapid numerical convergence of the iterative scheme.

Because our grids vary with dip angle α, it is important that calculated coil center magnetic field strengths remain constant and equal to theoretical values, if α is varied for different homogeneous medium, low-frequency, high resistivity runs. A closed form expression mapping I to Js for the custom grids used assures that this is always the case. Thus, rapidly convergent, accurate solutions are always obtained for all dip angles from 0° to 90°.

1.4.3 Analytical dipole validations.

Classical solutions for point magnetic dipoles in uniform media were used to compute the ratio Bimag/Breal at several receiver locations for different ranges of frequencies and resistivities. The known decay of Breal with distance, and the value of Breal itself, provided additional check points. All validations were successfully performed.

Similar comparisons using two-layer vertical and horizontal dipole solutions, e.g., Banos (1966) or Kaufman and Keller (1983), were also successfully undertaken.

Because the above closed form solutions are useful in grid calibration, these apps have been incorporated in pull-down menus in the software.

1.4.4 Fully three-dimensional solutions.

Our coarse mesh solutions give 25 and 31 inch receiver amplitudes and phase differences that agree with dipole results to 5% over four orders of magnitude in resistivity, simultaneously with near-field, nonzero diameter, coil alone transmitter B field results to 10%, and at the same time, agree with the Coope-Shen (1984) and Shen (1991) exact solutions, the latter incorporating steel mandrel effects, to 5–10%.

The great majority of calculations here was performed using a 21 × 21 × 21 mesh system. A number were re-run with finer 31 × 31 × 31 grids, tripling the number of unknowns and quintupling the computing time. Results showed improved accuracy, but importantly, that trends can be adequately determined using coarser meshes.

1.5 Practical Applications

1.5.1 Example 1. Granularity transition to coil source.

In this suite of calculations, the implications of the source model shown in Figure 1.2 were assessed. First, we considered a six inch coil alone tool in a homogeneous 1 Ω–m, unit ε and μ medium. Figure 1.7a gives the magnitude of the complete electric field at the transmitter plane. The eight points where current is prescribed show all expected symmetries, but display a degree of granularity not characteristic of real coils. However, this disappears quickly. Within four grids, or 70% of diameter, the smoothed field in Figure 1.7b is obtained. Thus, transmitter length scale (which is comparable to layer thicknesses) is modeled accurately. This is important in thinly laminated zones, where resistivity and NMR deconvolution can be used to unsmear averaged layer results. In our model, we assume that a constant current can be imposed through the transmitter coil and implicitly that typical excitation wavelengths greatly exceed geometric dimensions.

Figure 1.7a. Granularity in transmitter coil plane.

Figure 1.7b. Smooth field achieved at one diameter.

Next, consider a three-layer, 1,000-1-1,000 Ω–m formation, where our coil alone transmitter resides across all layers. In Figure 1.8a, high electric field attenuation occurs in the middle layer, with 0 at the center due to antisymmetry. This calculation was repeated for a coil with drill collar. Figure 1.8b shows a wider dead zone internal to the steel mandrel. These results must be contrasted with Figure 1.7b for uniform media.

Figure 1.8a. Coil alone in three-layer rock.

Figure 1.8b. Coil with drill collar run.

We now reverse the resistivity sequence, and consider 1-1,000-1. In our coil alone run, Figure 1.9a accents the strength of the middle layer. Next, Figure 1.9b shows similar results assuming a steel mandrel. For completeness, dipole results are shown in Figure 1.10, obtained by assuming a 1 in coil in our model. The length scales of the background grid are different from those used above, so direct comparisons can be misleading. But it is clear that no resemblance is found between this and foregoing results.

Figure 1.9a. Coil alone simulation.

Figure 1.9b. Coil with drill collar run.

Figure 1.10. Dipole simulation.

We emphasize that, throughout the color graphics used in this book and in the software, the incremental spaces shown at the blue bottom do not represent equal distances, but instead, numerical indexes describing a physically expanding grid. Plotting to actual scale would have, of course, provided few features to display visibly.

1.5.2 Example 2. Magnetic field, coil alone.

We study the magnetic field for a coil alone tool. In uniform media, the axial symmetry in Figure 1.11 is expected. For an anisotropic medium with Rh = 1 Ω–m, and Rv taking unit values except in a middle 1,000 Ω–m layer that contains the coil, Figures 1.12 and 1.13 show distortions in B and interfacial charge for a horizontal well application (here, Rv/Rh = 1,000 >> 1).

Figure 1.11. Axial magnetic field.

Figure 1.12. B field contour plot.

Figure 1.13. Interfacial surface charge.

1.5.3 Example 3. Steel mandrel at dip.

Here a six inch diameter, 2 MHz coil with steel mandrel is oriented at 45° dip in a three-layer 1-1,000-1 Ω–m isotropic formation, with the coil wholly contained in the middle layer. For illustrative purposes, consider the Ax field sketched and computed in Figure 1.14a and 1.14b. Perfect antisymmetries are recovered, and expected zero (green) fields are obtained in the x = 0 plane. Figures 1.14c,d display alternative views, made possible by integrated graphics capabilities allowing zoom, rotate, and translate. These may be useful in evaluating sensitive volume penetration in steel collar NMR applications.

Figure 1.14a. Deviated 45° run.

Figure 1.14b. Imaginary Ax solution.

Figure 1.14c. Contour ceiling plot.

Figure 1.14d. Underside projection.

1.5.4 Example 4. Conductive mud effects in wireline and MWD logging.

In the following simulations, the same homogeneous medium is taken throughout, and contour plots for total electric field strength near the transmitter plane are given. Figures 1.15a,b show results for coil alone and coil with drill collar in the absence of borehole mud (note the hole left by the steel mandrel). Figures 1.15c,d introduce very conductive mud, which severely attenuates the produced field; here a large σμωR² ≈ O(1) renders invalid the method of Druskin et al (1999). Observe how the red glow in Figures 1.15a,b is no longer evident in Figures 1.15c,d. Finally, in Figures 1.15e,f, we eccenter both coil alone and coil with collar tools in our hole containing conductive mud. These results show that receiver responses (which measure azimuthal averages) arise mainly from one side of the formation, a dangerous situation which is suggestive of incorrect well log interpretation.

Figure 1.15a. Coil alone, no mud.

Figure 1.15b. Steel mandrel, no mud.

Figure 1.15c. Coil alone, conductive mud.

Figure 1.15d. Steel mandrel, conductive mud.

Figure 1.15e. Eccentered coil alone, conductive mud.

Figure 1.15f. Eccentered steel mandrel, conductive mud.

1.5.5 Example 5. Longitudinal magnetic fields.

In this computational suite, we consider five horizontal well results for the inhomogeneous formations shown in Figure 1.16 with isotropic layers. A six inch, 20 KHz coil alone induces strong interfacial charge polarization in all cases. But we are less interested here in resistivity logging, as we are in NMR sensitive volume imaging; as such, we study the longitudinal magnetic fields associated with our circular transmitter coils.

Figure 1.16. Five baseline logging scenarios.

The results in Figures 1.17a to 1.17f, characterizing the formation one diameter away from our coil, are easily interpreted using Figure 1.11. Red indicates the strong fields associated with the coil core. On the other hand, yellow, green, and blue represent decreases in strength as the coil wire itself is approached, while depressions in the blue fabric describe the turning of field lines as they form closed loops. In this work, transmitter and receiver coils are circular with planes perpendicular to the tool axis, however, they may in general take any closed or open form; minor source code changes to the Js vector which excites A will be required.

Figure 1.17a. Axial magnetic field, Case 1.

Figure 1.17b. Axial magnetic field, Case 2.

Figure 1.17c. Axial magnetic field, Case 3.

Figure 1.17d. Axial magnetic field, Case 4.

Figure 1.17e. Split magnetic field, Case 5.

Figure 1.17f. Split, underside, Case 5.

Noticeable differences are seen, but the most pronounced is shown for Case 5, where the transmitter coil resides across two adjacent layers. Here, the magnetic field in Figure 1.17e is split by strong heterogeneities, the way woodcutters split fire logs. Its underside, shown in Figure 1.17f, highlights this bias. This result bears strong implications in both resistivity and NMR logging. Instrument readings, and possibly, interpretation software, are likely to smear individual layer results, suggesting instead, averaged formation properties. Consequently, saturation estimates can be calculated incorrectly, and oil and gas producing zones can be bypassed unintentionally.

1.5.6 Example 6. Elliptical coils.

In the above runs, axial magnetic fields in horizontal wells with circular coils were studied. High conductivity contrasts, it was seen, can prevent logging tools from reaching equally into adjacent layers, thus prompting the question, Are there transmitter shapes that are ‘optimal’ for penetration depth?

In our model, transmitter coils are easily reconfigured to any desired geometry (source code access is required). For instance, an eight point ellipse is shown in Figure 1.18a, with a major-to-minor axis aspect ratio of two. Such cross-sections are viable candidates for modern resistivity tool and NMR design. Computed focused magnetic fields, shown in Figures 1.18b to 1.18e, are difficult to anticipated a priori.

Figure 1.18a. Eight point ellipse.

Figure 1.18b. Axial magnetic field.

Figure 1.18c. Contour plot.

Figure 1.18d. Underside plot.

Figure 1.18e. Side view.

1.5.7 Example 7. Calculating electromotive force.

We have seen how transmitters can be reconfigured, so that they might optimally excite formations consistently with the heterogeneities. For any given transmitter, differences in formation properties can also be interrogated by custom designed receivers. Special modeling options permit rewiring of nodal outputs in order to evaluate novel antenna concepts. In general, the electromotive force between a and b obtains from emf = ∫ab E • dl. For closed circuits, Stokes’ theorem yields Lenz’s law emf= − ∂/∂t ∫∫ B • dS, which leads to the usual "- iωSBaxial" used in induction logging where S is enclosed area. But a direct application of

(1.6)

equation

where a to b is open, or perhaps, one or more unconventionally directed closed loops which may or may not wrap around the collar, may uncover more petrophysical information than is possible with existing coil arrays.

For example, "V(a) − V(b)" is an indicator of interfacial charge intensity, that is, formation dip and conductivity contrast, while the integral measures the usual eddy current effects. The antenna of the future may well be an intelligent array that reconfigures itself as it logs a formation, perhaps, in response to similar scenarios stored in a dynamic database.

We are only now embarking on transmitter and receiver studies. To support such endeavors and the needs of petrophysical log analysts, the three-dimensional field results emphasized thus far are also post-processed to provide log-oriented data. For example, Figures 1.19a,b show real and imaginary (that is, in and out-of-phase) voltage responses for standard receiver coils distributed at nodal positions along the tool axis. For illustrative purposes, the calculations shown assume a uniform medium; the anticipated symmetry with respect to maxima at the transmitter plane is clearly evident. In actual multilayer calculations, this symmetry is broken, of course, and the resulting response curves may contain considerable formation information. But even so, much detail, is lost since coil readings always average individual azimuthal responses. To assist log analysts in evaluating their models, we also provide incremental circumferential responses to coil readings, as explained in the next simulation suite.

Figure 1.19a. In-phase coil voltage.

Figure 1.19b. Out-of-phase coil voltage.

1.5.8 Example 8. Detailed incremental readings.

The above results emphasize total voltage, but our receiver models actually sum incremental readings at button nodes, which make up the emf integral given earlier. For example, voltages from each of the eight coil segments in Figure 1.20a are shown in Figures 1.20b,c.

Figure 1.20a. Three-layer model.

Figure 1.20b. In-phase increments.

Figure 1.20c. Out-of-phase increments.

1.5.9 Example 9. Coil resting along bed interface.

In Figures 1.17e,f, we presented a horizontal well example where a standing transmitter coil resides across two layers. Finally, let us consider a two-layer, high resistivity contrast (1:1,000) vertical well problem where the sleeping coil lies entirely on the interface itself, as shown in Figure 1.21a.

Figure 1.21a. Vertical well problem.

This simple problem is challenging because both current excitation and interfacial matching logic must be implemented along the same z-layer plane. Figures 1.21b,c display axial receiver voltage responses on each side of the transmitter coil. The first, almost symmetric, represents the expected in-phase component. The second, strongly dependent on conductivity, shows the bias resulting from the assumed resistivity contrast; for example, one side displays the large voltage drops associated with high conductivity.

Figure 1.21b. In-phase receiver voltage.

Figure 1.21c. Out-of-phase receiver voltage.

1.6 Closing Remarks

Although we have focused on the powerful simulation capabilities of the new model, our ultimate objectives are really directed more at modern logging than they are at numerical analysis, and particularly, electromagnetic behavior. Electromagnetics provides the common bond linking resistivity logging to NMR, and ultimately, to fluid imaging. In both cases, it provides the external stimulus. And in both cases, the response is smeared across layers, insomuch as existing tools lack the ability to penetrate deeply or resolve details finely.

Simulation provides the key to deconvolving measured signals, unlocking hidden details of the formation. And good deconvolution requires a modern capability unrestricted by the limitations of old analytical approaches or cumbersome numerical procedures. The present model fulfills this need and our three-dimensional emXplorer™ simulator is available for advanced applications, e.g., triaxial array induction, full field azimuthal imaging, and so on.

As is readily apparent, the results in Figures 1.15a to 1.15f and 1.17a to 1.17f imply that once difficult-to-answer questions related to electric and magnetic field depth of investigation can now be answered routinely and quickly. Others are just as interesting. How much power does a new tool require, say, to probe ten feet into a formation? What types of transmitters and receivers are needed to extract anisotropy, true conductivity and actual dip? Are inclined coils beneficial to formation evaluation? Can NMR sensitive volumes be focused narrowly to enhance resolution? While Figures 1.18a to 1.18e emphasized noncircular transmitters, a special user interface permits custom receiver design, allowing closed and opened coils as a post-processing option.

All of the foregoing questions can

Reviews for Electromagnetic Well Logging

[Anonymous mRMtxD6KE8 · Rating: 5 out of 5 stars]

Best