SIMULATION OF RESISTIVITY LOGGING-WHILE-DRILLING (LWD) MEASUREMENTS USING A SELF-ADAPTIVE GOAL-ORIENTED HP FINITE ELEMENT METHOD. D. PARDO† ,‡ , L. DEMKOWICZ† , C. TORRES-VERD´IN‡ , AND M. PASZYNSKI† ,§ Abstract. We simulate electromagnetic (EM) measurements acquired with a Logging-WhileDrilling (LWD) instrument in a borehole environment. The measurements are used to assess electrical properties of rock formations. Logging instruments as well as rock formation properties are assumed to exhibit axial symmetry around the axis of a vertical borehole. The simulations are performed with a self-adaptive goal-oriented hp-Finite Element Method (FEM) that delivers exponential convergence rates in terms of the quantity of interest (for example, the difference in the electrical current measured at two receiver antennas) against the CPU time. Goal-oriented adaptivity allows for accurate approximations of the quantity of interest without the need of obtaining an accurate solution in the entire computational domain. In particular, goal-oriented hp-adaptivity becomes essential to simulate LWD instruments, since it reduces the computational cost by several orders of magnitude with respect to the global energy-norm based hp-adaptivity. Numerical results illustrate the efficiency and high-accuracy of the method, and provide physical interpretation of resistivity measurements obtained with LWD instruments. These results also describe the advantages of using magnetic buffers in combination with solenoidal antennas for strengthening the measured EM signal so that the ’signal-to-noise’ ratio is minimized. Key words. hp-finite elements, exponential convergence, goal-oriented adaptivity, computational electromagnetics, Maxwell’s equations, through casing resistivity tools (TCRT). AMS subject classifications. 78A25, 78A55, 78M10, 65N50
1. Introduction. A plethora of energy-norm based algorithms intended to generate optimal grids have been developed throughout the last decades (see, for example, [7, 18] and references therein) to accurately solve a large class of engineering problems. However, the energy-norm is a quantity of limited relevance for most engineering applications, especially when a particular objective is pursued, for instance, to simulate the electromagnetic response of geophysical resistivity logging instruments in a borehole environment. In these instruments, the amplitude of the measurement (for example, the electric field) is typically several orders of magnitude smaller at the receiver antennas than at the transmitter antennas. Thus, small relative errors of the solution in the energy-norm do not imply small relative errors of the solution at the receiver antennas. Indeed, it is not uncommon to construct adaptive grids delivering a relative error in the energy-norm below 1% while the solution at the receiver antennas still exhibits a relative error above 1000% (see [13]). Consequently, in order to accurately simulate LWD resistivity measurements in this paper, we develop a self-adaptive strategy to approximate a specific feature of the solution. Refinement strategies of this type are called goal-oriented adaptive algorithms [11, 17], and are based on minimizing the error of a prescribed quantity of interest mathematically expressed in terms of a linear functional (see [3, 9, 12, 11, 17, 19] for details). † Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin, Austin TX 78712 ‡ Department of Petroleum and Geosystems Engineering, The University of Texas at Austin, Austin TX 78712 § On leave from AGH University of Science and Technology, Department of Computer Methods in Metallurgy, Cracow, Poland
1
2
D. PARDO, L. DEMKOWICZ, C. TORRES-VERD´IN, M. PASZYNSKI
In this paper, we formulate, implement, and study (both theoretically and numerically) a self-adaptive hp goal-oriented algorithm intended to solve electrodynamic problems. This algorithm is an extension of the fully automatic (energy-norm based) hp-adaptive strategy described in [7, 18], and a continuation of concepts presented in [14, 20] for elliptic problems. We apply the self-adaptive hp goal-oriented algorithm to accurately simulate induction LWD instruments in a borehole environment with axial symmetry. These instruments are widely used by the geophysical logging industry, and their simulation requires resolution of EM singularities generated by the LWD geometry and rock formation materials [22], as well as resolution of high material constrasts that occur between the mandrel and the borehole. The organization of this document is as follows: In Section 2, we describe the main characteristics of induction logging instruments. We also describe our problem of interest, composed of an induction LWD instrument in a borehole environment, and used for the assessment of the rock formation electrical properties. In Section 3, we introduce Maxwell’s equations, governing the electromagnetic phenomena and explaining the physics of resistivity measurements. We also derive the corresponding variational formulation for axisymmetric problems. A self-adaptive goal-oriented hp algorithm for electrodynamic problems is described in Section 4. The corresponding details of implementation are discussed in the same section. Simulations and numerical results concerning the response of LWD instruments in a borehole environment are shown in Section 5. Section 6 draws the main conclusions, and outline future lines of research. Finally, in the Appendix, we compare numerical results with a semi-analytical solution obtained using Bessel functions for a simplified LWD model problem. The comparison is intended to verify the code as well as to illustrate the high accuracy results obtained with the self-adaptive goal-oriented hp-FEM. 2. Alternate Current (AC) Logging Applications. In this article, we consider an induction1 LWD instrument operating at 2 Mhz. The instrument makes use of one of the following two types of source antennas/coils: • solenoidal coils (Fig. 2.1, left panel), and • toroidal coils (Fig. 2.1, right panel). 2.1. Induction LWD Instruments Based on Solenoidal Coils. For axisymmetric problems, these logging instruments generate a T Mφ field, i.e., the only non-zero components of the electromagnetic (EM) fields are Eφ , Hρ , and Hz , where (ρ, φ, z) denote the cylindrical system of coordinates. A solenoidal coil (Fig. 2.1) produces an impressed current Jimp that we mathematically describe as (2.1)
ˆ Jimp (r) = φIδ(ρ − a)δ(z) ,
where I is the electric current measured in Amperes (A), δ is the Dirac’s delta function, and a is the radius of the solenoid. In the numerical computations, we replace function δ(ρ − a)δ(z) with an approximate function UF that considers the finite dimensions of R the coil, and such that UF dρdz = 1. The analytical electric far-field solution excited by a solenoidal coil of radius a radiating in homogeneous media is given in terms of the electric field by (see [10]) 1 Induction logging instruments are characterized by the fact that impressed current J imp is divergence free (i.e., ∇ · Jimp = 0).
HP -FEM: ELECTROMAGNETIC APPLICATIONS
3
Fig. 2.1. Two coil antennas: a solenoid antenna (left panel) composed of a wire wrapped around a cylinder, and a toroid antenna (right panel) composed of a wire wrapped around a toroid.
(2.2)
2 ˆ E = φωµkIπa
e−jkd j ρ [1 − ] , 4πd kd d
p √ where k = ω 2 ² − jωσ is the wave number, j = −1 is the imaginary unit, ω is angular frequency, ², µ, and σ stand for dielectric permittivity, magnetic permeability, and electrical conductivity of the medium, respectively, and d is the distance between the source coil and the receiver coil. In order to avoid the dependence upon the dimensions of the solenoid, we impose a current on the solenoidal coil equal to 1/(πa2 ) A, i.e., equivalent to that of 1 A with a Vertical Magnetic Dipole (VMD). The corresponding far-field solution in homogeneous media is given by (see [10]) (2.3)
e−jkd j ρ ˆ E = φωµkI [1 − ] . 4πd kd d
Thus, solution (2.3) is independent of the dimensions of the coil2 . 2.2. Induction LWD Instruments Based on Toroidal Coils. For axisymmetric problems, these logging instruments generate a T Eφ field, i.e., the only nonzero components of the EM fields are Hφ , Eρ , and Ez . A toroidal coil induces a magnetic current IM in the azimuthal direction. If we place a toroid of radius a radiating in homogeneous media, the resulting magnetic far-field is given by (see [10]) (2.4)
−jkd j ρ ˆ + jω²)πa2 IM jk e [1 − ] . H = φ(σ 4πd kd d
In order to avoid the dependence upon the dimensions of the toroid, we impose a magnetic current on the toroidal coil equal to that induced by a (σ + jω²) A electric current excitation with a Vertical Electrical Dipole (VED), also known as Hertzian 2 In resistivity logging applications, it is customary to consider solutions that have been divided by the geometrical factor (also called K-factor) [1], so that results are independent (as much as possible) of the logging instrument’s geometry. Thus, solutions obtained from different logging instruments can be readily compared.
4
D. PARDO, L. DEMKOWICZ, C. TORRES-VERD´IN, M. PASZYNSKI
dipole. The corresponding magnetic far-field solution in homogeneous media is given by (see [10]) (2.5)
−jkd j ρ ˆ + jω²)Ijk e [1 − ] . H = φ(σ 4πd kd d
In this case, IM = I/(πa2 ). 2.2.1. Goal of the Computations. We are interested in simulating the EM response of an induction LWD instrument in a borehole environment. For a solenoidal coil, the main objective of our simulation is to compute the first difference of the voltage between the two receiving coils of radius a divided by the (vertical) distance ∆z between them, i.e., ¶ µI I 2πa V1 − V 2 E(l) dl /(∆z) = (2.6) = (E(l1 ) − E(l2 )) , E(l) dl − ∆z ∆z l2 l1 where l1 and l2 are the first and second receiving coils, respectively, and l1 ∈ l1 , l2 ∈ l2 are two arbitrary points located at the receiving coils. Notice that due to the axisymmetry of the electric field, E(lij ) = E(lik ) for all lij , lik ∈ li . This quantity of interest (first difference of voltage) is widely used in resistivity logging applications. Indeed, a first-order asymptotic approximation of the electric field response at low frequencies (Born’s approximation) shows that the voltage at a receiver coil is proportional to the rock formation resistivity in the proximity of such a coil (see [10] for details). At higher frequencies (> 20 Khz), asymptotic approximations (see [1] for details) also indicate the dependence of the voltage upon the rock formation conductivity. Thus, an adequate approximation of the rock formation conductivity (which is unknown a priori in practical applications) can be estimated from the voltage measured at the receiving coils. Computing the first difference of the voltage between two receivers (rather than the voltage at one receiver) is convenient for improving the vertical resolution of the measurements. This well-known fact among well-logging practitioners will be illustrated here with numerical experiments. For a toroidal coil, the main objective of these simulations is to compute the first difference of the electric current at the two receiving coils of radius a divided by the (vertical) distance ∆z between them, i.e., ¶ µI I 2πa I1 − I 2 H(l) dl /(∆z) = H(l) dl − (2.7) = (H(l1 ) − H(l2 )) . ∆z ∆z l2 l1 Notice that the main difference between a toroidal and a solenoidal coil is that the former generates an impressed magnetic current, while the latter produces an impressed electric current. This fact leads to the physical consideration that, if the voltage due to a solenoidal coil is proportional to the rock formation conductivity, then the electric current enforced by a toroidal coil is also proportional to the rock formation resistivity. Thus, the selection of the quantity of interest for toroidal coils (first difference of electric current) is dictated by the physical relation between solenoidal and toroidal coils, and the previous choice of a quantity of interest for solenoidal coils (first difference of voltage). 2.3. Description of a LWD Instrument in a Borehole Environment. We consider a LWD instrument composed of the following axisymmetric materials (all dimensions are given in cm): • one transmitter and two receiver coils defined on
5
HP -FEM: ELECTROMAGNETIC APPLICATIONS
1. ΩC1 = {(ρ, φ, z) : 7.1 < ρ < 7.3 , −2.5 < z < 2.5}, 2. ΩC2 = {(ρ, φ, z) : 7.1 < ρ < 7.3 , 98.75 < z < 101.25}, and, 3. ΩC3 = {(ρ, φ, z) : 7.1 < ρ < 7.3 , 113.75 < z < 116.25}, respectively, • three magnetic buffers with resistivity 104 Ω · m and relative permeability 104 , defined on 1. ΩB1 = {(ρ, φ, z) : 6.675 < ρ < 6.985 , −5 < z < 5}, 2. ΩB2 = {(ρ, φ, z) : 6.675 < ρ < 6.985 , 97.5 < z < 102.5}, and, 3. ΩB3 = {(ρ, φ, z) : 6.675 < ρ < 6.985 , 112.5 < z < 117.5}, respectively, and • a metallic mandrel with resistivity 10−6 Ω · m defined on ΩM = {(ρ, φ, z) : ρ < 7.6}−({(ρ, φ, z) : 6.675 < ρ < 7.6 , −5 < z < 5}∪{(ρ, φ, z) : 6.675 < ρ < 7.6 , 97.5 < z < 102.5} ∪ {(ρ, φ, z) : 6.675 < ρ < 7.6 , 112.5 < z < 117.5}). This LWD instrument is moves along the vertical direction (z-axis) in a subsurface borehole environment composed of: • a borehole mud with resistivity 0.1 Ω · m defined on 1. ΩBH = {(ρ, φ, z) : ρ < 10.795} − (∪i ΩBi ∪ ΩM ), and, • three formation materials of resistivities 100 Ω · m, 10000 Ω · m, and 1 Ω · m, defined on 1. ΩM1 = {(ρ, φ, z) : ρ ≥ 10.795 , (z < −50 or z > 100)}, 2. ΩM2 = {(ρ, φ, z) : ρ ≥ 10.795 , −50 ≤ z < 0}, and, 3. ΩM3 = {(ρ, φ, z) : ρ ≥ 10.795 , 0 ≤ z ≤ 100}, respectively. Fig. 2.2 shows the geometry of the described logging instrument and borehole environment. 6.675 cm
100 Ohm − m
5 cm
Mandrel 0.000001 Ohm m
100 cm
10000 Ohm − m
10 cm
1 Ohm − m
50 cm
Magnetic Buffer 10000 Ohm − m 10000 Relative Permeability
0.1 Ohm − m
15 100 cm cm 0.000001 Ohm − m
100 Ohm − m
Borehole 0.1 Ohm − m Radius = 10.795 cm
Radius 7.6 cm
Fig. 2.2. 2D cross-section of the geometry of an induction LWD problem composed of a metallic mandrel, one transmitter and two receiver coils equipped with magnetic buffers, a borehole, and four layers in the rock formation (with different resistivities). The right panel is an enlarged view of the geometry (left panel) in the vicinity of the transmitter antenna.
3. Maxwell’s Equations. In this section, we first introduce the time-harmonic Maxwell’s equations in the frequency domain. They form a set of first-order Partial Differential Equations (PDE’s). Then, we describe boundary conditions needed for
6
D. PARDO, L. DEMKOWICZ, C. TORRES-VERD´IN, M. PASZYNSKI
the simulation of our logging applications of interest. Finally, we derive a variational formulation in terms of either the electric or the magnetic field, and we reduce the dimension of the computational problem by considering axial symmetry. 3.1. Time-harmonic Maxwell’s equations. Assuming a time-harmonic dependence of the form ejωt , where t denotes time, and ω 6= 0 is angular frequency, Maxwell’s equations can be written as ∇×H = (σ + jω²)E + Jimp Ampere’s Law, ∇×E = −jωµ H − Mimp Faraday’s Law, (3.1) ∇ · (²E) = ρ Gauss’ Law of Electricity, and ∇ · (µH) = 0 Gauss’ Law of Magnetism.
Here H and E denote the magnetic and electric field, respectively, Jimp is a prescribed, impressed electric current density, Mimp is a prescribed, impressed magnetic current density, ², µ, and σ stand for dielectric permittivity, magnetic permeability, and electrical conductivity of the medium, respectively, and ρ denotes the electric charge distribution. We assume µ 6= 0. The equations described in (3.1) are to be understood in the distributional sense, i.e. they are satisfied in the classical sense in subdomains of regular material data, and they also imply appropriate interface conditions across material interfaces. Energy considerations lead to the assumption that the absolute value of both electric field E and magnetic field H must be square integrable. Mimp is assumed to be divergence free due to physical considerations. Maxwell’s equations are not independent. Taking the divergence of Faraday’s Law yields the Gauss’ Law of magnetism. By taking the divergence of Ampere’s Law, and by utilizing Gauss’ Electric Law we arrive at the so called continuity equation, (3.2)
∇ · (σE) + jωρ + ∇ · Jimp = 0 .
3.2. Boundary Conditions (BC’s). There exist a variety of BC’s that can be incorporated into Maxwell’s equations. In the following, we describe those BC’s that are of interest for the logging applications discussed in this paper. At this point, we are considering general 3D domains. A discussion on boundary terms corresponding to the axisymmetry condition is postponed to Section 3.4. 3.2.1. Perfect Electric Conductor (PEC). Maxwell’s equations are to be satisfied in the whole space minus domains occupied by a PEC. A PEC is an idealization of a highly conductive media. Inside a region where σ → ∞, the corresponding electric field converges to zero3 by applying Ampere’s law. Faraday’s law implies that the tangential component of the electric field E must remain continuous across material interfaces in the absence of impressed magnetic surface currents. Consequently, the tangential component of the electric field must vanish along the PEC boundary, i.e., (3.3)
n×E = 0 ,
where n is the unit normal (outward) vector. 3 This result is true under the physical consideration that impressed volume current J imp and σE should remain finite, i.e., hJimp , ψi, hσE, ψi < ∞ for every test function ψ. See [16] for details.
HP -FEM: ELECTROMAGNETIC APPLICATIONS
7
Since the electric field vanishes inside a PEC, Faraday’s law implies that the magnetic field should also vanish inside a PEC in the absence of magnetic currents. The same Faraday’s law implies that the normal component of the magnetic field premultiplied by the permeability must remain continuous across material interfaces. Therefore, the normal component of the magnetic field must vanish along the PEC boundary, i.e., (3.4)
n·H=0.
The tangential component of magnetic field (surface current) and normal component of the electric field (surface charge density) need not be zero, and may be determined a-posteriori. 3.2.2. Source Antennas. Antennas are modeled by prescribing an impressed volume current Jimp . Using the equivalence principle (see, for example, [8]), we can substitute the original impressed electric volume current Jimp by the equivalent electric surface current (3.5)
Jimp = [n×H]S , S
defined on an arbitrary surface S enclosing the support of Jimp , where [n×H]S denotes the jump of n×H accross S. Similarly, an impressed magnetic volume current M imp can be replaced by the equivalent magnetic surface current (3.6)
= −[n×E]S , Mimp S
defined on an arbitrary surface S enclosing the support of Mimp . 3.2.3. Closure of the Domain. We consider a bounded computational domain Ω. A variety of BC’s can be imposed on the boundary ∂Ω such that the difference between solution of such a problem and solution of the original problem defined over R 3 is small. For example, it is possible to use an infinite element technique (as described in [4]). Also, since the electromagnetic fields and their derivatives decay exponentially in the presence of lossy media (non-zero conductivity), we may simply impose a homogeneous Dirichlet or Neumann BC on the boundary of a sufficiently large computational domain. In the field of geophysical logging applications, it is customary to impose a homogeneous Dirichlet BC on the boundary of a large computational domain (for example, 2-20 meters in each direction from a 2 Mhz source antenna in the presence of a resistive media). We will follow the same approach. 3.3. Variational Formulation. From Maxwell’s equations and the BC’s described above, we derive the corresponding standard variational formulation in terms of the electric or magnetic field as follows. First, we notice from Faraday’s law that ∇×E ∈ (L2 (Ω))3 if and only if Mimp ∈ 2 (L (Ω))3 . Since our objective is to find a solution E ∈ H(curl; Ω) = {F ∈ (L2 (Ω))3 : ∇×F ∈ (L2 (Ω))3 }, we shall assume in the case of the electric field formulation (E-formulation) derived below that Mimp ∈ (L2 (Ω))3 . If the prescribed Mimp ∈ / (L2 (Ω))3 , we may still solve Maxwell’s equations with H(curl)-conforming finite elements for the magnetic field by using the H-formulation (3.3.2), or simply by pre˜ imp such that Mimp − M ˜ imp does not radiate outside scribing an equivalent source M the antenna [21]. Similarly, for the H-formulation, we will assume that Jimp ∈ (L2 (Ω))3 .
D. PARDO, L. DEMKOWICZ, C. TORRES-VERD´IN, M. PASZYNSKI
8
3.3.1. E-Formulation . By dividing Faraday’s law by magnetic permeability ¯ where F ∈ H(curl; Ω) is an arbitrary µ, multiplying the resulting equation by ∇×F, test function, and integrating over the domain Ω, we arrive at the identity Z Z Z 1 imp 1 ¯ ¯ ¯ H · (∇×F)dV − (3.7) (∇×E) · (∇×F)dV = −jω M · (∇×F)dV. µ µ Ω Ω Ω Z ¯ dV by parts, and applying Ampere’s law, we obtain Integrating H · (∇×F) Ω
(3.8)
Z
Z Z ¯ ¯ ¯ t dS = H · (∇×F) dV = (∇×H) · F dV − [n×H]ΓN · F ΓN Ω Ω Z Z Z ¯ dV + ¯ dV − ¯ t dS, (σ + jω²)E · F Jimp · F [n×H]ΓN · F Ω
Ω
ΓN
¯ denoting a surface contained in the closure of Ω where an impressed with ΓN ⊂ Ω electric surface current Jimp ΓN may be prescribed. Ft = F − (F · n) · n is the tangential component of vector F on ΓN , and n is the unit normal outward (with respect Ω if ΓN ⊂ ∂Ω) vector. Substitution of (3.8) into (3.7), and use of equation (3.5) yields to the following variational identity, valid for any test function F ∈ H(curl; Ω): Z Z Z 1 ¯ dV ¯ dV = −jω ¯ Jimp · F k2 E · F (∇×E) · (∇×F)dV − Ω Ω Ω µZ Z (3.9) 1 imp ¯ t dS − ¯ dV , · F Jimp +jω M · (∇×F) ΓN Ω µ ΓN where k 2 = ω 2 ² − jωσ.
Finally, in order to obtain a unique solution E ∈ H(curl; Ω) for problem (3.9), we introduce a Dirichlet boundary condition on a part ΓD of the boundary ∂Ω of the computational domain Ω. Thus, we obtain the following variational formulation: E ∈ ED + HD (curl; Ω) such that: Find Z Z Z 1 ¯ dV ¯ dV = −jω ¯ dV − Jimp · F k2 E · F (∇×E) · (∇×F) (3.10) µ Ω Ω Ω Z Z 1 imp imp ¯ ¯ dV ∀ F ∈ HD (curl; Ω) , · F dS − J +jω M · (∇×F) t ΓN µ Ω ΓN
where ED is a lift (typically ED = 0) of the essential boundary condition data ED (denoted with the same symbol), and HD (curl; Ω) = {F ∈ H(curl; Ω) : (n×F)|ΓD = 0} is the space of admissible test functions associated with problem (3.10). Conversely, we can derive (3.1), (3.3), and (3.5) from variational problem (3.10). Remark 1. At this point, an impressed magnetic surface current Mimp defined S on a subset of ΓD may be introduced into the formulation by using equation (3.6). It follows that ED = n×E|ΓD = −Mimp ΓD . 3.3.2. H-Formulation . By dividing Ampere’s law by σ + jω², multiplying the ¯ where F ∈ H(curl; Ω) is an arbitrary test function, and resulting equation by ∇×F, integrating over the domain Ω, we arrive at the identity
HP -FEM: ELECTROMAGNETIC APPLICATIONS
(3.11)
−jω
Integrating
(3.12)
Z
Z
Z
Ω
Ω
9
Z 1 ¯ ¯ dV (∇×H) · (∇× F)dV = E · (∇×F) k2 Ω Z 1 imp ¯ dV . −jω J · (∇×F) 2 Ω k
¯ dV by parts, and applying Faraday’s law, we obtain E · (∇×F)
Z Z ¯ t dS = ¯ ¯ [n×E]ΓN · F E · (∇×F) dV = (∇×E) · F dV − ΓN Ω Ω Z Z Z ¯ dV − ¯ dV − ¯ t dS, −jω µH · F Mimp · F [n×E]ΓN · F Ω
Ω
ΓN
¯ denoting a surface contained in the closure of Ω where an impressed with ΓN ⊂ Ω magnetic surface current Mimp ΓN may be prescribed. Substitution of (3.12) into (3.11), and use of equation (3.6) yields the following variational identity, valid for any test function F ∈ H(curl; Ω): Z Z 1 ¯ dV = ¯ µH · F (∇×H) · (∇×F)dV + jω −jω 2 Ω Ω k Z Z Z (3.13) 1 imp ¯ dV . ¯ dV + ¯ J · (∇×F) − Mimp · F Mimp ΓN · Ft dS − jω 2 Ω ΓN Ω k Finally, in order to obtain a unique solution H ∈ H(curl; Ω) for problem (3.13) we introduce a Dirichlet boundary condition on a part ΓD of the boundary ∂Ω of the computational domain Ω. Thus, we obtain the following variational formulation: Find H ∈ HD + HD (curl; Ω) such that: Z Z Z 1 ¯ ¯ ¯ (∇×H) · (∇×F)dV + jω µH · FdV =− Mimp · FdV (3.14) σ + jω² Ω Ω ΩZ Z 1 ¯ ¯ Jimp · (∇×F)dV ∀ F ∈ HD (curl; Ω), Mimp + ΓN · Ft dS + Ω σ + jω² ΓN where HD is a lift (typically HD = 0) of the essential boundary condition data HD (denoted with the same symbol). At this point, an impressed electric surface current Jimp defined on a subset of ΓD may be introduced into the formulation by using S equation (3.5). It follows that HD = n×H|ΓD = Jimp ΓD .
3.4. Cylindrical Coordinates and Axisymmetric Problems . We consider cylindrical coordinates (ρ, φ, z). For the geophysical logging applications considered in this article, we assume that both the logging instrument and the rock formation properties are axisymmetric (invariant with respect to the azimuthal component φ) around the axis of the borehole. Under this assumption, we obtain that for any vector ˆ φ + zˆAz , field A = ρˆAρ + φA (3.15)
∇×A = −ˆ ρ
∂Aφ ˆ ∂Aρ − ∂Az ) + zˆ 1 ∂(ρAφ ) . + φ( ∂z ∂z ∂ρ ρ ∂ρ
3.4.1. E-Formulation. Next, we consider the space of all test functions F ∈ HD (curl; Ω) such that F = (0, Fφ , 0). According to (3.15), (3.16)
∇×F = −ˆ ρ
∂Fφ 1 ∂(ρFφ ) + zˆ . ∂z ρ ∂ρ
10
D. PARDO, L. DEMKOWICZ, C. TORRES-VERD´IN, M. PASZYNSKI
Variational formulation (3.10) reduces to a formulation in terms of the scalar field Eφ , namely,
(3.17)
˜ 1 (Ω) such that: Find Eφ ∈ Eφ,D + H D µ ¶ Z Z ¯ 1 ∂Eφ ∂ Fφ 1 ∂(ρEφ ) ∂(ρF¯φ ) + k 2 Eφ F¯φ dV = dV − Ωµ ∂z ∂z ρ2 ∂ρ ∂ρ Ω Z Z imp ¯ imp ¯ Jφ Fφ dV + jω −jω Jφ,ΓN Fφ dS ΓN Ω · Z ¯ ¯ ¸ 1 imp 1 ∂(ρFφ ) imp ∂ Fφ 1 ˜D − + Mz −Mρ dV ∀ Fφ ∈ H (Ω) , µ ∂z ρ ∂ρ Ω
˜ 1 (Ω) = {Eφ : (0, Eφ , 0) ∈ HD (curl; Ω)} = {Eφ ∈ L2 (Ω) : 1 Eφ + ∂Eφ ∈ where H D ρ ∂ρ ∂E φ 2 2 ∈ L (Ω), Eφ |ΓD = 0}. Similarly, for a test function F = (Fρ , 0, Fz ), L (Ω) , ∂z variational problem (3.10) simplifies to: ˜ D (curl; Ω) such that: Find E = (Eρ , 0, Ez ) ∈ ED + H µ ¶µ ¯ ¶ Z Z 1 ∂Eρ ∂ Fρ ∂Ez ∂ F¯z dV − k 2 (Eρ F¯ρ + Ez F¯z ) dV = − − Ω µ ∂z ∂ρ ∂z ∂ρ Ω Z Z (3.18) imp ¯ imp ¯ imp ¯ imp ¯ Jρ,Γ Fρ + Jz,Γ Fz dS F dV + jω F + J J −jω z ρ z ρ N N Γ Ω N · ¸ Z ∂ F¯z 1 imp ∂ F¯ρ ˜ D (curl; Ω) , − Mφ − dV ∀ F = (Fρ , 0, Fz ) ∈ H ∂z ∂ρ Ω µ
˜ D (curl; Ω) = {(Eρ , Ez ) : E = (Eρ , 0 , Ez ) ∈ L2 (Ω) , (∇×E)|φ = ∂Eρ − where H ∂z ∂Ez ∈ L2 (Ω) , (n×E)|ΓD = 0}. ∂ρ In summary, problem (3.10) decouples into a system of two simpler problems described by (3.17) and (3.18). 1 ˜D It has been shown in [2] (Lemma 4.9) that space H (Ω) can also be 1 ˜ 1 (Ω) = {Eφ ∈ L2 (Ω) : Eφ ∈ L2 (Ω) , ∇(ρ,z) Eφ ∈ L2 (Ω)}. expressed as H D ρ
Remark 2.
3.4.2. H-Formulation. Using the same decomposition of test functions (i.e., F = (0, Fφ , 0) , and F = (Fρ , 0, Fz )) for variational problem (3.14), we arrive at the following two decoupled variational problems in terms of (0, Hφ , 0) (3.19), and (Hρ , 0, Hz ) (3.20), respectively:
(3.19)
˜ 1 (Ω) such that: Find Hφ ∈ Hφ,D + H D µ ¶ Z ¯φ 1 ∂ F 1 ∂(ρHφ ) ∂(ρF¯φ ) ∂H φ + 2 dV Ω σ + jω² ∂z ∂z ρ ∂ρ ∂ρ Z Z Z imp ¯ Fφ dS Mφ,Γ Mφimp F¯φ dV + µHφ F¯φ dV = − +jω N Γ Ω Ω N · ¸ Z 1 1 ∂(ρF¯φ ) ∂ F¯φ 1 ˜D + + Jzimp −Jρimp dV ∀ Fφ ∈ H (Ω) . ∂z ρ ∂ρ Ω σ + jω²
11
HP -FEM: ELECTROMAGNETIC APPLICATIONS
˜ D (curl; Ω) such that: Find H = (Hρ , 0, Hz ) ∈ HD + H µ ¶ µ ¶ Z 1 ∂Hz ∂ F¯z ∂Hρ ∂ F¯ρ − − dV ∂z ∂ρ ∂z ∂ρ Ω σ Z+ jω² Z Mρimp F¯ρ + Mzimp F¯z dV µ(Hρ F¯ρ + Hz F¯z ) dV = − +jω Ω Ω · ¯ ¸ Z Z 1 ∂ F¯z imp ¯ imp ¯ imp ∂ Fρ Mρ,ΓN Fρ + Mz,ΓN Fz dS + dV Jφ − + ∂z ∂ρ ΓN Ω σ + jω² ˜ D (curl; Ω) . ∀ F = (Fρ , 0, Fz ) ∈ H
(3.20)
From the formulation of problems (3.17) trough (3.20), we remark the following: • Physically, solution of problems (3.18), and (3.19) correspond to the T E φ mode (i.e. Eφ = 0), and solution of problems (3.17), and (3.20) correspond to the T Mφ -mode (i.e. Hφ = 0). • The axis of symmetry is not a boundary of the original 3D problem, and therefore, a boundary condition should not be needed to solve this problem. Nevertheless, formulations of problems (3.17) through (3.20) require the use ˜ 1 (Ω) and H ˜ D (curl; Ω) described above. The former space inof spaces H D volves the singular weight ρ1 , which implicitly requires a homogeneous Dirichlet boundary condition along the axis of symmetry. The latter space can be considered as it is (by using 2D edge elements), and no BC is necessary 4 to solve the problem. 4. Self-Adaptive Goal-Oriented hp-FEM . We are interested in solving variational problems (3.10) and (3.14) (or alternatively, (3.17), (3.18), (3.19), and (3.20)), that we state here in terms of sesquilinear form b, and antilinear form f : ( Find E ∈ ED + V (4.1) b(E, F) = f (F) ∀F ∈ V , where • • • •
ED is a lift of the essential (Dirichlet) BC. V is a Hilbert space. f ∈ V0 is an antilinear and continuous functional on V. b is a sesquilinear form. More precisely, we have: 1 a(E, F) − k 2 c(E, F) E-Formulation µ (4.2) b(E, F) = 1 a(E, F) − µc(E, F) H-Formulation k2
,
where sesquilinear forms a and c are assumed to be Hermitian, continuous and V-coercive. We define an “energy” inner product on V as: 1 a(E, F) + |k 2 |c(E, F) E-Formulation µ (E, F) := (4.3) , 1 2 a(E, F) + µc(E, F) H-Formulation |k | with the corresponding (energy) norm denoted by kEk.
4 From the computational point of view, this effect can be achieved by artificially adding a homogeneous natural (Neumann) BC.
12
D. PARDO, L. DEMKOWICZ, C. TORRES-VERD´IN, M. PASZYNSKI
4.1. Representation of the Error in the Quantity of Interest. Given an hp-FE subspace Vhp ⊂ V, we discretize (4.1) as follows: ( Find Ehp ∈ ED + Vhp (4.4) b(Ehp , Fhp ) = f (Fhp ) ∀Fhp ∈ Vhp . The objective of goal-oriented adaptivity is to construct an optimal hp-grid, in the sense that it minimizes the problem size needed to achieve a given tolerance error for a given quantity of interest L, with L denoting a linear and continuous functional. By recalling the linearity of L, we have: (4.5)
Error of interest = L(E) − L(Ehp ) = L(E − Ehp ) = L(e) ,
where e = E − Ehp denotes the error function. By defining the residual rhp ∈ V0 as rhp (F) = f (F) − b(Ehp , F) = b(E − Ehp , F) = b(e, F), we look for the solution of the dual problem: ( ¯ ∈V Find W (4.6) b(F, W) = L(F) ∀F ∈ V . Using the Lax-Milgram theorem we conclude that problem (4.6) has a unique solution in V. The solution W, is usually referred to as the influence function. By discretizing (4.6) via, for example, Vhp ⊂ V, we obtain: ( ¯ hp ∈ Vhp Find W (4.7) b(Fhp , Whp ) = L(Fhp ) ∀Fhp ∈ Vhp . Definition of the dual problem plus the Galerkin orthogonality for the original problem imply the final representation formula for the error in the quantity of interest, namely, L(e) = b(e, W) = b(e, W − Fhp ) = ˜b(e, ²) . | {z } ²
At this point, Fhp ∈ Vhp is arbitrary, and ˜b(e, ²) = b(e, ²¯) denotes the bilinear form corresponding to the original sesquilinear form. Notice that, in practice, the dual problem is solved not for W but for its complex ¯ utilizing the bilinear form and not the sesquilinear form. The linear conjugate W system of equations is factorized only once, and the extra cost of solving (4.7) reduces to only one backward and one forward substitution (if a direct solver is used). Once the error in the quantity of interest has been determined in terms of bilinear form ˜b, we wish to obtain a sharp upper bound for |L(e)| that depends upon the mesh parameters (element size h and order of approximation p) only locally. Then, a selfadaptive algorithm intended to minimize this bound will be defined. First, using a procedure similar to the one described in [7], we approximate E and W with fine grid functions E h , p+1 , W h , p+1 , which have been obtained by 2 2 solving the corresponding linear system of equations associated with the FE subspace V h , p+1 . In the remainder of this article, E and W will denote the fine grid solutions 2 of the direct and dual problems (E = E h , p+1 , and W = W h , p+1 , respectively), and 2 2 we will restrict ourselves to discrete FE spaces only. Next, we bound the error in the quantity of interest by a sum of element contributions. Let bK denote a contribution from element K to sesquilinear form b. It then follows that
HP -FEM: ELECTROMAGNETIC APPLICATIONS
(4.8)
|L(e)| = |b(e, ²)| ≤
X K
13
|bK (e, ²)| ,
where summation over K indicates summation over elements. 4.2. Projection based interpolation operator. Once we have a representation formula for the error in the quantity of interest in terms of the sum of element contributions given by (4.8), we wish to express this upper bound in terms of local quantities, i.e. in terms of quantities that do not vary globally when we modify the grid locally. For this purpose, we introduce the idea of projection-based interpolation operators. First, in order to simplify the notation, we define the following three spaces of admissible solutions: • V = HD (curl; Ω), ˜ D (curl; Ω), and, • V2D = H ˜ 1 (Ω). • V 1D = H D 2D 1D The corresponding hp-Finite Element spaces will be denoted by Vhp , Vhp , and Vhp , respectively. At this point, we introduce three projection-based interpolation operators that have been defined in [6, 5], and used in [7, 18] for the construction of the fully automatic energy-norm based hp-adaptive algorithm: curl,3D • Πhp : V −→ Vhp ,
curl,2D 2D : V2D −→ Vhp , and, • Πhp 1D 1D 1D −→ Vhp . • Πhp : V We shall also consider three projection operators in the energy-norm: curl,3D • Php : V −→ Vhp ,
curl,2D 2D : V2D −→ Vhp • Php , and, 1D 1D 1D −→ Vhp . • Php : V To further simplify the notation, we will utilize the unique symbol Πcurl hp to denote all projection based interpolation operators mentioned above. Depending upon the problem formulation (and corresponding space of admissible solutions), Π curl hp should curl,2D curl,3D for problems (3.18) for problems (3.10) and (3.14), Πhp be understood as Πhp and (3.20), or Π1D hp for problems (3.17) and (3.19). Similarly, we will use the unique curl,2D curl,3D 1D curl . , or Php , Php symbol Php to denote either Php
We denote Ehp = Pcurl hp E. Equation (4.8) then becomes X X curl curl |L(e)| ≤ |bK (E, ²)| = |bK (E − Πcurl (4.9) hp E, ²) + bK (Πhp E − Php E, ²)| . K
K
Given an element K, it is expected that |bK (Πcurl hp E − Php E, ²)| will be negligible compared to |bK (E − Πcurl E, ²)|. Under this assumption, we conclude that: hp X |bK (E − Πcurl (4.10) |L(e)| ¹ hp E, ²)| . K
In particular, for ² = W − Πcurl hp W, we have: X curl (4.11) |L(e)| ¹ |bK (E − Πcurl hp E, W − Πhp W)| . K
D. PARDO, L. DEMKOWICZ, C. TORRES-VERD´IN, M. PASZYNSKI
14
By applying Cauchy-Schwartz inequality, we obtain the next upper bound for |L(e)|: (4.12)
|L(e)| ¹
X K
k˜ ekK k˜²kK ,
˜ = E − Πcurl where e ² = W − Πcurl hp E, ˜ hp W, and k · kK denotes energy-norm k · k restricted to element K. 4.3. Fully Automatic Goal-Oriented hp-Refinement Algorithm. We describe an hp self-adaptive algorithm that utilizes the main ideas of the fully automatic (energy-norm based) hp-adaptive algorithm described in [7, 18]. We start by recalling the main objective of the self-adaptive (energy-norm based) hp-refinement strategy, which consists of solving the following maximization problem: ˜ Find an optimal hp-grid in the following sense: 2 curl 2 X kE − Πcurl (4.13) hp EkK − kE − Πhp c EkK ˜ , hp = arg max ∆N c hp K
where • E = E h , p+1 is the fine grid solution, and 2
c • ∆N > 0 is the increment in the number of unknowns from grid hp to grid hp. Similarly, for goal-oriented hp-adaptivity, we propose the following algorithm based on estimate (4.12): ˜ Find an optimal hp-grid in the following sense: " curl X kE − Πcurl hp EkK · kW − Πhp WkK hp ˜ = arg max ∆N c (4.14) hp K # curl kE − Π c EkK · kW − Πcurl c WkK hp hp − , ∆N
where: • E = E h , p+1 and W = W h , p+1 are the fine grid solutions corresponding to 2 2 the direct and dual problems, and c • ∆N > 0 is the increment in the number of unknowns from grid hp to grid hp. Implementation of the goal-oriented hp-adaptive algorithm is based on the optimization procedure used for energy-norm hp-adaptivity [7, 18]. 4.4. Implementation details. In what follows, we discuss the main implementation details needed to extend the fully automatic (energy-norm based) hp-adaptive algorithm [7, 18] to a fully automatic goal-oriented hp-adaptive algorithm. 1. First, the solution W of the dual problem on the fine grid is necessary. This goal can be attained either by using a direct (frontal) solver or an iterative (two-grid) solver (see [13]). 2. Subsequently, we should treat both solutions as satisfying two different partial differential equations (PDE’s). We select functions E and W as the solutions of the system of two PDE’s. 3. We proceed to redefine the evaluation of the error. The energy-norm error evaluation of a two dimensional function is replaced by the product k E − curl Πcurl hp E k · k W − Πhp W k.
15
HP -FEM: ELECTROMAGNETIC APPLICATIONS
4. After these simple modifications, the energy-norm based self-adaptive algorithm may now be utilized as a self-adaptive goal-oriented hp algorithm. 5. Numerical Results. In this Section, we apply the goal-oriented hp selfadaptive strategy described in Section 4 to simulate the response of the induction LWD instrument operating at 2 Mhz considered in Section 2.3, using formulation (3.17) for solenoidal coils, and (3.19) for toroidal coils. Exactly the same results are obtained with formulations (3.20) and (3.18), respectively, as predicted by the theory. Thus, formulations (3.20) and (3.18) have been used as an extra verification of the simulations, and the corresponding results have been omitted in this article to avoid duplicity. Fig. 5.1 displays the first vertical difference of the electric field (divided by the distance between the two receivers) for the described LWD instrument equipped with solenoidal coils. The three curves correspond to: 1. the rock formation with no mud-filtrate invasion, 2. the rock formation with a 2 Ω·m 40 cm. horizontal mud layer invading the 1 Ω·m rock formation layer, and a 5Ω·m 90cm. horizontal mud layer invading the 10000 Ω·m rock formation layer, and 3. the previous (mud invaded) rock formation, using a mandrel with relative magnetic permeability of 100. For toroidal antennas, we display in Fig. 5.2 the first vertical difference of the magnetic field (divided by the distance between the two receivers). The three displayed curves correspond to the three situations discussed above. Invasion Study
Vertical Position of Receiving Antenna (m)
2.5 2
Solenoid
100 Ohm−m
1.5 1 0.52 Ohm−m (0.4m)
1 Ohm−m
0
5 Ohm−m (0.9m)
10000 Ohm−m
−0.5 −1
3
No Invasion 0.4/0.9 m Invasion 0.4/0.9m Invasion, Perm. Mandrel=100
100 Ohm−m
−1.5 0.04 0.08 0.12 0.16 0.2 Amplitude First Vert. Diff. Electric Field (V/m2)
2.5 Vertical Position of Receiving Antenna (m)
3
No Invasion 0.4/0.9 m Invasion 0.4/0.9m Invasion, Perm. Mandrel=100
100 Ohm−m
2 1.5 1
0.52 Ohm−m (0.4m)
1 Ohm−m
0
5 Ohm−m (0.9m)
10000 Ohm−m
−0.5
100 Ohm−m
−1 −1.5 −200
−100
0 100 Phase (degrees)
200
Fig. 5.1. LWD problem equipped with a solenoidal source. Amplitude (left panel) and phase (right panel) of the first vertical difference of the electric field (divided by the distance between receivers) at the receiving coils. Results obtained with the self-adaptive goal-oriented hp-FEM. The spatial distribution of electrical resistivity is also displayed to facilitate the physical interpretation of results.
D. PARDO, L. DEMKOWICZ, C. TORRES-VERD´IN, M. PASZYNSKI
16
Invasion Study
Toroid
3
2
No Invasion 0.4/0.9 m Invasion 0.4/0.9m Invasion, Perm. Mandrel=100
100 Ohm−m
1.5 1
1 Ohm−m
0.52 Ohm−m (0.4m) 0
5 Ohm−m (0.9m)
10000 Ohm−m
−0.5 −1
100 Ohm−m
−1.5 −5 −3 −1 10 10 10 2 Amplitude First Vert. Diff. Magnetic Field (A/m )
2.5 Vertical Position of Receiving Antenna (m)
Vertical Position of Receiving Antenna (m)
2.5
3 No Invasion 0.4/0.9 m Invasion 0.4/0.9m Invasion, Perm. Mandrel=100
100 Ohm−m
2 1.5 1
0.52 Ohm−m (0.4m)
1 Ohm−m
0
5 Ohm−m (0.9m)
10000 Ohm−m
−0.5
100 Ohm−m
−1 −1.5 −200
−100
0 100 Phase (degrees)
200
Fig. 5.2. LWD problem equipped with a toroidal source. Amplitude (left panel) and phase (right panel) of the first vertical difference of the magnetic field (divided by the distance between receivers) at the receiving coils. Results obtained with the self-adaptive goal-oriented hp-FEM. The spatial distribution of electrical resistivity is also displayed to facilitate the physical interpretation of results.
These results illustrate the strong dependence of the LWD response on the rock formation resistivity. We observe that solenoidal antennas are more sensitive to highly conductive formations as well as to the electrical permeability of the mandrel, while toroidal antennas are more sensitive to highly resistive formations. Fig. 5.3 illustrates the effect of the magnetic buffers. By removing the magnetic buffers from the logging instrument’s design, the amplitude of the received signal decreases by a factor of up to 200 in the case of a solenoidal source. For practical applications, a strong signal on the receivers is desired to minimize the noise-to-signal ratio. Thus, it is appropriate to use magnetic buffers in combination with solenoidal antennas. On the contrary, the use of magnetic buffers with toroidal antennas is not advisable since they weaken the received signal. In both cases, the phase and shape of the solution is not sensitive to the presence (or not) of magnetic buffers, and the corresponding results have been omitted. The exponential convergence obtained using the self-adaptive goal-oriented hpFEM is shown in Fig. 5.4 (left panel), by considering an arbitrary fixed position of the logging instrument for a solenoid antenna. The final grid delivers a relative error in the quantity of interest below 0.00001%, i.e. the first 7 significant digits of the quantity of interest are exact. In Fig. 5.4 (right panel), we display the exponential convergence of the energy-norm based hp-FEM. The final hp-grid delivers an energynorm error below 0.01%. Nevertheless, the quantity of interest still contains a relative error above 15%. A final goal-oriented hp-grid delivering a relative error in the quantity of interest of 0.1% is displayed in Fig. 5.5.
17
HP -FEM: ELECTROMAGNETIC APPLICATIONS
Solenoid
Toroid
100 Ohm−m
1.5 1 0.5
1 Ohm−m
0
10000 Ohm−m −0.5 −1
100 Ohm−m
−1.5 −4 −2 0 10 10 10 Amplitude First Vert. Diff. of Electric Field (V/m2)
3
3 WithWith Magnetic Magnetic Buffers Buffers Without Without Magnetic Magnetic Buffers Buffers
2.5 2.5 Vertical Position of Receiving Antenna (m)
2
With Magnetic Buffers Without Magnetic Buffers Vertical Position of Receiving Antenna (m)
Vertical Position of Receiving Antenna (m)
2.5
3
2
100100 Ohm−m Ohm−m
2
1.5 1.5 1
1
1 Ohm−m 1 Ohm−m
0.5 0.5 0
0
10000 10000 Ohm−m Ohm−m
−0.5−0.5 −1 −1 −1.5−1.5 −5 −200 10 −100
100100 Ohm−m Ohm−m
0 100 100 200 105 Phase (degrees) Amplitude First Vert. Diff. of Magnetic Field (A/m2)
With Magnetic Buffers Without Magnetic Buffers
2.5 Vertical Position of Receiving Antenna (m)
3
100 Ohm−m
2 1.5 1
1 Ohm−m
0.5 0
10000 Ohm−m −0.5
100 Ohm−m
−1 −1.5 −200
−100
Fig. 5.3. LWD problem equipped with a solenoidal source. Results obtained with the selfadaptive goal-oriented hp-FEM correspond to the use of solenoidal antennas (left panel), and toroidal antennas (right panel), respectively. The spatial distribution of electrical resistivity is also displayed to facilitate the physical interpretation of results.
6. Summary and Conclusions. We have successfully applied a self-adaptive goal-oriented hp-FE algorithm to simulate the axisymmetric response of an induction LWD instrument in a borehole environment. These simulations would not be possible with energy-norm adaptive algorithms. Numerical results illustrate the exponential convergence of the method (allowing for high accuracy simulations), the suitability of the presented formulations for axisymmetric electrodynamic problems, and the main physical characteristics of the presented induction LWD instrument. These results suggest the use of solenoidal antennas for the assessment of highly conductive rock formation materials, and toroidal antennas for the assessment of highly resistive materials. Solenoidal antennas should be used in combination with magnetic buffers to strengthen the measured EM signal, while the use of magnetic buffers with toroidal antennas should be avoided. Both types of antennas can be used to study mud-filtrate invasion. Since the influence function used by the self-adaptive goal-oriented hp-adaptive algorithm is approximated via Finite Elements, the numerical method presented in this article is problem independent, and it can be applied to 1D, 2D, and 3D FE discretizations of H 1 -, H(curl)-, and H(div)-spaces. Appendix. A Loop-Antenna Radiating in a Homogeneous Lossy Medium in the Presence of a Highly Conductive Metallic Mandrel. In this appendix, we consider a problem with known analytical solution. We use this problem as an additional mechanism to verify the code, as well as to provide comparative results between analytical and numerical solutions. We consider a solenoid (or a toroid) of radius a radiating at a frequency of 2 Mhz
0 100 Phase (degrees)
200
D. PARDO, L. DEMKOWICZ, C. TORRES-VERD´IN, M. PASZYNSKI
18
Goal−Oriented hp−Adaptivity
3
10 Upper bound for |L(e)|/|L(u)| |L(e)|/|L(u)|
2
10
1
1
10
10
0
Relative Error in %
Relative Error in %
Energy−norm error |L(e)|/|L(u)|
2
10
10
−1
10
−2
10
−3
−1
10
−2
10
10
−4
−4
10
10
−5
−5
0
0
10
−3
10
10
Energy−norm hp−Adaptivity
3
10
1000 8000 27000 64000 Number of Unknowns N (scale N1/3)
10
0
1000 8000 27000 64000 Number of Unknowns N (scale N1/3)
Fig. 5.4. LWD problem equipped with a solenoidal source. Left panel: convergence behavior obtained with the self-adaptive goal-oriented hp-FEM shows exponential convergence rates for estimate (4.8) (solid curve) used for optimization. The dashed curve describes the relative error in the quantity of interest. Right panel: convergence behavior obtained with the self-adaptive energy-norm hp-FEM shows exponential convergence rates for the energy-norm. The dashed curve describes the relative error in the quantity of interest.
2Dhp90: D A Fully automatic m hp-adaptive Finite EElement m code
p=8 p=7 Rece ver II −→
p=6
Rece ver I −→
p=5 p=4 p=3 p=2
Transm tter −→
p=1 y x x 200 cm) o F g 5 5 LWD ns rumen equ pped w h a so eno da source Por on (120z cm he fina hp-gr d D fferen co ors nd ca e d fferen po ynom a orders o approx ma on rang ng rom 1 ( gh grey) o 8 (wh e)
HP -FEM: ELECTROMAGNETIC APPLICATIONS
19
in a homogeneous lossy medium (with resistivity equal to 1 Ω · m), in the presence of an infinitely large cylindrical mandrel (with resistivity equal to 10−6 Ω · m) of radius b < a. The coil and the mandrel exhibit axial symmetry (see Fig. A.1). z
COIL (Toroid/Solenoid)
b
a
MANDREL
Fig. A.1. Geometry of a loop-antenna radiating in a homogeneous lossy medium in the presence of a highly conductive metallic mandrel.
For a solenoidal coil located at z = 0, the resulting solution for a ≤ ρ ≤ b is given by [10, 15]: Z −ωµ ∞ (1) (1) (A.1) Eφ (ρ, z) = [J1 (kρ a) + ΓH1 (kρ a)]H1 (kρ ρ)eikz z dkρ , 4πa −∞ (1)
(1)
where Γ = −J1 (kρ b)/H1 (kρ b), Jp and Hp
are qthe Bessel and Hankel functions of the first type of order p, respectively, and kz = k 2 − kρ2 . For a toroidal coil located at z = 0, the resulting solution for a ≤ ρ ≤ b is given by: Z ∞ −i (1) (1) (A.2) Hφ (ρ, z) = [J1 (kρ a) + ΓH1 (kρ a)]H1 (kρ ρ)eikz z dkρ , 4πa −∞ (1)
where Γ = −J0 (kρ b)/H0 (kρ b). In figures A.2 and A.3, we display a comparison between analytical and numerical results (obtained using the self-adaptive hp goal-oriented algorithm) for the solenoidal and toroidal coils, respectively. We selected b = 0.0254 cm, and a = 0.03048 cm. The numerical results accurately reproduce the analytical ones, both in terms of amplitude and phase. When considering a solenoid, the logging instrument response using a mandrel of resistivity 10−5 Ω · m or a PEC mandrel are indistinguishable in terms of amplitude. A similar situation occurs for a toroid. In terms of phase, induction instruments equipped with solenoidal coils appear to be more sensitive to the mandrel resistivity than those equipped with toroidal coils. Acknowledgments. This work was financially supported by Baker-Atlas and the Joint Industry Research Consortium on Formation Evaluation supervised by Prof. C. Torres-Verdin. We would also like to acknowledge the expertise and technical advise
20
D. PARDO, L. DEMKOWICZ, C. TORRES-VERD´IN, M. PASZYNSKI
3.5
3.5 Analytical solution −7 Mandrel Resisitivity: 10
Analytical solution −7 Mandrel Resisitivity: 10
−5
Mandrel Resisitivity: 10
2.5
2
1.5
1
0.5 −8 10
Mandrel Resisitivity: 10−3 Mandrel Resisitivity: 1
3 Distance in z−axis from transmitter to receiver (in m)
Distance in z−axis from transmitter to receiver (in m)
−5
Mandrel Resisitivity: 10
Mandrel Resisitivity: 10−3 Mandrel Resisitivity: 1
3
2.5
2
1.5
1
−6
10
−4
−2
10 10 Amplitude (V/m)
0
10
0.5 −180
−90
0 Phase (degrees)
90
180
Fig. A.2. Solution (electric field) along the vertical axis passing through a solenoid radiating in a homogeneous medium in the presence of a metallic mandrel. Analytical solution (mandrel is a PEC) against the numerical solution for different mandrel resistivities (10 −7 , 10−5 , 10−3 , and 1 Ω · m) obtained with the self-adaptive goal oriented hp-FEM.
received from L. Tabarovsky, A. Bespalov, T. Wang, and other members of the Science Department of Baker-Atlas. REFERENCES [1] B. I. Anderson, Modeling and Inversion Methods for the Interpretation of Resistivity Logging Tool Response, PhD thesis, Delft University of Technology, 2001. [2] F. Assous, C.(Jr.) Ciarlet, and S Labrunie, Theoretical tools to solve the axisymmetric Maxwell equations., Math. Meth. Appl. Sci., 25 (2002), pp. 49–78. [3] R. Becker and R. Rannacher, Weighted a posteriori error control in FE methods., in ENUMATH 97. Proceedings of the 2nd European conference on numerical mathematics and advanced applications held in Heidelberg, Germany, September 28-October 3, 1997. Including a selection of papers from the 1st conference (ENUMATH 95) held in Paris, France, September 1995. Singapore: World Scientific. 621-637 , Bock, Hans Georg (ed.) et al., 1998. [4] W. Cecot, W. Rachowicz, and L. Demkowicz, An hp-adaptive finite element method for electromagnetics. III: A three-dimensional infinite element for Maxwell’s equations., Int. J. Numer. Methods Eng., 57 (2003), pp. 899–921. [5] L. Demkowicz, Finite element methods for Maxwell equations., Encyclopedia of Computational Mechanics, (eds. E. Stein, R. de Borst, T.J.R. Hughes), Wiley and Sons, 2003, in review, (2003). [6] L. Demkowicz and A. Buffa, H 1 , H(curl), and H(div) conforming projection-based interpolation in three dimensions: quasi optimal p-interpolation estimates., Comput. Methods Appl. Mech. Eng., 194 (2005), pp. 267–296. [7] L. Demkowicz, W. Rachowicz, and Ph. Devloo, A fully automatic hp-adaptivity., J. Sci. Comput., 17 (2002), pp. 117–142. [8] R. F. Harrington, Time-harmonic electromagnetic fields., McGraw-Hill, New York, 1961. [9] V. Heuveline and R. Rannacher, Duality-based adaptivity in the hp-finite element method., J. Numer. Math., 11 (2003), pp. 95–113.
21
HP -FEM: ELECTROMAGNETIC APPLICATIONS
3.5
3.5 Analytical solution
Analytical solution
−7
−7
Mandrel Resisitivity: 10
−5
Mandrel Resisitivity: 10
Mandrel Resisitivity: 10
−5
Mandrel Resisitivity: 10
−3
2.5
2
1.5
1
0.5 −10 10
Mandrel Resisitivity: 10−3 Mandrel Resisitivity: 1
3
Mandrel Resisitivity: 10 Mandrel Resisitivity: 1
Distance in z−axis from transmitter to receiver (in m)
Distance in z−axis from transmitter to receiver (in m)
3
2.5
2
1.5
1
−5
0
10 10 Amplitude (A/m)
5
10
0.5 −180
−90
0 Phase (degrees)
90
180
Fig. A.3. Solution (magnetic field) along a vertical axis passing through a toroid radiating in a homogeneous medium in the presence of a metallic mandrel. Analytical solution (mandrel is a PEC) against the numerical solution for different mandrel resistivities (10 −7 , 10−5 , 10−3 , and 1 Ω · m) obtained with the self-adaptive goal oriented hp-FEM.
[10] J. R. Lovell, Finite Element Methods in Resistivity Logging, PhD thesis, Delft University of Technology, 1993. [11] J.T. Oden and S. Prudhomme, Goal-oriented error estimation and adaptivity for the finite element method., Comput. Math. Appl., 41 (2001), pp. 735–756. [12] M. Paraschivoiu and A. T. Patera, A hierarchical duality approach to bounds for the outputs of partial differential equations., Comput. Methods Appl. Mech. Eng., 158 (1998), pp. 389– 407. [13] D. Pardo, Integration of hp-adaptivity with a two grid solver: applications to electromagnetics., PhD thesis, The University of Texas at Austin, April 2004. [14] D. Pardo, L. Demkowicz, and C. Torres-Verdin, A Goal Oriented hp-Adaptive Finite Element Method with Electromagnetic Applications. Part I: Electrostatics., ICES Report 04-57. Submitted to Int. J. Numer. Methods Eng., 0 (2005), p. 0. [15] M. Paszynski, L. Demkowicz, and D. Pardo, Verification of Goal-Oriented hp-Adaptivity., ICES Report 05-06, (2005). [16] C. R. Paul and S. A. Nasar, Introduction to Electromagnetic Fields, McGraw-Hill, New York, 1982. [17] S. Prudhomme and J.T. Oden, On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors., Comput. Methods Appl. Mech. Eng., 176 (1999), pp. 313–331. [18] W. Rachowicz, D. Pardo, and L. Demkowicz, Fully automatic hp-adaptivity in three dimensions., Tech. Report 04-22, ICES Report, 2004. [19] R. Rannacher and F.T. Suttmeier, A posteriori error control in finite element methods via duality techniques: application to perfect plasticity., Comput. Mech., 21 (1998), pp. 123– 133. [20] P. Solin and L. Demkowicz, Goal-oriented hp-adaptivity for elliptic problems., Comput. Methods Appl. Mech. Eng., 193 (2004), pp. 449–468. [21] J. Van Bladel, Singular Electromagnetic Fields and Sources., Oxford University Press., New York, 1991. [22] T. Wang and J. Signorelli, Finite-difference Modeling of Electromagnetic Tool Response for Logging While Drilling., Geophysics, 69 (2004), pp. 152–160.
