Characterization of borehole fractures by the body and interface waves

(1)

Characterization of borehole fractures

by the body and interface waves

(2)

(3)

Characterization of borehole fractures

by the body and interface waves

Proefschrift

ter verkrijging van de graad van doctor

aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof.dr.ir. J.T. Fokkema,

voorzitter van het College voor Promoties,

in het openbaar te verdedigen op woensdag 26 januari 2005 om 10.30 uur

door

Florence HENRY

Ingénieur en Géophysique et Géotechniques

de l’Institut de Sciences et Technologie (IST)

de l’Université Pierre et Marie Curie (Paris 6)

(4)

Dit proefschrift is goedgekeurd door de promotor Prof.dr.ir. J.T. Fokkema

Toegevoegd promotor: Dr. C.J. de Pater

Samenstelling promotiecommissie:

Rector Magnificus, voorzitter

Prof.dr.ir. J.T. Fokkema Technische Universiteit Delft Dr. C.J. de Pater Technische Universiteit Delft

Prof.dr. F.H. Cornet Institut de Physique du Globe, Parijs Prof.dr. P.K. Currie Technische Universiteit Delft Prof.dr.ir. C.P.A. Wapenaar Technische Universiteit Delft Prof.dr. S.M. Luthi Technische Universiteit Delft Dr.ir. M.D. Verweij Technische Universiteit Delft

Martin Verweij heeft als begeleider in belangrijke mate aan de totstandkoming van het proefschrift bijgedragen.

ISBN: 90-8559-027-2

(5)

Sur l’avenir, tout le monde se trompe. L’homme ne peut être sûr que du moment présent. Mais est-ce bien vrai? Peut-il vraiment le connaître le présent? Est-il capable de le juger? Bien sur que non. Car comment celui qui ne connaît pas l’avenir pourrait-il comprendre le sens du présent? Si nous ne savons pas vers quel avenir le présent nous mène, comment pourrions nous dire que ce présent est bon ou mauvais, qu’il mérite notre adhésion, notre méfiance ou notre haine?

Milan Kundera, L’ignorance

On the future, anybody makes a mistake. Human can only be sure of the present. Nevertheless, can we really know the present? Are we able to judge it? Of course, not. Because how could the person who does not know the future, understand the meaning of the present? If we do not know through which future the present brings us, how could we say that this present is good or bad, and deserves our support, suspicion or hatred?

(6)

Financial Support

The research reported in this thesis has been financially supported by the Dutch Technoloy Foundation (STW) and the Delfrac Consortium.

(7)

Table of Contents

( )

f : signal in quadrature with ˆx f

( )

141

B.1.2 Analytic signal z(t) associated to x(t) 141

B.1.3 Properties for causal functions 141 B.2 Application: signal envelope 142

References 143 Acknowledgements ...145

(9)

Summary

The success of the fracturing process in the oil and gas industry depends on our ability to define the geometry of the hydraulic fracture. To have a method of measurement, which allows to characterize completely a fracture (i.e. width, length and height for an axial fracture) from a single well in reliable way, will be a primordial importance, firstly in term of economical cost and secondly, in term of a better oil recovery. That is why this research has focused on the optimisation of the sonic log measurements in terms of frequency, source-receivers configuration and on the improvement and development of modeling to interpret field data issued from the interaction of the Stoneley and body waves with borehole fractures.

The Stoneley wave measurement is currently a common tool used for detecting the presence of fractures and for evaluating the thickness of a single fracture or the permeability of a close-spaced fracture zone. The fracture characterization from the Stoneley wave measurements could be firstly improved by defining an optimum frequency range; this frequency range must allow to keep a Stoneley wave reflection highly sensitive to the fracture thickness and to have a sufficient wave interaction with the fracture tip. Following the predictions of the existing infinite plane fracture model, we have suggested to use the Stoneley wave measurement in the frequency range of 20 – 50 Hz for the characterization of a transverse fracture and 20 – 100 Hz for an axial fracture for in-situ conditions. In the perspective of optimizing the interaction of the Stoneley wave with the fracture tip and the detection of the reflected tip wave, the propagation of the Stoneley waves in a fluid-filled borehole intersected by a fracture with a finite length has been investigated theoretically and experimentally. The analytical model developed for the finite transverse fracture case has taken into account the formation elasticity and the outgoing and incoming wave motion inside the fracture. The calculated transmission and reflection coefficients have been influenced by the resonant frequency related to the finite fracture length and by the formation elasticity. In general, the experimental results obtained from our laboratory scaled model have shown that the finite dimensions of an axial or transverse near-wellbore fracture have affected the Stoneley wave attenuation. For the case of a transverse fracture, it has been established that for a ratio of fracture length to Stoneley wavelength greater than 4-5, the fracture can be considered to be of infinite extent. The existence of a reflected tip wave has been experimentally emphasized; its time arrival and its amplitude could respectively provide information on the fracture length and the fracture thickness. The position of the receiver in the well has had an importance on the possibility to extract additional information on the fracture from the Stoneley wave; a receiver located at the fracture opening has optimized the detection of the reflected tip wave and the recording of the resonant frequency of the fracture.

To characterize fracture larger than 4-5 Stoneley wavelengths from the wellbore, the interaction of the elastic waves with fractures can be used as an alternative method. The in-situ sonic measurement have shown a potential to image the extension of the fracture at a distance from the single well by the detection and the migration of the mode converted sonic waves at the fracture. In this thesis, a semi-analytical model has been

(10)

Summary

developed to include fracture properties in the sonic imaging processes and thus, to investigate the effects of changing fracture properties on the amplitudes and phases of the diffracted events in order to improve the interpretation of field data. This problem has been treated as an application of direct scattering approach by using the integral representation. By associating the fracture to a distribution of equivalent surface source density, the stress and displacement response of the fracture having finite dimensions is approximated by the linear slip model. The characteristics of the scattered wavefield computed from our semi-analytical model have been sensitive to fracture size, fracture profile and fracture compliances. The modeling results have shown that wave diffractions are generated at the fracture tip and present a significant difference of amplitude according that the fracture compliance falls abruptly to zero or tapers off at the fracture tip. Nevertheless, the question of the validation of the accuracy of the linear slip condition and the determination of the slip variation near the fracture tip is still unanswered. Finally, the angle dependency of the results may offer the possibility to model the diffraction of the elastic wave with an inclined circular plane fracture.

(11)

Samenvatting

Het succes van de stimulatie van olie- en gasputten met hydraulische scheuren hangt af van ons vermogen om de geometrie van zo’n scheur te bepalen. Een meetmethode die in staat is om vanuit het boorgat de wijdte, lengte en hoogte van de scheur te meten, is van fundamenteel belang voor een verbeterde oliewinning. Daarom heeft dit onderzoek zich geconcentreerd op de optimalisatie van akoestische logs wat betreft de frequentie, positie van bron en ontvangers en de verbetering en ontwikkeling van modellen voor de interpretatie van signalen die ontstaan door de interactie tussen de Stoneley-golven en de golven in het gesteente met een scheur bij een boorgat.

De meting van Stoneley-golven is tegenwoordig een algemene methode voor het ontdekken van de aanwezigheid van scheuren en voor de evaluatie van de wijdte van een scheur of de permeabiliteit van een breuk-zone. De karakterisering van een scheur met metingen van de Stoneleygolf zou ten eerste kunnen worden verbeterd door een optimaal frequentie-spectrum te kiezen. Dit frequentie-spectrum moet het mogelijk maken om zowel de transmissie van de Stoneley golf waar te nemen zodat we de scheurwijdte kunnen meten, en om voldoende interactie van de golf te verkrijgen met de scheurtip. Op grond van voorspellingen met een bestaand model voor een oneindige vlakke scheur, hebben wij voorgesteld om onder praktijkcondities de Stoneley golfmeting in het frequentie-bereik van 20 - 50 Herz te gebruiken voor de karakterisering van een transversale scheur en 20 - 100 Herz voor een axiale scheur. We hebben zowel theoretisch als experimenteel de voortplanting van de Stoneley golven onderzocht in een vloeistof-gevuld boorgat dat door een scheur met een eindige lengte wordt doorsneden. Het doel was om een optimale interactie van de Stoneley golf met de scheurtip te verkrijgen en om de detectie van de gereflecteerde golf te optimaliseren.

Het analytische model dat ontwikkeld is voor het geval van een eindige transversale scheur heeft rekening gehouden met de formatiestijfheid en de uitgaande en inkomende golf die langs de scheur loopt. De berekende transmissie- en reflectie-coëfficiënten worden beïnvloed door de resonantie-frequentie in verhouding tot de scheurlengte en de stijfheid van de formatie. De experimentele resultaten die werden verkregen met ons schaalmodel hebben aangetoond dat de eindige afmeting van een axiale of transversale scheur de verzwakking van de Stoneley-golf beïnvloedt. Voor het geval van een transversale scheur, is er vastgesteld dat de scheur als oneindig beschouwd kan worden als die groter is dan 4 tot 5 maal de Stoneley-golflengte. Het bestaan van een gereflecteerde golf van de tip is experimenteel bevestigd; haar aankomsttijd en amplitude konden respectievelijk informatie over de scheurlengte en de scheurwijdte geven. De positie van de ontvanger in het boorgat is van belang om extra informatie over de scheur te krijgen met behulp van de Stoneley-golf; een ontvanger die dichtbij de scheur wordt geplaatst heeft de beste kans om de gereflecteerde golf van de scheurtip waar te nemen.

Een alternatief om de dimensies van scheuren groter dan 4 tot 5 Stoneley-golflengtes vanuit het boorgat te bepalen, is de interactie tussen elastische golven en scheuren te gebruiken. In de praktijk hebben akoestische metingen aangetoond dat ze de grootte

(12)

Samenvetting

van scheuren in beeld kunnen brengen op enige afstand van het boorgat, door het detecteren en migreren van de gereflecteerde golven van de scheur. In dit proefschrift, is een semi-analytisch model ontwikkeld om scheureigenschappen toe te voegen aan de akoestische interacties en zodoende de gevolgen van varierende scheur-eigenschappen op de amplitude en fase van gediffracteeerde golven te onderzoeken, met het doel om de interpretatie van praktijkgegevens te verbeteren. Dit probleem is behandeld als toepassing van de directe verstrooiingsbenadering door het proces te beschrijven met integraalvergelijkingen. Door de scheur te associëren met een distributie van bronnen op het scheuroppervlak, kan de spanning en de verplaatsings-responsie van een eindige scheur worden benaderd met het lineaire slip model. De kenmerken van het verstrooide golfveld dat we berekenen met dit semi-analytische model zijn gevoelig voor scheurgrootte, scheurprofiel en scheur-stijfheid. De modelleringsresultaten hebben aangetoond dat diffracties die worden geproduceerd door de scheurtip significant verschillen als de breuk-stijfheid abrupt naar nul daalt dan wel geleidelijk afneemt bij de scheurtip. Niettemin, blijft het een probleem om de nauwkeurigheid van de lineaire slip voorwaarde te bevestigen en het is ook een open vraag hoe de slipvariatie dichtbij de scheurtip verloopt. Tenslotte biedt de hoekafhankelijkheid van de resultaten misschien de mogelijkheid om de diffractie van de elastische golf met een schuin breukvlak te modelleren.

(13)

Nomenclature and convention

Nomenclature

Units: SI (m = metre, s = second, Pa = Pascal, Hz=Hertz) Dimensions: m = mass, L = length, t = time

Variable Description Units t : time coordinate [s]

f : frequency [Hz]

w : angular frequency [Hz] E : Young's modulus [Pa]

n : Poisson ratio [-] l : Lamé coefficient [Pa] m : Lamé coefficient [Pa] k : compressibility of fluid [Pa-1]

r : volume density of mass [kg.m-3] VP : P-wave velocity [m.s-1]

VS : S-wave velocity [m.s-1]

Vf : fluid wave velocity [m.s-1]

k : angular wavenumber [m-1] p : acoustic pressure [Pa]

ui : particle displacement [m]

vi : particle velocity [m.s-1]

tij : stress in elastic medium [Pa] eij : strain in elastic medium [-]

fi : volume source density of force [Pa.m-1]

hij : volume source density of strain [-]

q : source density of volume injection rate [s-1]

Cijpq : stiffness tensor of an elastic solid [Pa.m-1]

Smnij : compliance tensor of an elastic solid [m.Pa-1]

a : radius of the wellbore [m] d : fracture width [m] L : fracture length, measured from borehole wall [m] H : fracture height [m] c : phase velocity of the interface wave [m.s-1]

(14)

Nomenclature and convention

Convention for the temporal and spatial Fourier Transforms

Temporal Fourier transforms:

(

)

( )

0 ˆ , , exp( ) y z w =¥

_ò

y z t i t dtw

( )

, 1 ˆ

(

,

)

exp( ) 2 y z t y z w i t dw w p +¥ -¥ =

_ò

-Spatial Fourier transforms in z-direction:

(

_z,

)

ˆ

(

,

)

exp( _z ) y k w +¥y z w ik z dz -¥ =

_ò

-%

(

)

1

(

) ( )

ˆ , , exp 2 z z z y z w y k w ik z dk p +¥ -¥ =

_ò

%

Conversion Factors

1 inch = 0.0254 m = 25.4 mm 1 ft = 0.3048 m 1 cP = 0.001 Pa s 1 lbf s ft -2 _{= 47.88 Pa}_s 1 lb = 453.59 g 1 gallon = 3.7853 l 1 barrel = 159.98 l =0.15998 m3 1 b/min = 0.00267 m 3_s -1 1 m = 39.3701 inch = 3.2808 ft 1 Pa s = 1000 cP = 0.0209 lbf s ft -2 1 kg = 2.3051 lb 1 l = 0.2642 gallon 1 m3₌_6.2508_bl

(15)

CHAPTER 1 Introduction

1.1 General introduction and purposes of the thesis

In low permeability reservoirs, production is often enhanced by fracturing. For lack of a presence of natural fractures in reservoir layer, hydraulic fracturing has become common practice in oil and gas production as a technique to increase the recovery of hydrocarbons. The operation consists of injecting a fluid in a wellbore under high pressure to generate hydraulic fractures in the reservoir layer in order to allow hydrocarbons to be drained. The geometry of the induced fractures is critically important for their efficiency in production. The probability of success of the fracturing process depends on our ability to define the shape and the size of the fracture. The detection and characterization of permeable fractures allow controlling fluid movement through the fractures. The two main fracture parameters which are important for the production are the fracture permeability and the fracture dimensions. Efficient methods for detecting and characterizing natural or hydraulic fractures are still a subject of investigation. The detection of the presence of the fracture can be done from a variety of methods such as borehole televiewer (Zemanek et al., 1969) and electrical borehole scans ( Luthi and Souhaite, 1990). Besides the fracture detection, the permeability and the extension of the fracture are essential to define the resource production. The passive monitoring, based on the mapping of induced seismicity generated during hydraulic fracturing operation, represents an efficient way of mapping hydraulic fracture extensions (Boadu, 1997; Brady et al., 1994; Jupe et al., 1998; Phillips et al., 1996; Sarda et al., 1988). However, this method cannot currently provide information on the fracture thickness or permeability. Another practical approach to investigate the characteristics of the fracture is by means of active monitoring methods.

(16)

CHAPTER 1 - Introduction

These methods must use frequencies that are high enough to provide spatial resolution at the scale of individual fractures, and low enough to penetrate into the formation. Four main borehole measurements are currently used to characterize borehole fractures: 1) Hydraulic Impedance Testing (HIT) (Holzhausen and Gooch, 1985; Paige et al., 1995), used at a frequency lower than 10 Hz, can theoretically characterize completely the fracture; however, this method appears to be relatively insensitive to large field fracture dimensions. 2) Vertical Seismic Profiling (VSP) (Beydoun et al., 1985; Cicerone and Toksöz, 1990; Toksöz et al., 1992; Meadows and Winterstein, 1994) acts in a higher frequency range: 50 – 100 Hz. VSP measurements provide an estimate of the fracture length; however the technique suffers from ambiguous interpretation of fracture size and requires difficult and expensive acquisition geometry. 3) Cross-well seismic ( Majer et al., 1997) represents an alternative to VSP. The measurements can be done in a higher frequency range: 50 – 250 Hz but they present the same limitations as VSP measurements. 4) Sonic Log (Brie et al., 1998; Gelinsky and Cheng, 1998) uses a source and a receiver array located in the borehole and works in a frequency range: 1 – 80 kHz. The sonic log is known to be mainly sensitive to fracture width or permeability.

Although these borehole methods have shown potential to measure fracture parameters, they generally exhibit uncertainty in the interpretation of the data. Among the methods which provide information of the fracture, the sonic log has the advantage that it can combine the generation of surface waves along the borehole wall like the Stoneley wave and of elastic waves in the formation and that it has a better resolution than seismic waves. The Stoneley wave is sensitive to open fractures that can contribute most to the production (Paillet and White, 1982; Hsu et al., 1985). In the literature, several authors ( Hornby et al., 1989; Tang, 1990; Tang and Cheng, 1993; Kostek et al., 1998) have analytically treated the interaction of the Stoneley wave with an infinite plane fracture and have provided a quantitative association between the amplitude attenuation of the Stoneley wave and fracture width by means of its reflection and transmission coefficients. An estimate of the fracture permeability can be defined from the fracture width. On the other hand, the effect of a fracture with a finite length on borehole wave propagation has not attracted so much attention. Although the Stoneley wave is a guided wave, it can be assumed that for a certain ratio of the Stoneley wavelength to fracture radius, the fracture tip by acting as a secondary source of wave, can modify the Stoneley wave attenuation, especially for near-wellbore fractures. Moreover, when the Stoneley wave passes the fracture, a fracture wave is generated and propagates radially away from the wellbore at approximately the speed of sound in the fluid. If this wave was reflected at the fracture tip and was observable in the borehole, it would provide a measure of the fracture size. The detection of the reflected tip wave can have some applications in the near-wellbore fracture characterization since a significant attenuation of the fracture wave can be expected. By investigating theoretically and experimentally the case of a fracture with a finite radial extent, our study aims to quantify the effect of such fracture geometry on the borehole wave propagation.

Besides a sensitivity of the Stoneley wave attenuation to the fracture properties, a frequency dependence of the reflection and transmission coefficients of the Stoneley wave was established by dynamic wave excitation (Hornby et al., 1989; Tang, 1990). The reflection coefficient at an open infinite fracture is close to one at low frequency

(17)

and decreases for higher frequencies and, inversely, for the transmission coefficient. Therefore, an improvement of the borehole measurement can be expected by defining an optimum frequency range for which reflection and transmission coefficients are more sensitive to the fracture characteristics under field conditions and for which a reasonable amount of energy is transmitted into the fracture.

Since fracture permeability can be estimated by the Stoneley wave log method, it is important to see the extension and the profile of such fractures at a distance from the borehole to assure fluid production through the fractures. The sonic measurements have shown a potential to characterize the fracture at a distance from the well. The monitoring of the fracture extent from the borehole has been investigated by many authors (Yamamoto et al., 1998). The method needs a sonic imaging technique whose resolution is better than the seismic one in order to allow a detection of thin fractures. Many authors (Yamamoto et al., 1998; Gelinsky and Cheng, 1998) managed to image the extension of the fractures by the migration of the mode converted sonic waves. However, the response of the fracture to the sonic wave by studying the amplitudes and phases of the events was still not well exploited to have access to the physical properties of the fracture. In the previous phases of the project, by using the ultrasonic imaging lab data, a theoretical interpretation of diffraction and transmission measurements of the elastic waves to determine the extent and width of fractures has been developed (Savic, 1995; Groenenboom, 1998). However, the application of this method was limited to interpret cross-hole measurements in the field. Extending this result for sonic measurements done from a borehole allows improving the interpretation of the field data. Such a problem will be treated by the direct scattering approach for a thin transverse fracture in order to investigate in which proportion the interaction of the body waves with the fracture and especially, their diffraction at the fracture tip can have some effects on the borehole signal and can represent an additional source of information on the fracture characteristics for the sonic log configuration.

1.2 Outline of the thesis

This research project focuses on the optimization of the sonic log measurements by using an appropriate frequency range and on the improvement of the interpretation of field data by quantifying the effect of a finite fracture on the borehole wave propagation; the existence of a reflected tip wave in waveform is also investigated. This thesis is divided into six chapters. This chapter has defined the subject of the thesis. Chapter 2 concerns the determination of the acoustic field due to a source in a fluid-filled borehole surrounded by an elastic formation and will complement chapter 4 about the direct scattering problem approach. Chapter 3 gives an overview of the existing analytical models of the interaction of the Stoneley wave with a borehole fracture. Then, a study about the sensitivity of the Stoneley wave to the borehole fracture parameters as function of the frequency measurement is realized for optimizing the frequency of the borehole measurements. This sensitivity analysis is also used as basis for a dimensional study in order to be able to make an extrapolation of the laboratory results to the in-situ conditions. Finally, an analytical model is developed to investigate the interaction of the Stoneley wave with a finite transverse fracture. This

(18)

model shows that the transmission and reflection coefficients of the Stoneley wave are sensitive to the resonance frequency of the fracture and to the fracture radial extent at very low frequencies. The fact to take into account the elasticity of the formation affects the resonant frequency of the fracture and the value of the transmission and reflection coefficient at low frequency. In chapter 4, the interaction of the body waves with a thin transverse borehole fracture is investigated by the direct scattering problem approach. The formulation of the scattered wavefield is obtained by introducing the linear slip theory to model the fracture. Chapter 5 concerns the experimental studies on the wave propagation in a borehole intersected by a finite transverse fracture and a finite axial fracture, respectively. Chapter 6 summarizes the important results and major conclusions of this thesis.

References

Beydoun W.B., Cheng C.H. and Toksöz M.N. (1985), “Detection of open fractures with vertical seismic profiling”, J. Geoph. Res., vol.90, p.4557-4566.

Boadu F.K. (1997), “Relating the hydraulic properties of a fractured rock mass to seismic attributes: theory and numerical experiments”, Int. J. Rock Mech. Sci., vol.34, p885-895.

Brady J.L., Withers R.J., Fairbanks T.D. and Don Dressen (1994), “Microseismic monitoring of hydraulic fractures in Prudhoe Bay”, SPE28553.

Brie A., Endo T., Ito H. And Badri M. (1998), “Fracture and permeability in a fault zone from Sonic waveform data”, proceedings 4th SEGJ Tokyo, p191-197.

Cicerone R.D. and Toksöz M.N. (1990), “A dynamic model of tube wave generation at fracture in hydrophone VSP data”, Reservoir delineation – Vertical Seismic Profiling Consortium, Annual Report.

Gelinski S. and Cheng C.H. (1998), “Integrated borehole acoustic fracture characterization”, proceedings, 4th SEGJ Tokyo, p111-116.

Groenenboom J. (1998), “Acoustic monitoring of hydraulic fracture growth”, PhD thesis, Delft University of Technology.

Holzhausen G.R. and Gooch R.P (1985), “Impedance of hydraulic fracture: its measurement and use for estimating fracture closure pressure and dimensions”, SPE13892.

Hornby B.E., Johnson D.L., Winkler K.W. and Plumb R.A. (1989), “Fracture evaluation using reflected Stoneley wave arrivals”, Geophysics, vol. 54, p.1274-1288. Hsu K., Brie A. and Plumb R. (1985), “A new method for fracture identification using array sonic tools”, SPE 14397.

(19)

Jupe A., Jones R., Dyer B. and Wilson S. (1998), “Monitoring and management of fractured reservoirs using microearthquake activity”, SPE/ISRM47315.

Kostek S., Johnson D.L., Winkler K.W. and Hornby B.E. (1998), “The interaction of the tube waves with borehole fractures, part 2: Analytical model”, Geophysics, vol.63, p809-815.

Luthi S.M. and Souhaite P. (1990), “Fracture aperture from electrical borehole scans”, Geophysics, vol. 55, p.821-833.

Majer E.L., Peterson J.E., Daley T., Kaelin B., Myer L., Queen J., D’Onfro P. and Rizer P. (1997), “Fracture detection using crosswell and single well surveys”, Geophysics, vol.62, p.495-504.

Meadows M.A. and Winterstein D.F. (1994), “Seismic detection of a hydraulic fracture from shear-wave VSP data at Lost Hills field, California”, Geophysics, vol. 59, p.11-26.

Paige R.W., Roberts J.D.M., Murray L.R. and Mellor D.W. (1992) “Fracture measurement using hydraulic impedance testing”, SPE 24824.

Paige R.W., Murray L.R. and Roberts J.D.M. (1995) “Field application of hydraulic impedance testing for fracture measurement”, SPE Production & Facilities.

Paillet F.L. and White J.E. (1982), “Acoustic modes of propagation in the borehole and their relationship to rock properties”, Geophysics, vol.47, p.1215-1228.

Phillips W., Rutledge J., Fairbanks T. Gardner T., Miller M. and Schuessler B. (1996), “Reservoir fracture mapping using microearthquakes: Austin Chalk, Gidding Field”, SPE36651.

Sarda J-P., Perreau P-J. and Deflandre J-P. (1988), “Acoustic emission interpretation for estimating hydraulic fracture extent”, SPE17723.

Savić M. (1995), “Ultrasonic scattering from a hydraulic fracture: Theory, Computation and experiment”, PhD thesis, Delft University of Technology

Tang X.M. (1990), “Acoustic logging in fractured and porous formations”, PhD thesis, MIT.

Tang X.M. and Cheng C.H. (1993), “Borehole Stoneley wave propagation across permeable structure”, Geoph. Prospecting, vol. 41, p.165-187.

Toksöz M.N., Cheng C.H. and Cicerone R.D. (1992), “Fracture detection and characterization from hydrophone vertical seismic profiling data”, Fault mechanics and transport properties of rocks, Academic Press, chap.16.

Yamamoto H., Watanabe S., Mikada H., Endo T. and Brie A. (1998), “Fracture imaging using borehole acoustic reflection survey”, proceedings, 4th SEGJ Tokyo, p375-382.

(20)

Zemanek J., Caldwell R.L., Glenn E.E., Holcomb S.V., Norton L.J. and Straus A.J.D. (1969), “The borehole televiewer – A new concept for fracture location and other types of borehole inspection”, J. Petr. Tech., vol. 25, p.762-774.

(21)

CHAPTER 2 General formulation of the

source-borehole configuration

In this chapter, the basic equations of the elastic wave theory derived from the conservation laws are firstly given. These basic equations are combined to obtain the wave equations for an ideal acoustic medium and for an elastic medium from which the formal solutions can be derived. By using the Fourier method and the separation of the variables, we have formulated the expression of the acoustic field due to a volume source of deformation rate in a fluid-filled borehole surrounded by an elastic solid formation. The total acoustic field is decomposed into a direct field due to the action of the source alone and a reflected field due to the interaction of the source with the formation. Then, the characteristic equation for a fluid-filled borehole is derived and some properties of the waves propagating in the borehole are also discussed.

2.1 Introduction

In this chapter, we formulate an expression for the acoustic field due to a volume injection source in a fluid-filled borehole surrounded by an elastic formation. First, the form of the acoustic field in the borehole and then the form of the elastic field in the formation are derived by the Fourier method and the separation of the variables as described in the work of Zou (1993).

The acoustic field in the borehole fluid due to a source is a superposition of two components. The first component is a contribution due to the volume injection source alone and given by the field that the source would radiate if the medium was an infinite homogeneous fluid. This solution is referred to the “direct” component of the acoustic

(22)

CHAPTER 2 – General formulation of the source-borehole configuration

field. The second component is due to an interaction of the source with the formation and is referred to as the “reflected field”. The solution in the solid is derived from the source free equations and is obtained after the separation of compressional and shear waves. By matching boundary conditions at the fluid-solid cylindrical interface, the representations of the wavefield for the source-borehole configuration are achieved. Then, some properties of the waves propagating in the borehole are discussed.

The geometry of the borehole is shown in figure (2.1). This consists of a fluid-filled borehole of infinite length surrounded by an infinite elastic medium, with an axially symmetric source centered on the z-axis of the borehole. In the figure, the radius of the borehole is denoted by a. The properties of the borehole rock formation and fluid are specified by subscript s and f, respectively.

The origins of both Cartesian and cylindrical coordinates are placed at the center of the source. As a consequence of the axial symmetry of our problem, the resulting wave field will be described with cylindrical coordinates.

2.2 Basic Equations

In this part, we give the basic equations of the theory of the elastic wave and of the acoustic wave derived from the conservation laws as described in Achenbach (1973). These equations serve as the point of departure for the determination of the wave equation and the formal solutions.

a eX eY e_Z 0 r q

(23)

For an elastic solid, the linearized equation of motion is written as

j ijt r t iv fi

¶ - ¶ = - , (2.1)

where fi is the volume source density of force.

In the equation of deformation, the general form of the deformation rate is composed of an “induced” part, which is denoted by ėij and is related to the stress in the solid, and an

“external” part, which is denoted by hij and is representative for the action of external

sources that apply a deformation rate to the solid. Thus, the equation of deformation for an elastic solid is expressed as

(

)

1

2 ¶ + ¶i jv j iv - =e&ij hij, (2.2)

where hij is also called the source density of strain.

In equation (2.3), the stress tensor tij is by definition linearly related to the strain tensor

epq by Hooke’s law for a solid. In its general form, the stress strain relation is written as

ij C eijpq pq

t = , (2.3)

where Cijpq is the stiffness tensor of the elastic solid. The stiffness tensor C has the

following symmetry properties

ijpq jipq jiqp ijqp

C =C =C =C . (2.4)

For an isotropic and perfectly elastic solid, the stiffness can be expressed by

(

)

ijpq ij pq ip jq iq jp

C =ld d +m d d +d d , (2.5)

in which _{l and m are Lamé coefficients and d}ij is symmetrical unit tensor.

In our further analysis, we also need to introduce the compliance tensor S that is interrelated to the stiffness tensor C by

(

)

1 2

mnij ijpq mp nq mq np

S C = d d +d d . (2.6)

The compliance tensor satisfies the symmetry relations

ijpq jipq jiqp ijqp

(24)

For a perfect fluid, no shear stresses can exist, which means that the Lamé coefficient m is equal to zero. As a consequence, all normal stresses are identical and equal to the negative pressure according to

ij p ij

t = - d . (2.8)

The equation of motion for a fluid is obtained by substituting (2.8) into equation (2.1)

ip r t iv fi

¶ + ¶ = , (2.9) where fi remains the volume source density of volume force, which is representative of

the action of acoustic source of the “dipole” type.

The standard form of the deformation rate equation is written as

i iv q q

¶ - =& , (2.10)

in which q&= - ¶k _tp and represents the induced part of the cubic dilatation rate. The second term q is the volume source density of injection rate, which is representative for the action of acoustic source of the “monopole” type.

In terms of particle displacement, by using the fact that the particle velocity and the particle displacement are interrelated by

i t i

v = ¶u . (2.11)

The equation of motion and the equation of deformation are written as

2 _, , i t i i i i p u f u p h r k ¶ + ¶ = ¶ + = (2.12)

in which k is the compressibility of the fluid and h is the source density of bulk strain that is related to q by

t

q= ¶h. (2.13)

2.3 Formal solution for the wave in a source free fluid-filled

borehole

We consider firstly the case of a source free fluid-filled borehole in order to determine the acoustic field referred to as “reflected” field.

(25)

The acoustic wave equation follows by eliminating vi from the equation of motion (2.9)

and the pressure strain relation (2.10) (Achenbach, 1973), yielding

2 2 2 ₀

tp Vf p

¶ - Ñ = , (2.14)

where V_f =

( )

r k_f -1/ 2. (2.15) In our problem, the azimuthal symmetry is assumed therefore we have an invariance followingq, i.e. ¶q = 0. The fluid pressure pf satisfies the following homogeneous scalar

wave equation

(

)

2 2 1 2 ₀ t p Vf r r pr zp r é ù ¶ - _ê ¶ ¶ + ¶ _ú= ë û . (2.16)

In terms of pf, the fluid displacement is given by:

2 1 1 tur rpf r ¶ = - ¶ , (2.17) . 2 1 1 tuz zpf r ¶ = - ¶ (2.18)

In order to get the formal solutions, both temporal and spatial transformations are applied to equations (2.16), (2.17) and (2.18). The spatial Fourier transformation is performed along the z-axis in which the geometry of the medium properties is infinite and invariant. The transformed equations are given as

2 2 2 f2 ₀ r f r f r f r ¶ p% + ¶r p% -r k p% = , (2.19) 2 1 rf r f f u p r w = ¶ % % , (2.20)

( )

2 1 zf z f f u ik p r w = % % , (2.21)

in which the radial fluid wavenumber is defined as 1/ 2 2 2 2 f r z f k k V w æ ö =ç_ç - ÷_÷ è ø . (2.22)

(26)

Equation (2.19) is recognized as a modified Bessel equation of zeroth order. For a finite fluid domain i.e. r < a, the solution of the equation that remains bounded as r®0, is

(

,

)

0

( )

refl f

f z r

p% =A k w I k r , (2.23)

where kz is the axial wavenumber, k_r( )f is the fluid radial wavenumber, I0 is the modified

Bessel function of the first kind of order 0 and A is a constant to be determined by using the boundary conditions at the borehole wall.

The corresponding expressions for the displacements ur and uz in the fluid are given as

(

)

1

( )

2 , f refl r f rf z r f k u A k w I k r r w = % , (2.24)

(

)

0

( )

2 , refl z f zf z r f ik u A k w I k r r w = % (2.25)

where I1 is the modified Bessel function of the first kind of first order.

Equations (2.23), (2.24) and (2.25) are the formal solutions of the pressure and the displacements in a source free bounded fluid medium.

2.4 Formal solution for the wave in the solid

2.4.1. Axisymmetric solutions for the displacements in the elastic medium

The homogeneous equation of motion in terms of displacement in an elastic medium is expressed as (Achenbach, 1973)

( )

2 2 2 . t VP VS ¶ u = Ñ Ñ - Ñ´Ñ´u u, (2.26) with _P 2 s V l m r + = and _S s V m r = , (2.27)

where VP and VS are the compressional and shear velocities, respectively.

For a problem with a cylindrical geometry, the solution of equation (2.26) is expected as a superposition of two Bessel functions. One of these functions is a function of the shear wavenumber and the other is a function of the compressional wavenumber in the solid (Zou, 1993).

(27)

The particle displacement field can be considered as having two components: a compressional field _{u and a shear field}( )P _{u , defined by:}( )S

( )P ( )S

= +

u u u , (2.28)

where _Ñ´_u( )P ₌0 _and _Ñ_._u( )S ₌₀

.

In assuming that u and ( )P u are linearly independent, we can deduce from equations ( )S (2.26) and (2.28) that ( )

(

( )

)

2 P 2 _. P t VP ¶ u = Ñ Ñu , (2.29) ( ) ( ) 2 S 2 S t VS ¶ u = - Ñ´Ñ´u (2.30) Due to the axial symmetry of our wavefield, the tangential component of the wave field does not exist and the radial and axial components are independent of q. Therefore, we have r z u u = r+ z u e e , (2.31) with ( )P ( )S r r r u =u +u , (2.32) ( )P ( )S z z z u =u +u . (2.33)

From both vectorial equations (2.29) and (2.30), we obtain a pair of scalar equations for the radial and axial components for the compressional and shear wave fields

( )

( ) 2 ( )P 2 1 P P tur VP r r rur r z zu r ì é ù ü ¶ = _í¶ _ê ¶ _ú+ ¶ ¶ _ý ë û î þ, (2.34) ( )

( )

( ) ( ) 2 P 2 1 P 2 P tuz VP z r rur z zu r ì é ù ü ¶ = _í¶ _ê ¶ _ú+ ¶ _ý ë û î þ, (2.35) ( )

(

( ) ( )

)

2 S 2 S S tur VS z z ru r zu ¶ = ¶ ¶ - ¶ , (2.36) ( )

(

( ) ( )

)

2 S 21 S S tuz VS r r z ru r zu r é ù ¶ = - ¶ _ë ¶ - ¶ _û. (2.37)

(28)

Temporal and spatial Fourier transformations are applied to the equations (2.34), (2.35), (2.36) and (2.37) ( )

( )

( ) ( ) 2 P 2 1 P P r P r r r z r z u V ru ik u r w ì é ùü - = _í¶ _ê ¶ + ¶ _ú_ý ë û î þ % % % , (2.38) ( ) ( )2

( )

( ) 1 P z P z P r r r ik u ru r k = ¶ % % , (2.39) ( ) ( )2 ( ) S z S r S r z r ik u u k = - ¶ % % , (2.40) ( )

(

( )

)

2 S 21 S S z S r z r r z u V r ik u u r w % = ¶ é_ë % - ¶ % ù_û, (2.41) with ( ) 1 2 ₂ 2 2 P r z P k k V w æ ö =_ç - _÷ è ø , (2.42) ( ) 1 2 ₂ 2 2 S r z S k k V w æ ö =_ç - _÷ è ø , (2.43) where ( )P r

k is the radial compressional wavenumber and ( )S r

k is the radial shear wavenumber.

Determination of the compressional displacement field components: ( )P r

u% and ( )P z

u%

The formal solutions for the displacements in elastic medium are derived from the coupled equations (2.38) and (2.39). By replacing the expression of ( )P

z

u% in equation (2.38), we obtain the following modified Bessel equation of the first order

( ) ( ) ( )2 ( )

2 2 P P P 2 ₁ P ₀

r r r r r r

r ¶ u% + ¶r u% -é_ëk r + ù_ûu% = . (2.44) In the infinite elastic formation around the borehole, the only solution of this differential modified equation allowed by the radiation condition is a function which is bounded when r ® ¥ ( )

₍

₎

( )

( ) 1 , P P r z r u% =B k w K k r . (2.45)

(29)

The expression for ( )P z

u% can be obtained from equation (2.39):

( ) ( )

(

,

)

0

( )

( ) P z P z P z r r ik u B k K k r k w = -% , (2.46)

where B(kz,w) is a coefficient and Kn(.) is the modified Bessel function of the second

kind of the n-th order.

Determination of the shear displacement field components: ( )S r

u% and ( )S z

u%

We proceed in the same way as previously to determine the shear displacement field components. In substituting ( )S

r

u% in equation (2.41) by its expression given in (2.40), we obtain the following differential equation:

( ) ( )

( )

( ) 2 ( )

2 2 S S S S ₀

r z r z r z

r ¶ u% + ¶r u% - k r u% = . (2.47) The solution of this differential modified Bessel equation has to remain bounded when

r® ¥ , thus ( )

₍

₎

( )

( ) 0 , S S z z r u% =C k w K k r , (2.48)

where C(kz,w) is a coefficient. Therefore, the expression of u%( )rS is derived from (2.40):

( ) ( )

(

,

)

1

( )

( ) S z S r S z r r ik u C k K k r k w = % . (2.49)

The particle displacement u% and _r u%zin the elastic solid are given by

( ) ( )

(

)

( )

( ) ( )

(

)

( )

( ) 1 1 , , P S r r r P _z S z r S z r r u u u ik B k K k r C k K k r k w w = + = + % % % (2.50) ( ) ( ) ( )

(

,

)

0

( )

(

,

)

0

( )

( ) . P S z z z P S z z r z r P r u u u ik B k K k r C k K k r k w w = + = - + % % % (2.51)

2.4.2. Formal solutions for the stresses in an elastic medium

In this section, the components of the stress t ,_rr t ,_qq t ,zz t ,rq t and rz t that are useful in qz

(30)

in the definition of the scattered pressure wave field of chapter 3, are calculated. The relations between the stresses and the displacements in an elastic medium (Achenbach, 1973) are given as

(

2

)

1 rr r ru ur z zu r t = l+ m ¶ +læ_ç + ¶ ö_÷ è ø,

(

2

)

r

(

)

r r z z u u u r qq t = l+ m + ¶l + ¶ ,

(

2

)

r zz r r z z u _u u r t =læ_ç¶ + ö_÷+ l+ m ¶ è ø , (2.52) 0 rq qz t =t = and

(

)

rz z ru r zu t =m ¶ + ¶ .

By applying the temporal and z-spatial Fourier transforms, the equations (2.52) become

(

2

)

1 rr r ru _rur ik uz z t = l+ m ¶ +l æ_ç + ö_÷ è ø % % % % ,

(

2

)

r

(

)

r r z z u u ik u r qq t% = l+ m % + ¶l % + % ,

(

2

)

r zz r r z z u _{ik u} u r t =læ_ç¶ + ö_÷+ l+ m è ø % % % % , (2.53) 0 rq qz t% =t% = , and

(

)

rz ik uz r r zu t% =m % + ¶ % .

The expression of t% ,_rr t% , qq t% ,zz t% ,rq t% and rz t% can be obtained by substituting (2.50) and qz

(31)

(

)

( )( ) ( )

( )

(

)

( )

( ) ( )

( )

( ) 2 2 0 1 0 1 2 , 2 1 , 2 , P P P P z r rr z P r r r r S S z z r S r r k k B k k K k r K k r r k C k ik K k r K k r k r m t w l m w m éæ æ - ö ö ù = êç_ç ç_ç ÷_÷- ÷_÷ - ú ê_è _è _ø _ø ú ë û é ù - _ê + _ú ë û % (2.54) 2 ( ) ( ) ( ) 0 1 ( ) ( ) 1 ( ) 2 ( , ) ( ) ( ) ( ) 2 ( , ) , P P P z r z P r r r S r z z S r k k B k K k r K k r k r K k r ik C k k r qq m t w l m w é æ _- ö ù = _ê _ç _÷ + _ú è ø ë û é ù + _ê _ú ë û % (2.55) 2 ( ) ( ) 0 ( ) ( ) 0 ( 2 ) ( , ) ( ) ( , ) 2 ( ) , P z P r zz z P r r S z z r k k B k K k r k C k ik K k r l l m t w w m éæ ö ù - + + = _ê_ç _÷ _ú è ø ë û + é_ë ù_û % (2.56) 0 rq qz t% =t% = , (2.57) and

(

)

( )

₍

₎

( )

( ) ( ) ( ) 2 2 1 1 2 , , S P S z r rz z z r z r S r k k ik B k K k r C k K k r k t =mêé w - w _çæç + ÷ö_÷úù ê è øú ë û % . (2.58)

2.5 Source specification and direct acoustic field

In this section, we define the solution that is referred to as the “direct” component of the acoustic field. It corresponds to the field that originates from the source and propagates in an infinite homogeneous fluid.

We consider a finite source term equivalent to a volume injection source, also called volume source of deformation rate. As a consequence, the term representing the force source fi in the basic equations for the fluid equals to zero. From this assumption, the

inhomogeneous wave equation of pressure is expressed as

2 2 2 1

tp Vf p _k tq

¶ - Ñ = ¶ . (2.59)

(32)

CHAPTER 2 – General formulation of the source-borehole configuration 2 2 2 2 f2 r r r i r r p r p r k p w q k ¶ + ¶ -% % % = - %. (2.60)

Firstly, we define the general form of the solution for the homogeneous wave equation of pressure. The type of this solution is identical to the source free wave equation in an infinite fluid

(

)

( ) ( ) 0 0 , , ( ) ( ) direct f f z v r v r p% r k w =D K k r +E I k r , (2.61) where Dv and Ev are constant and linearly related to the magnitude of the volume source.

The solution has to satisfy the radiation condition therefore it has to remain bounded when r® ¥

(

)

( ) 0 , , ( ) direct f z v r p% r k w =D K k r . (2.62)

By deriving in time (2.20), the form of the incident velocity field can be defined from the pressure field given in (2.62) by

(

)

( )

(

)

(

( )

)

1 , , , f direct r f r z v z r f ik v r k w D k w K k r r w = % . (2.63)

In order to determine the constant Dv, we consider the following volume injection

source

(

)

( )

(

)

ˆ , , ˆ , ,

q x y z =m w d x y z , (2.64)

where _{m w}ˆ_{( )} represents the source spectrum and d

(

x y z, ,

)

is the delta pulse function. In cylindrical coordinates, equation (2.64 ) becomes

(

)

( )1

(

)

ˆ , , ˆ , ,

q r z m r z

r

q = _w d q . (2.65)

A spatial Fourier transform is applied to (2.65) according to the z-axis and gives

(

, , _z

)

ˆ( )1

( )

,

q r k m r

rd q

q = w

% . (2.66)

Next, we consider the basic equation (2.10) expressed in the frequency domain

ˆ ˆ ˆ

. - iwkp q

(33)

which after transformation to cylindrical coordinates and the application of the Fourier transformation with respect to z, becomes

( )

1

r rvr ik vz z i p q

r¶ % + % - wk % %= . (2.68)

Both sides of this equation are now integrated over an infinitesimally small disk with a radius D ®r 0. For the right-hand side of equation (2.68), this gives

( )

( ) 2 2 1 ₀ ₀ 0 0 0 ˆ ˆ lim lim , r r r I qrdrd m r drd m p p q w d q q w D D ® D ® =

_{ò ò}

% =

_ò

= . (2.69)

In order to catch the singularity of the field in the origin, we perform the integration of the left-hand side of equation (2.68) as follows

( )

(

)

( )

2 ₀ 0 1 lim r r r z z r r I rv r ik v r i p r r drd r p wk q D D ® -D ì ü = _í ¶ + - _ý î þ

ò ò

% % % . (2.70)

Since r v%_z

( )

r and r p%

( )

r tend to zero for r ®0, these terms do not contribute to the integral I2 and we find

( ) 2 lim 2_r ₀ r

I p rv r

D ®

= D % D . (2.71)

In view of the fact that I1 must be equal to I2 and by considering the asymptotic behavior

of the Bessel function in the expression of the incident velocity field for

0

r® (Abramovitz and Stegun, 1970)

(

( )

)

( )

1 ( ) ( ) 0 1 2 1 lim 1 2 f r f f r r r K k r k r k r D ® D = G D = D . (2.72)

It follows from equations (2.63), (2.69), (2.71) and (2.72) that the constant Dv is defined

as

(

,

)

ˆ( ) 2 f v z i m D k w r w w p -= . (2.73)

The pressure field and the displacement field due to the volume injection source are given by

(

)

( )

(

( )

)

0 ˆ , , 2 f direct f z r i m p r k w r w w K k r p -= % (2.74) and

(34)

(

)

( ) ( )

(

( )

)

1 2 ˆ 1 , , 2 f direct direct r f r z r r f im k u r k w p w K k r r w pw = ¶ = % % . (2.75)

2.6 Boundary conditions at the solid-fluid interface

In the fluid-filled borehole, the total field is the sum of the “reflected” field and the “direct” field due to the source:

(

, ,

)

(

, ,

)

(

, ,

)

total refl direct

f z f z f z

p% r k w = p% r k w +p% r k w , (2.76)

(

, ,

)

(

, ,

)

(

, ,

)

total refl direct

f r kz w = f r kz w + f r kz w

u% u% u% . (2.77)

In the surrounding formation, the total fields are only composed of the outgoing solutions defined in section 2.4

(

, ,

)

(

, ,

)

total ijs r kz ijs r kz t% w =t% w , (2.78)

(

, ,

)

(

, ,

)

total s r kz w = s r kz w u% u% (2.79)

The superscripts f and s describe the fluid medium and the solid medium, respectively. At the borehole wall (r = a), there are three boundary conditions:

· Continuity of the radial displacement

(

, ,

)

(

, ,

)

total total

rf z rs z

u% a k w =u% a k w , (2.80)

· Continuity of the normal stress

(

, ,

)

(

, ,

)

total total

f z rrs z

p a k w t a k w

-% = % , (2.81)

· Vanishing of the shear stress

(

, ,

)

0

total rzs a kz

t% w = . (2.82)

By considering the formal solutions of the acoustic field and elastic field, and the boundary conditions (2.80) – (2.82) at the borehole wall (r = a), the unknown constants A, B and C are determined by the following matrix equation:

(35)

CHAPTER 2 – General formulation of the source-borehole configuration 11 12 13 1 21 22 23 2 31 32 33 3 m m m A s m m m B s m m m C s æ öæ ö æ ö ç ÷ç ÷ ç ÷₌ ç ÷ç ÷ ç ÷ ç ÷ç ÷ ç ÷ è øè ø è ø , (2.83)

in which the terms due to the direct field are given by

(

)

(

)

1 2 2 3 1 , , , , , , 0. direct r z f direct z s p a k s p a k s w r w w = ¶ = -= % % (2.84)

The terms of the coefficient matrix are given by ( ) ( )

( )

( ) ( )

( )

( ) ( ) ( )

( )

( ) ( )

( )

( ) ( )

( )

( ) ( )

( )

( ) 11 2 1 12 1 13 1 21 0 2 2 22 0 1 23 0 1 31 32 1 2 2 33 1 , , , , 2 2 , 1 2 , 0, 2 , . f f r r f P r S z r S r f r P P P P z r r r r P r S S z r S r r P z r S S z r r S r k m I k a m K k a ik m K k a k m I k a k k m k K k a K k a a k m i k K k a K k a ak m m i k K k a k k m K k a k r w m l m m m m = -= = = æ é - ù ö =ç_ç ê ú- ÷_÷ -ë û è ø é ù = - _ê + _ú ë û = = é ₊ ù = - ê ú ë û (2.85)

We can rewrite m32 and m33 as followed:

( ) ( )

( )

( ) ( )

(

)

( ) ( )

( )

( ) 2 32 1 2 2 33 1 2 , . S P P z r r r S S P S z r r r r m ik k k K k a m k k k k K k a = = - +

From the system of equations and by using the substitution method, we can define the different constants A, B and C

(36)

(

)

(

)

1 22 33 23 32 2 12 33 13 32 21 33 1 11 33 2 21 32 1 11 32 2 , , . s m m m m s m m m m A m m s m m s B m m s m m s C - - -= D - + = D -= D (2.86)

where D is the determinant and is given by

11 22 33 13 21 32 11 23 32 12 21 33

m m m m m m m m m m m m

D = + - - (2.87)

By substituting (2.84) and (285) into equations (2.86), (2.87), the solutions of the constants A, B and C are

(

)

(

)

(

) (

₍

₎

)

₍

(

₎

)

( ) ( ) ( ) ( ) 1 1 2 2 2 ( )2 ( ) 2 2 ( ) ( ) ( ) 2 0 0 2 ( )2 2 ( )2 2 ( ) ( ) 1 1 1 , ( ) ( ) 2 2 4 , P S P S z r r r r P direct s S direct z r P r P r S f r S P S r r S z r z r P S S r r A k K k a K k a k k V k k p p k V k K K k a k a k k k k aV K K k a k a w r w l r w m w = é_ë ù_û D éæ ö ì _- _¶ - -êç ÷ í ê î ëè ø ü ùï + + + _úý úï ûþ % % (2.88)

(

)

(

) ( )

(

)

(

)

2 ( )2 ( ) ( ) ( ) 1 2 ( ) ( ) 0 1 1 , [ ], S S P z r S r r z r f direct direct f f r r r k k B k K k a k k I k a p I k a p w r w + = D ¶ % - % (2.89)

(

)

(

)

(

)

(

)

( ) ( )2 ( ) 1 2 ( ) ( ) ( ) 0 1 1 2 , [ ] P S z P r r z r f direct f direct f f r r r r ik C k K k a k k I k a p I k a k p w r w = D ¶ % - % (2.90) with ( )

( )

( ) ( ) ( ) ( )

( )

( ) ( ) ( ) ( ) ( ) ( )

(

)

( )

( )( ) ( ) ( )

( )

1 1 2 2 2 2 0 1 2 2 2 0 0 2 2 2 2 1 1 2 2 4 . P S P S r r r r P z r f _{s S} f f P r r r P r S f r P S r r S S z r _P z r _S S r r K k a K k a k k k k V I k a I k a k k V k K k a K k a k k k k a V K k a K k a r w l r w m w é ù D = ë û éæ é - ù ö ì _êç _ë _û _÷ - + -í _êç _÷ î _êè_ë _ø ü ùï ú + + + _ý úï ûþ (2.91)