A Numerical Mode-Matching Method Based on Multiple-Knot Cubic B-Splines Applied to the Analysis of Well-Logging Tools

The numerical mode-matching (NMM) technique is a very efficient approach to analyze well-logging tools used for hydrocarbon exploration. This problem is typically modeled as a cylindrically layered medium (subterranean formation) including both horizontal and vertical stratifications. In this paper, we present a new NMM formulation based on the use of B-splines expansion functions with variable knot multiplicity to represent the dependency of fields on the vertical direction, and cylindrical Bessel functions to represent the fields along the radial (horizontal) direction. We show that this type of expansion allows to better capture the spatial variation of the electromagnetic field near vertical layer interfaces. Illustrative examples elucidate the advantage of knot multiplicity in the modeling of earth formations with anisotropic and lossy layers. A case study on compensated dual resistivity (CDR) well-logging tool sensors is also presented based on the proposed method.


I. INTRODUCTION
OGGING-WHILE-DRILLING (LWD) tools are electromagnetic sensors routinely used to estimate constitutive parameters (permittivity and conductivity) of subterranean formations [1].This information is useful is assessing the potential for hydrocarbon bearing zones.
Continual advances in computing power have enabled the use of brute-force numerical techniques such as finite-differences (FD) [2][3][4][5] and finite volumes (FV) [7] to model LWD tools operating in complex geophysical environments.However, it is still impractical to employ brute-force methods in inversion efforts where a large number of forwards simulation is required for a parameter space sweep.On the other hand, methods that combine both numerical and analytical steps [8][9][10][11][12][13], such as the numerical mode-matching (NMM) [14][15][16][17], are able to model complex geophysical formations while requiring only a modest amount of computational resources compared to brute-force methods.
In this paper, we present a new NMM formulation based on the use of cubic B-splines with variable knot multiplicity to expand the fields along the vertical direction.The knot multiplicity of these basis functions allow a better representation of the field behavior across interfaces while using a coarser computational grid, leading to a lesser demand on computational resources.
The remaining of this work is organized as follows.Section II summarizes the proposed formulation of NMM using cubic Bsplines.Section III shows numerical results to validate the technique, drawing attention to some advantages of the proposed modal expansion based on multiple knots.Section IV presents a case study on compensated dual resistivity (CDR) well-logging tools.Finally, Section V provides some concluding remarks.

II. NUMERICAL MODE-MATCHING FORMULATION
The geometry of this problem can be adequately described in terms of a cylindrically layered medium composed of both radial (horizontal) and vertical stratifications [17].In the NMM formulation presented in [15], the electromagnetic field is expanded in terms of basis functions that depend on the vertical coordinate z, while reflection and transmission matrices (obtained by enforcing the appropriate boundary conditions across interfaces) are used to represent the dependency along the radial coordinate ρ in terms of Bessel and Hankel functions.

A. Radial NMM
Assuming a time-harmonic dependence for the fields and sources in the form exp( ) it   , we obtain the two vector wave equations for the electric and magnetic fields, respectively, as [6]  

12
, and where the (anisotropic) constitutive properties of the medium are represented by the permittivity and permeability tensors,  and  [12,15].The impressed electrical current J is produce by horizontal coil antennas wrapped around a metallic mandrel with radius T  .In addition, we assume the antennal to be axially located at T z .Using a cylindrical coordinate system ( , , ) z  aligned with the tool axis, the fields given in ( 1) and ( 2) can be represented by means of the Bessel transforms [18] as 1 ( , , ) ( , ), 2 where () The substitution of (6) in (5) (7) where in practice, the series are truncated using a sufficient number of terms for a given grid resolution.The resulting problem can be arranged in a matrix form as Equation ( 8) represents a generalized eigenvalue problem where 2 k  are the eigenvalues and a are the eigenvectors.Note that that matrices L  and p  are symmetric in reciprocal media (i.e., with symmetric  and  ).
In regions with sources, the z-component of the fields obey the following equation: where the electric flux density is given by the expansion:  .The expansion described in ( 10) and ( 11) can be used to rearrange (9): by taking the product of both sides of (9) with by   1 pv z    and then applying an inner prod- uct along z, we find where Tq j  represents the spectral contribution of the source for the TE modes.Note that only TE modes are excited by hori- zontal coil antennas as considered here.

B. Transimpedance Calculation
With the expressions of the fields in the source region, we can now obtain the values of the induced voltages at a given receiver coil due to a current in the transmitter antenna.The transimpedance is defined as Where, since the receiver coil antennas are horizontal, we only need to consider the azimuthal electric field component.The latter can be written in the following form [6]: where is a factor depending on the reflection coefficients across layers.See [15] for more details on the evaluation of (13).

C. Cubic B-Spline Eigenfunctions
In the NMM formulations described in [6] and [19], the electromagnetic fields were approximated by triangular, sinusoidal and quadratic B-splines basis functions.In the present work, we adopt cubic B-splines as base (expansion) functions [20] because of their high-order nature and ability to handle knots of higher multiplicity.The function   n Sz employed in (6) can be each written as a combination of four polynomial terms:       Each B-spline has five knots (nodes).When these knots are equally spaced, we obtain a uniform B-spline, as shown in in Fig. 1.When two or more knots are collocated in space, we obtain a knot with higher multiplicity.Fig. 2 illustrates the behavior of a cubic B-spline function for different knot multiplicities.The top curve shows a case with single knots only.The mid curve shows a case with one knot with double multiplicity.The bottom curve shows a case with one knot with triple multiplicity.By increasing the multiplicity, one decreases the regularity of the function.At the knot with double multiplicity the cubic B-spline has discontinuous second derivative, whereas at the knot with triple multiplicity the cubic B-spline has discontinuous first derivative.This properties allows for properly capturing the behavior of the electromagnetic fields from the boundary conditions across layer interfaces.

III. THE KNOTS-MULTIPLICITY EFFECTS
In order to first examine the effect of the multiplicity of knots of the cubic B-splines functions in a simpler setting, we consider the problem of a parallel plate guide filled by two layers in the z-direction, as shown in Fig. 3.This problem admits an analytical solution as described in [21].The eigenvalues for TM and TE (to z ) waves satisfy 2 2 2 0 , and respectively.In the examples that follow, we compare the eigenvalues obtained by the NMM using the cubic B-splines described   .The relative error in the solution of the radial wavenumber k  for the first 9 modes with expansions using multiple knots (MK) or only single knots (SK) are depicted in Fig. 4 and Fig. 5, for the TM and TE fields, respectively.We can observe that lower order modes present a significant reduction by using MK.The differences between the MK and SK results are more pronounced for TM modes.This behavior can be explained by the presence of a double derivative with respect to z in the first term within the brackets in (7).In this case, the MK expansion provides, of course, a more accurate result compared to the SK expansion.

B. Isotropic and Lossy Layers
As a second example, consider now a parallel plate waveguide filled with lossy layers.Layer 1 has 10 .Fig. 6 and Fig. 7 present the relative error in the solution of k  for the first 10 modes with MK and SK.We again observe smaller errors by using MK compared to SK, for both TM and TE modes.   .Fig. 8 and Fig. 9 show results for the first ten modes using SK and MK, for both TM and TE modes, respectively.Once more, the advantage of the MK approach is clear compared to SK.

IV. CASE STUDY: BOREHOLE WITH RUGOSITY
Next, we consider a compensated dual-resistivity (CDR) LWD tool traversing a rugose borehole.This problem was considered before in [22,23] and is examined here to validate our formulation and illustrate the ability of the proposed algorithm based on MK cubic B-splines to efficiently model this kind of problem.The CDR sensor operates at 2MHz and is composed by two transmitters and two receivers as depicted in Fig. 10.Typically, this sensor yields two apparent resistivities, ps R and da R , asso- ciated with the phase difference and the amplitude ratio of the voltage inducted at the receiver antennas.The ps R resistivity pro- vides information of radial regions close to the antennas while the da R resistivity provides information from regions farther from the antennas.
The transmitters and receivers are modelled as coil antennas with radius of 3.5 in.The conductor mandrel is modeled as a PEC cylinder with radius of 3.25 in.The sensor antennas are placed concentrically around the mandrel, and their axial positions are represented in Fig. 10.The CDR sensor moves along a 2.5-in-rugose borehole wall, as depicted in Fig. 11, inside a 3.25-inborehole with resistivity 0.1 Ω.m, surrounded by a soil formation with resistivity 2 Ω.m.For further details on this geometry, see 671 [22,23].
To eliminate spurious reflections arising from the truncation of the domain along z, we employed a cylindrical perfectly matched layer with properties similar to that used in [24].Initially, we use a uniform grid along z with 0.020 m z  within the domain of interest and based on SK.The resistivity results due to transmitter 2 excitation are shown in Fig. 12 (note that, in this uncompensated case, transmitter 1 is not active).The abscissa axis represents the midpoint between the two receivers, i.e.,

 
12 /2 rr zz  . We can clear see that the rugosity in the borehole wall disturbs the response of the sensor resistivity compared to the smooth wall solution.This can lead to misinterpretations when reading the LWD sensor response.In addition, Fig. 13 shows the compensated CDR response obtained via the proposed NMM with MK cubic B-splines.The compensated response is obtained by using both transmitters individually, and subsequently averaging the phase and amplitude results.The discretization grid along z in this case has z=0.03048m in the washout region, and z=0.0217m out of it.Very good agreement is observed versus the results reported in [22].It should be noted that the MK approach requires 55% less discretization points near the vertical interfaces compared to the (SK based) uniform grid.

( 11 )
Equation(11) can represent the electric flux density or the magnetic flux density.For TE modes, we have the  -dependence for each modal field, and similarly   qq zz   

4 Fig. 1 .
Fig.1.The four polynomials that compose a uniform cubic B-spline with five equally space knots.

Fig. 4 .
Fig. 4. Relative error on TM eigenvalues, in the case of isotropic and lossless layers.

Fig. 5 .
Fig. 5. Relative error on TE eigenvalues, in the case of isotropic and lossless layers.

Fig. 6 .
Fig. 6.Relative error on TM eigenvalues, in the case of isotropic and lossy layers.

Fig. 7 .
Fig.7.Relative error on TE eigenvalues, in the case of isotropic and lossy layers.

670C.
Anisotropic and Lossy LayersIn a third example, consider uniaxially anisotropic and lossy layers.The upper layer has 1

Fig. 8 .
Fig. 8. Relative error on TM eigenvalues, in the case of anisotropic and lossy layers.

Fig. 9 .
Fig. 9. Relative error on TE eigenvalues, in the case of anisotropic and lossy layers.

Fig. 12 .
Fig.12.Amplitude and phase apparent resistivities of a CDR sensor tool due the transmitter 2 alone (uncompensated).
is a Bessel function of order v , and( , ) n Sz :