1. Field of the Invention
The invention disclosed herein relates to analysis of resistivity data collected in geophysical exploration and, in particular, to techniques for improved speed and precision of numerical simulations.
2. Description of the Related Art
Exploration for oil, gas and other hydrocarbons uses a variety of technologies embodied in various tooling and instrumentation. One important technology makes use of resistivity measurements. In resistivity measurements, an electric current is directed into formations surrounding a wellbore and various measurements are performed. Data collected from measurements is then interpreted in a variety of ways. Techniques include use of data inversion and simulation.
Data inversion is used to determine a property of an earth formation from the data collected by the resistivity measurements. To generate a data inversion algorithm, a simulation of a resistivity tool is used to estimate the data that would be obtained from the tool using known properties of the formation. Both data inversion and simulation require an accurate and precise method of simulating measurements performed by the resistivity tool.
It is known that in existing methods, (such as for two-dimensional electromagnetic simulation, as may be used in studies and data inversion), simulation precision deteriorates when a grid cell is one of heterogeneous or filled with a highly conductive material (such as metal), so that its size is much greater than the skin depth.
Thus, what are needed are techniques for overcoming drawbacks of the classical “parallel-serial” way of material averaging (inside a heterogeneous discretization cell), both by theoretical analysis and numerical experiments. Preferably, the technique provides precise numerical results over existing solutions. Among other things, the techniques should enable accurate simulation of problems with thin skin layers on coarse grids.
Disclosed is a method for simulating a conductivity measured by a resistivity tool disposed in a borehole penetrating an earth formation, the method including: constructing a main grid extending from the borehole into the formation; constructing an auxiliary grid that is local to the main grid and fitted to a boundary between a first region having a first conductivity and a second region having a second conductivity; discretizing a problem operator on the auxiliary grid wherein the problem operator applies Maxwell's Equations to a volume represented by the auxiliary local grid and the Maxwell's Equations include a term representing the conductivity of a material in the volume; calculating one-dimensional functions for points on the auxiliary grid; using the one-dimensional functions for construction of basic functions on the main grid; discretizing by the basic functions for points on the main grid; solving the resulting set of equations to obtain the simulated measured conductivity for the resistivity tool; and providing the simulated measured conductivity as output.
Also disclosed is a method for estimating a property of an formation penetrated by a borehole, the method including: conveying a logging tool in the borehole; performing resistivity measurements with the logging tool; simulating resistivity measurements performed by the logging tool by implementing the following steps: constructing a main grid extending from the borehole into the formation; constructing an auxiliary grid that is local to the main grid and fitted to a boundary between a first region having a first conductivity and a second region having a second conductivity; discretizing a problem operator denoted as Ai+1/2h on the auxiliary grid wherein the problem operator applies Maxwell's Equations to a volume represented by the auxiliary local grid and the Maxwell's Equations include a term representing the conductivity of a material in the volume; calculating one-dimensional functions Vi(r) and Vi+1(r) for points r on the auxiliary grid such that
A
i+1/2
h
V
i(r)=0,Vi(ri)=1,Vi(ri+1)=0;
A
i+1/2
h
V
i+1(r)=0,Vi+1(ri)=0,Vi+1(ri+1)=1;
discretizing by basic functions Fi,jφ=Vi(r)Lj(z) for points r and z on the main grid, where Lj(z) is a one dimensional continuous and piecewise linear function such that Lj(z)=1 at a grid node zj and Lj(z)=0 at all other grid nodes; solving the resulting set of equations to obtain the simulated measured conductivity for the resistivity tool; comparing the resistivity measurements to the simulated measurements; and
estimating the property from the comparing.
Further disclosed is a machine-readable medium having machine-executable instructions for simulating a conductivity measured by a resistivity tool disposed in a borehole penetrating an earth formation, the instructions including the following steps: constructing a main grid extending from the borehole into the formation; constructing an auxiliary grid that is local to the main grid and fitted to a boundary between a first region having a first conductivity and a second region having a second conductivity; discretizing a problem operator denoted as Ai+1/2h on the auxiliary grid wherein the problem operator applies Maxwell's Equations to a volume represented by the auxiliary local grid and the Maxwell's Equations include a term representing the conductivity of a material in the volume; calculating one-dimensional functions Vi(r) and Vi+1(r) for points r on the auxiliary grid such that
A
i+1/2
h
V
i(r)=0,Vi(ri)=1,Vi(ri+1)=0;
A
i+1/2
h
V
i+1(r)=0,Vi+1(ri)=0,Vi+1(ri+1)=1;
discretizing by basic functions Fi,jφ=Vi(r)Lj(z) for points r and z on the main grid, where Lj(z) is a one dimensional continuous and piecewise linear function such that Lj(z)=1 at a grid node zj and Lj(z)=0 at all other grid nodes; and solving the resulting set of equations to obtain the simulated measured conductivity for the resistivity tool.
The subject matter, which is regarded as the invention, is particularly pointed out and distinctly claimed in the claims at the conclusion of the specification. The foregoing and other features and advantages of the invention are apparent from the following detailed description taken in conjunction with the accompanying drawings in which:
One instrument for making resistivity measurements is available from Baker Hughes, Incorporated of Houston, Tex. The instrument, referred to as an “Earth Imager,” has provided for a variety of resistivity images.
Reference may be had to
Turning now to
As is known in the art, at least one of the imaging instrument 20 and the service truck 15 include electronics and support equipment as is needed for operating the imaging instrument 20. Included with the electronics and support equipment may be a power supply for providing power to the imaging instrument 20, processing capabilities, data storage, memory and other such components. The power provided to the imaging instrument 20 may be delivered over a broad range of frequencies, f, and currents, I.
In some embodiments, the power supply for the instrument 20 provides alternating current (AC) that is in a relatively high frequency, f, range (for example, of about 1 MHz to about 10 MHz). However, the instrument 20 may be operated at frequencies above or below this range. For example, the instrument 20 may be operated in ranges from about 100 kHz to 100 MHz. Alternatively, the instrument 20 may be used with direct current (DC) if desired.
As known to those skilled in the art of resistivity imaging, there are two major classes of instruments. In a first class, “induction” instruments generally provide primary excitation by a magnetic field generated by an induction transmitter. The induction transmitter transmits electromagnetic energy into a formation. The electromagnetic energy induces circulating currents in the formation. The magnitudes of the currents are related to conductivities of the formation. The currents in turn induce secondary electromagnetic energy in a receiver at the induction instrument. The magnitude of the secondary electromagnetic energy is related to the conductivities. Thus, by measuring the magnitude of the induced secondary electromagnetic energy, the conductivities can be determined. The second class, “galvanic” instruments use a set of electrodes or electrical dipoles to generate an electrical field due to a potential between them and thus drive currents into the materials surrounding the instrument. The amount of driven current is measured and related to conductivities of the materials.
Signal measurement may include known techniques for analog signal processing and digital signal processing as appropriate. As these and other aspects are known in the art, such aspects are generally not discussed in greater detail herein. It is noted here that the term “conductivity” may be used in place of resistivity recognizing that conductivity is the inverse of resistivity.
With regard to the prior art, numerical solution of a differential problem by finite-difference or finite-element methods requires discretization, (i.e. division of a computational domain into elementary cells (also referred to as “grid generation”)) with succeeding construction of a corresponding approximate discrete problem. In the majority of practical cases, a discretized problem has discontinuous coefficients. That is, solution to the problem is disrupted on interfaces between subdomains. It is well known that the best precision is achieved when all these interfaces are fitted by the grid. But it is often impossible for different reasons. For example, the problem may be technically undoable, or lead to an excess quantity of unknown variables.
Thus, users of prior art techniques are provided with a problem of correcting approximations in a heterogeneous grid cell, for example, one that is intersected by at least one interface. Prior art attempts to address this problem make use of “material averaging” in such cells. That is, the prior art calls for replacing true varying coefficients by a constant effective tensor derived with a “parallel-serial” rule. It is demonstrated here, by theoretical analysis and numerical experiments, that this approach is inadequate, at least in some cases.
Accordingly, a new approach, which is applicable to two-dimensional (2D) simulation with “pancake” earth models is disclosed herein. The teachings provide special basic functions of the finite-element method. The functions are provided in such a way that they are conformal with varying problem coefficients in a heterogeneous cell. To do this, an auxiliary local discretization of the second level is used, with a grid fitted to all inner interfaces.
Advantageously, the technique also provides users with adequate resolution of thin skin layers on coarse grids. To achieve this, the auxiliary grid just needs to be chosen fine enough, (e.g., with a stepsize governed by the skin depth).
In order to provide some context, first consider that simulation in the frequency domain of measurements from induction resistivity tools and galvanic resistivity tools, calls for solving (three-dimensional, or 3D) Maxwell's equations, in the “electric” formulation:
∇×((iωμ)−1∇×E)−{tilde over (σ)}E=Jout
or (equivalently) in the “magnetic” one:
∇×({tilde over (σ)}−1∇×H)−iωμE=∇×({tilde over (σ)}−1Jout).
Here E and H represent an electric field and a magnetic field (complex-valued) respectively, ω represents an angular frequency, {tilde over (σ)}=σ−iωε represents a generalized conductivity tensor, μ represents a magnetic permeability tensor, and Jout represents an exciting electrical current.
For the sake of the finite element discretization, these differential problems are transformed into integral (Galerkin) forms by multiplying both sides on an arbitrary (but smooth enough) vector function F, integrating over the whole space, and applying the Green's formula:
for any F
or, correspondingly,
A finite element discretization technique is now presented. Let us search an approximate solution Eh of (1) as:
where Fih are some basic functions.
Substituting (5) into (3) and requiring the equality to be satisfied for the basic functions Fih, i=1, . . . , N, only, we get a system of linear equations for the coefficients ui:
Integral formulation (2) is discretized analogously.
Now, let us assume that a tool is placed into an axially symmetric earth model (see FIG. 3)—i.e. none of the electric parameters depends on the azimuth φ, and a source is axially symmetric as well. In this case the azimuthal dependence of a solution is known apriori:
E(r,φ,z)=E(r,z)exp(ikφ), H(r,φ,z)=H(r,z) exp(ikφ), so the problem becomes two-dimensional. In the majority of cases, we have k=0 (laterolog galvanic tools; induction tools with a longitudinal transmitter), but sometimes k=1 (induction tools with a transversal transmitter). Hereafter the case k=0 is considered only.
If a tool is induction excited by a longitudinal magnetic dipole then the electric formulation is used, where the azimuthal component Eφ is non-zero only. If a tool is galvanic excited by a longitudinal electric dipole then the magnetic formulation is used, where also the azimuthal component Hφ is non-zero only.
In the 2D case with k=0, the 3D basic functions (3) are reduced to the following construction. We build a grid:
(ri,zj),i=0, 1, . . . , nr, j=1, . . . , nz;r0=0. (5)
Then we define auxiliary one-dimensional function Li(x) where x is any of the gridded space directions from (5). Function Li(x) equals 1 at a grid node xi, 0 at all other nodes, and it is continuous and piecewise linear.
The basic functions for the azimuthal field component are defined as follows (using the double numeration as in (5)):
Let us further assume that the earth model consists of NL horizontal layers, and inside each layer, electric parameters do not depend on z (as in
˜NLnr3 ln(nzNL−1). (7)
A classical material averaging technique is now presented. The best precision of the simulation for the considered 2D model is achieved when all horizontal and vertical boundaries are fitted by the grid (5) (i.e. their positions are grid nodes). It is affordable for the horizontal ones (so it is always done in FEMAX), but, in many cases, it is not so for vertical boundaries. Namely, when there are a lot of layers each having its own vertical interfaces, an attempt to fit all of them leads to great increase of nr and, thus, to unacceptable arithmetic expenses (see the cubic dependence in estimate (7)).
Thus, some of grid cells in some layers can be intersected by interfaces between radial zones, and the problem of correct approximation in them appears. Consider first the classical “parallel-serial” material averaging. Let an inner interface {circumflex over (r)} be inside a radial step (ri,ri+1, ri<{circumflex over (r)}<ri+1, so that the generalized conductivity is:
{tilde over (σ)}(r)={tilde over (σ)}1 for ri<r<{circumflex over (r)},{tilde over (σ)}(r)={tilde over (σ)}2 for {circumflex over (r)}<r<ri+1
(see
where each cylindrical component is derived from the assumption that the electric current is either fully radial, or fully azimuthal, or fully vertical. So far the current has to be divergent-free:
its components are taken as:
J
r=const·r−1;Jφ=const;Jz=const;
which leads to the following formulas:
This averaging is easily generalized for a case of several interfaces. Formulas (8) are not defined at the first radial step (where ri=r0=0), so no averaging is done there.
The classical method is now analyzed. To analyze this approach, let us compare two cases, see
1) some galvanic tool is simulated;
2) the first zone is borehole filled with conductive mud;
3) the second zone is a resistive formation.
In this case, the current may gradually change its direction from almost radial to almost vertical inside this cell (when approaching the borehole wall). Obviously, the zone swapping depicted in
This theoretical analysis is supported by the numerical experiment presented in
borehole 8.5″, mud resistivity 0.1 Ohmm;
formation with resistivity of layers 0.1; 1; 10; 102; 103; 104 Ohmm;
each layer is of thickness 50 m except for the first and the last ones, which are infinite.
The geometrically stretching radial grid had 22 steps; the borehole radius 0.10795 m was not a node of the original grid (5), lying in between r3=0.07 m and r4=0.15 m. Two logs are compared: for this non-fitted grid with averaging (8) (the dashed line in
It is seen that these logs are truly close to each other for the formation-mud contrast not exceeding 100, but for its greater values the classical averaging produces huge error.
Now, a new approach for the considered problem is presented. Unlike in the classical one, no equivalent tensor is built; instead, the radial dependence of the basic functions is changed, from (6) to the form conformal with varying coefficients in a heterogeneous cell. The algorithm for this is the following one: for each horizontal layer and for each radial step (ri,ri+1):
r
i
=r
i,0
<r
i,1
< . . . <r
i,n−1
<r
i,n
=r
i+1; (9)
A
i+1/2
h
V
i(r)=0,Vi(ri)=1,Vi(ri+1)=0;
A
i+1/2
h
V
i+1(r)=0,Vi+1(ri)=0,Vi+1(ri+1)=1;
F
i,j
φ
=V
i(r)Lj(z) (10)
where functions Vi(r) have been constructed in stage 3.
Some remarks are presented concerning the new approach:
The calculation of Vi(r) is computationally cheap being just solving a set of 1D tridiagonal problems.
Functions Vi(r) are in general complex-valued (for nonzero frequencies).
The new functions (solid lines) are illustrated by
The result of application of the new approach to the benchmark from the previous Section is presented in the same
The new approach can simultaneously solve another important problem—correct resolution of the skin effect. It may so happen that the skin depth
in some problem zones (usually in metals) is significantly smaller than stepsizes of the main grid (5), but not so small that this zone could be considered as a perfect conductor. In this case, the conventional discretization (4)-(6) becomes imprecise. The situation can be remedied by corresponding grid refinement, to make the stepsizes to be comparable with the skin depth; but it can lead to significant increase of the computational expenses (see again estimate (7)).
Let us now apply the new approach in such a way that the auxiliary radial grid (9) satisfies the additional condition:
r
i,j
−r
i,j−1≈δ/2.
As to the vertical grid—there is no problem to make it very fine due to the logarithmical dependence of the expenses on nz in (7).
The approach was tested on the following benchmark time-domain problem (see
The problem is solved by the Fourier transformation from the frequency domain. Adequate resolution of the skin layer in the pipe is very important because it is the main source of the signal.
The simulation results are presented in
In all cases, the pipe boundary was fitted by the grid.
We see that, while the first curve is heavily erroneous, the third one practically coincides with the precise solution #2. Moreover, the calculation of the third curve was about 200 times faster than of the second one.
It has been demonstrated that the classical “parallel-serial” way of material averaging is not adequate in some cases. The method newly proposed for 2D simulation provides much better precision. Also, it enables accurate resolution of thin skin layers on coarse grids.
In support of the teachings herein, various analysis components may be used to implement the method of the new approach, including a digital and/or an analog system. The system may have components such as a processor, storage media, memory, input, output, communications link (wired, wireless, pulsed mud, optical or other), user interfaces, software programs, signal processors (digital or analog) and other such components (such as resistors, capacitors, inductors and others) to provide for operation and analyses of the apparatus and methods disclosed herein in any of several manners well-appreciated in the art. It is considered that these teachings may be, but need not be, implemented in conjunction with a set of computer executable instructions stored on a computer readable medium, including memory (ROMs, RAMs), optical (CD-ROMs), or magnetic (disks, hard drives), or any other type that when executed causes a computer to implement the method of the present invention. These instructions may provide for equipment operation, control, data collection and analysis and other functions deemed relevant by a system designer, owner, user or other such personnel, in addition to the functions described in this disclosure.
Further, various other components may be included and called upon for providing for aspects of the teachings herein. For example, a power supply (e.g., at least one of a generator, a remote supply and a battery), vacuum supply, pressure supply, cooling component, heating component, motive force (such as a translational force, propulsional force or a rotational force), magnet, electromagnet, sensor, electrode, transmitter, receiver, transceiver, antenna, controller, optical unit, electrical unit or electromechanical unit may be included in support of the various aspects discussed herein or in support of other functions beyond this disclosure.
Elements of the embodiments have been introduced with either the articles “a” or “an.” The articles are intended to mean that there are one or more of the elements. The terms “including” and “having” are intended to be inclusive such that there may be additional elements other than the elements listed. The conjunction “or” when used with a list of at least two terms is intended to mean any term or combination of terms. The terms “first” and “second” are used to distinguish elements and are not used to denote a particular order.
It will be recognized that the various components or technologies may provide certain necessary or beneficial functionality or features. Accordingly, these functions and features as may be needed in support of the appended claims and variations thereof, are recognized as being inherently included as a part of the teachings herein and a part of the invention disclosed.
While the invention has been described with reference to exemplary embodiments, it will be understood that various changes may be made and equivalents may be substituted for elements thereof without departing from the scope of the invention. In addition, many modifications will be appreciated to adapt a particular instrument, situation or material to the teachings of the invention without departing from the essential scope thereof. Therefore, it is intended that the invention not be limited to the particular embodiment disclosed as the best mode contemplated for carrying out this invention, but that the invention will include all embodiments falling within the scope of the appended claims.
This application is filed under 37 CFR §1.53(b) and claims priority under 35 U.S.C. §119(e) to U.S. Provisional Patent Application No. 61/043,057, entitled “Method of Finite Element Discretization in Heterogeneous and Highly Conductive Grid Cells,” filed Apr. 7, 2008, which is incorporated herein by reference in its entirety.
Number | Date | Country | |
---|---|---|---|
61043057 | Apr 2008 | US |