SPE 84246 A Compendium of Directional Calculations Based on the Minimum Curvature Method S.J. Sawaryn, SPE, J.L. Thorogood, SPE, BP plc.
Copyright 2003, Society of Petroleum Engineers Inc. This paper was prepared for presentation at the SPE Annual Technical Conference and Exhibition held in Denver, Colorado, U.S.A., 5 – 8 October 2003. This paper was selected for presentation by an SPE Program Committee following review of information contained in an abstract submitted by the author(s). Contents of the paper, as presented, have not been reviewed by the Society of Petroleum Engineers and are subject to correction by the author(s). The material, as presented, does not necessarily reflect any position of the Society of Petroleum Engineers, its officers, or members. Papers presented at SPE meetings are subject to publication review by Editorial Committees of the Society of Petroleum Engineers. Electronic reproduction, distribution, or storage of any part of this paper for commercial purposes without the written consent of the Society of Petroleum Engineers is prohibited. Permission to reproduce in print is restricted to an abstract of not more than 300 words; illustrations may not be copied. The abstract must contain conspicuous acknowledgment of where and by whom the paper was presented. Write Librarian, SPE, P.O. Box 833836, Richardson, TX 75083-3836, U.S.A., fax 01-972-952-9435.
Abstract The minimum curvature method has emerged as the accepted industry standard for the calculation of 3D directional surveys. Using this model, the well’s trajectory is represented by a series of circular arcs and straight lines. Collections of other points, lines and planes can be used to represent features such as adjacent wells, lease lines, geological targets and faults. The relationships between these objects have simple geometrical interpretations, making them amenable to mathematical treatment. The calculations are now used extensively in 3D imaging and directional collision scans, making them both business and safety critical. However, references for the calculations are incomplete, scattered in the literature and have no systematic mathematical treatment. These features make programming a consistent and reliable set of algorithms more difficult. Increased standardisation is needed. Investigation shows that iterative schemes have been used where explicit solutions are possible. Explicit calculations are preferred because they confer numerical predictability and stability. Though vector methods were frequently adopted in the early stages of the published derivations, opportunities for simplification were missed because of premature translation to Cartesian coordinates. This paper contains a compendium of algorithms based on the minimum curvature method (includes co-ordinate reference frames, toolface, interpolation, intersection with a target plane, minimum and maximum TVD in a horizontal section, point closest to a circular arc, survey station to a target position with and without the direction defined, nudges and steering runs). Consistent, vector methods have been used throughout with improvements in mathematical efficiency, stability and predictability of behaviour. The resulting algorithms are also simpler and more cost effective to code and test. This paper describes the practical context in which each of the algorithms
is applied and enumerates some key tests that need to be performed.
Introduction The first reference to the minimum curvature directional 1 survey calculation method is credited to Mason and Taylor in 2 1971. In the same year, Zaremba submitted an identical algorithm that he termed the circular arc method. In the minimum curvature method, two adjacent survey points are assumed to lie on a circular arc. The arc is located in a plane which orientation is defined by the known inclination and direction angles at the ends. By 1985 the minimum curvature method was recognised by the industry as one of the most accurate methods but it was regarded as cumbersome for hand calculation3,4. The emergence of well trajectory planning packages to help manage directional work in dense well clusters increased its popularity. It was natural to use the same model for both the surveys and the segments of the well plan trajectories. Today, with the wide spread use of computers, computational power is no longer an issue and the method has emerged as the accepted industry standard. Industry Requirements Over the years, various algorithms based on the minimum curvature method have been published for the construction of increasingly complex trajectories and tasks such as interpolation. Since these algorithms have emerged piecemeal, they have tended to use different nomenclatures and mathematical techniques for their solution. The result of this piecemeal development is duplicated and inefficient computer code and a poor understanding of the engineering integrity of the systems. Safety and Business Criticality An undetected fault in the coding or use of directional surveying and collision scanning software has been classified as having the potential to cause property damage, 5 environmental damage, personal injury or loss of reputation . The integrity of these business and safety critical drilling systems is therefore a concern. Modern 3D imaging and directional scanning packages execute thousands of calculations for each task. Increased automation of the workflows associated with these tasks means that most calculations must be taken for granted and will pass unchecked. Sawaryn6 et. al., al., have described a process for managing these systems and identified a number of requirements related to the equations they use. Specifically, equations must be traceable back to the source documentation
2
SPE 84246
that must clearly explain their purpose, limitations and use. The general characteristics of the published algorithms can be assessed against these requirements. Angular Change Like other survey calculation methods, the minimum curvature algorithm was originally developed to calculate a well’s position from directional surveys. The spacing between the survey stations was normally 30 to 500ft. At that time with typical build rates, the total angle change over a 100ft course length would rarely be allowed to exceed 5 deg and the final inclination of most of these early wells was below 90 deg. When creating directional well plans the total angle change between adjacent stations in the plan may be considerably larger. These days, in designer wells the angular change between two adjacent points on a well plan may exceed 90 deg and the final inclination often exceeds 90 deg. One well7 is recorded as having reached 164.7 deg inclination. Many of the published algorithms do not contain an explicit definition of the maximum permitted angle change. The multiple solutions arising from periodicity of the trigonometric equations involved makes this a serious concern. Mathematical Behaviour The possibility of multiple solutions means the results of the calculations may not always be as intended, unless great care has been paid to their implementation. Some algorithms employ iterative schemes so that even if the scheme coverges, there is no guarantee that it converges to the correct solution. Ideally, iterative schemes should be accompanied by proof of convergence. At the very least they should be thoroughly tested over some specified range of variables. Additionally, there are cases where no solution exists and extra code is needed to trap this condition. Explicit expressions are more predictable and usually confer advantages in speed and maintainability of the computer code.
For certain values, for example in geometrically straight hole, expressions may be indeterminate. One solution adopted by Zaremba is to define a suitably small number at which the expression jumps abruptly to the asymptotic value2. However, this can give rise to random differences between software 6 packages . A better method is to develop series expressions that enable a smooth transition to be maintained.
geometric interpretations and several examples of its use are highlighted in this paper. Inconsistent nomenclature also leads to implemenational difficulties. Review shows the nomenclatures used in the literature are neither consistent with each other, nor consistent with accepted mathematical practice. An example is the definition of the normal vector that mathematical convention has pointing towards the centre of curvature. This is opposite to the convention used in the earlier drilling literature. Because of the expansion in directional drilling applications, symbols inevitably conflict. The Industry Steering Committee on Wellbore Survey Accuracy (ISCWSA) has proposed some standards8, but with limited scope. We conclude the SPE documentation standards9 associated with this subject area are no longer adequate and need revising.
Directional Calculations A consistent vector notation is used throughoutthis paper. This simplifies the development of the 3D equations and improves the clarity and presentation of the results. For convenience, the main vector operations are summarized in the appendix. In some cases, series expansions have been used to ensure the smooth transition of an expression to what would otherwise be an indeterminate form. The thresholds at which the series approximations should be used depend on the machine precision. The constants used in this paper assume calculations are good to at least 9 significant digits. The angle α subtending the arc may assume values such that 0 ≤ < π . Throughout, it is also assumed the start and end points of the arc are not coincident. Until such time as the standards are officially revised, we have chosen to maintain commonality with earlier papers on this subject and use a vector b pointing away from the centre of the arc. For comparison with mathematical texts, the normal vector is therefore – b. Reference Frames Co-ordinate Reference Frame The traditional reference frame for directional work uses North, East and Vertical coordinates that comprise a right handed set, Fig. 1. A point N , E ,V can be represented by the vector p, equation (2).
Common Constructs The absence of consistent mathematical methods and nomenclature may hide common constructs and potential simplifications in the coding of the algorithms. For example, an equation of the form C = A sin(α ) + B cos(α ) appears in many of the geometric constructions associated with the minimum curvature method. This equation can be solved explicitly for α and several mathematically equivalent forms exist. Zaremba proposed the form (1) that is used throughout this paper.
A ± ( A 2 + B 2 − C 2 ) α = 2 tan −1 B + C
1 2
N
E t
…(1)
When presented in this way, it can be seen that no real solutions exist if C 2 > A2 + B2. This inequality has simple
V Fig. 1
SPE 84246
p
3
N = E V
…(2)
A unit direction vector t can be represented in terms of the local inclination θ and azimuth φ , equation (3). The inclination and azimuth values can be calculated from the vector’s components using the expressions (4) and (5).
t
∆ N sin θ cos φ = ∆ E = sin θ sin φ ∆V cosθ
(∆ N 2 + ∆ E 2 ) θ = tan ∆V −1
1 2
…(3)
…(4)
∆ E ∆ N
φ = tan −1
…(5)
By using this reference frame an implicit assumption is made that the earth is flat. For moderate distances from the origin this assumption holds. For larger distances, such as those encountered in extended reach wells Earth’s curvature is important and corrections to the coordinates must be made. 10 Williamson and Wilson discuss the matter in detail .
cosθ cos φ h = cosθ sin φ − sin θ
…(6)
− sin φ r = cos φ 0
…(7)
0 v = 0 1
…(8)
Dogleg Severity Dogleg severity is a measure of the change in inclination and/or direction of a borehole, Fig. 3. The change is usually expressed in degrees per 100 ft of course length in oilfield units3 and degrees per 30m in metric units. Dogleg severity is used to determine stress fatigue in drill pipe, casing wear and casing design loads. It can also be a limiting factor in casing running and directional drilling operations. For the minimum curvature method, the expression for the dogleg severity takes the form (18000*α /π )/( D2− D1) in oilfield units. The difference in measured depths D2 – D1 between the points is referred to as the course length S 12.
N
Borehole Reference Frames Two reference frames are associated with the borehole, Fig. 2. The first frame is formed by the highside, rightside and tangent unit vectors h, r and t respectively. These form a right-handed, mutually orthogonal set. In curved hole, the second frame comprises the normal, binormal and tangent unit vectors – b, n and t respectively. These also form a right-handed, mutually orthogonal set. The angle between the highside vector h and normal vector – b is the toolface angle τ .
E t 1 t 2
h -b
v
r
t
Fig. 2
The highside, rightside and vertical unit vectors are represented in equations (6), (7) and (8). Expressions for the normal and binormal vectors – b and n can be found in the appendix (A-8 to A-12).
V Fig. 3
Most expressions found in the literature involve the calculation of an arc cosine. These have been used even though it is recognised that cosines of small angles are more 3 difficult to handle accurately than sines of small angles . Equation (9) developed by Lubinski does not make any assumption about the actual path of the wellbore, yet it is mathematically equivalent to the expression traditionally used in the minimum curvature method. An expression for tan(α /2) is readily developed from it.
θ − θ φ − φ 1 α = 2 sin −1 sin 2 2 + sin θ 1 sin θ 2 sin 2 2 1 2 2 …(9)
1 2
4
SPE 84246
Since the trigonometrical identity and computational advantages of (9) were recognised it is surprising it has not been adopted earlier. The dogleg severity can be related to both the radius of curvature and the curvature κ of the arc using the relationships (10).
κ =
1 R
=
α S 12
=
α * S *
…(10)
Survey Calculation Accurate determination of wellbore position is critical to well placement, collision avoidance, reservoir modeling and equity determination. Though the accuracy of the minimum curvature method is acknowledged, Stockhausen and Lesso11 showed that modern drilling practices could introduce systematic errors even with survey intervals as frequent as 100ft.
n12 b1
po
R p1
b2
t 1
p1 Fig. 5
Interpolation It is often required to identify the coordinates of a particular point, say p* on a trajectory Fig. 6. In all cases, the problem reduces to one of interpolation or extrapolation on an arc defined by the positions p1 and p2 and directions t 1 and t 2 of its end points. The algorithms presented here may be used for both functions. The interpolation may be driven by one of several parameters such as measured depth, subtended angle, inclination, azimuth, northing, easting or vertical ordinate.
n12
R
t 2
p1
The position of the next survey point p2 is calculated from p1 using (11), Fig. 4. The shape factor f (α) equals tan(α /2)/(α /2). Details are summarised in the appendix (A-6 to A-18).
p
2
= p 1 +
t 1 p*
sin θ 1 cos φ 1 + sin θ 2 cos φ 2 S 12 f (α ) sin θ 1 sin φ 1 + sin θ 2 sin φ 2 2 cos θ 1 + cosθ 2
p2
S 12
p2 Fig. 4
situation is of no practical concern. This may not be the case for adjacent points in a well plan trajectory, which may be separated by considerable distances.
…(11)
Straight Hole Conditions When α equals zero the shape factor is mathematically indeterminate, so for α < 0.02 radians the series expansion (12) should used instead12. The series is 13 presented in Horner form to minimise both the number of arithmetic operations and the propagation of errors.
α 2 α 2 31α 2 1 + 1 + 1 + …(12) f (α ) ≈ 1 + 12 10 168 18 α 2
There is a second possible solution of (9), equal to (2π – α ). The measured depth between the survey stations must be the same in both cases, implying the second solution has a greater curvature, Fig. 5. When calculating directional surveys, the density of survey stations, behaviour of the bottom hole assemblies and knowledge of the toolface settings means this
t * p2 t 2 Fig. 6
Before discussing each of these cases it is worth reviewing the properties of a circular arc. The subtended angle, inclination and azimuth of a point on a circular arc can be determined solely by the attitude of the circle in the coordinate reference frame. Knowledge of the size of the circle enables course lengths to be determined. Finally, it is only if a north, east or vertical ordinate is needed that the absolute position of the circle in the reference frame must be defined. Interpolation on Measured Depth The association of events or observations such as formation tops and overpulls with points on the wellbore is a common requirement. Interpolations on measured depth determined from the pipe tally are therefore the most common interpolation mode. If S * is the course length along the arc at which the properties are required, the relationships α * = (S */S 12)α and (α - α *) = (1 -
SPE 84246
5
S */S 12)α can be used to reduce the interpolation on measured depth to the interpolation on subtended angle (13). Interpolation on Subtended Angle The expression (13) enables the direction vector t * at a point p* on an arc to be determined solely from the direction vectors of its start and end points and the angle subtended from the first point p1 to the point of interest. Refer to the appendix (A-19 to A-27).
=
*
t
(
sin α − α *
) t + sin α
*
1
sin α
sin α
t 2
…(13)
Many of the algorithms presented in this paper involve the determination of the subtended angle α * as a first step. The relationship (13) provides a convenient means of determining the other parameters at the point once the subtended angle has been found. For example, once t * is known, the corresponding point p* can be calculated using the minimum curvature equation (11). Straight Hole Conditions When the subtended angle α equals zero, both terms in (13) are indeterminate. For α < 0.02 radians these terms which contain factors of the form 14 sin(cα )/sinα may be expanded in another Horner series (14). sin (cα ) sin α
≈ c + α
2
1 c 2 − 1 c 2 c − + α 2 c 7 + c 2 36 + 120 6 6 360
31 − 7 1 c 2 + α 2 c + c 2 + c 2 − 15120 2160 720 5040 127 − 31 7 − 1 c 2 + cα 2 + c 2 + c 2 + c 2 + 604800 90720 30240 362880 43200 …(14) It should be noted this interpolation mode is impossible when the subtended angle is identically zero, unless the constant c can be defined in some other way, such as by using the ratio of the measured depths. Interpolation on Azimuth Occasionally, it is necessary to truncate a well plan trajectory at a depth so that it is lined up with some specified direction. Firstly check the condition sin θ 1 sin θ 2 sin (φ 2 − φ 1 ) ≠ 0 to determine that the arc does
not lie in the vertical plane and that a solution exists. The subtended angle is then determined using (15).
*
α
= tan
−1
sin (φ 1 − φ * )sin α sin θ 1 * * sin(φ − φ 2 )sin θ 2 + sin (φ 1 − φ )sin θ 1 cosα …(15)
In the appendix it is shown that the expression (A-29) used to calculate the inclination θ * from φ * and α * is of the form C = A sin(θ *)+ B cos(θ *). The solution for θ * is determined using (1). In this case the constants A, B and C are given by (16),
(17) and (18). Choose the smallest root unless it is less than or equal to both θ 1 and θ 2 when the largest root should be chosen.
A = sin θ 1 cos(φ *
− φ 1 )
…(16)
B = cosθ 1
…(17)
C = cosα *
…(18)
As the subtended angle was originally determined from an azimuth on the arc a solution must exist. Finally, the position p* is determined using the minimum curvature equation (11). 2 2 2 Generally, if C > A + B the orientation of the arc is too shallow for the inclination at any point on it to reach the desired value. For details, see the appendix (A-28 to A-30). Straight Hole Conditions Azimuth varies linearly with measured depth in near straight hole conditions. For angles α < 10-4 radians, expression (19) should be used to determine the corresponding measured depth. If the azimuth values in either the numerator or denominator straddle north, the minimum angular difference should be used. The remaining properties can be determined by interpolating on this depth. *
S
φ * − φ 1 ≈ S 12 φ φ − 2 1
…(19)
Interpolation on Inclination The determination of measured depths corresponding to some inclination range appears in the automation of calculations for hole cleaning and rock mechanics. Firstly check that both points θ 1 and θ 2 are not equal to π /2. In this condition the arc would lie in the horizontal plane and no solutions are possible. In the appendix (A-31) it is shown that the subtended angle α * can be calculated using (1). In this case the constants A, B and C are given by the expressions (20), (21) and (22).
A = cosθ 2
− cosα cosθ 1
…(20)
B = sinα cosθ 1
…(21)
C = sin α cos θ *
…(22)
If C 2 > A2 + B2 then the orientation of the plane containing the arc is too close to the horizontal to enable the desired inclination to be reached and no solutions exist. The corresponding azimuth value can be determined from (23).
φ *
sin θ 1 sin φ 1 sin (α − α * ) + sin θ 2 sin φ 2 sin α * = tan −1 * * sin θ 1 cosφ 1 sin (α − α ) + sin θ 2 cos φ 2 sin α …(23)
Straight Hole Conditions Inclination varies linearly with measured depth in near straight hole conditions. For the subtended angle α < 10-4 radians, expression (24) should be used to determine the course length. The remaining properties can be determined by interpolating on this depth.
6
SPE 84246
*
S
θ * − θ 1 ≈ S 12 θ θ − 2 1
…(24)
Interpolation at a Plane A plane can be used to represent many geological features such as formation horizons and faults. In 3D visualisation tools, collections of interlocking planes are used to represent complex geological features. Imaginary planes can be used to represent lease boundaries or the north, east and vertical coordinate limits on directional plots. The point p* at which the well meets or crosses these features is of great practical interest, Fig. 7. The plane is uniquely defined by its normal vector m and any point px on it.
If C 2 > A2 + B2 the curvature of the arc is too large to intersect the plane and no solutions exist. If C 2 = A2 + B2 then the arc just touches the plane and there is only one solution. For C 2 < 2 2 A + B there are two intersections and in this case the subtended angles to both of them may be less than π . The two solutions correspond to the plus and minus signs in (1). This situation may be encountered when landing in a pay zone and the assembly is not building at the desired rate. Using this interpolation mode it is possible to determine if and where the bottom of the zone will be breached and the point at which the well is expected to re-enter. The lost production interval can then be calculated directly. Straight Hole Conditions When the subtended angle α -4 equals zero the above solution is indeterminate. For α < 10 radians a series expansion (29) is used to determine the measured depth. See the appendix (A-37 to A-44). No solution is possible if the line is parallel to the plane, when (m • t 1 )
R p1
equals zero.
px
t 1 p2
ζ =
S *
Fig. 7
There are five possible relationships between the arc and the plane. • An infinite number of intersections, when the plane containing the arc and the target plane are parallel. This condition should be tested first by establishing if both (m • t 1 ) and (m • t 2 ) are equal to zero.
• • •
Two intersections, when the arc completely cuts the plane. One intersection, when the arc just touches the plane. No intersections, when the curvature of the arc is too large. Lastly, if the curvature is so small then for practical purposes the problem is reduced to the intersection of a straight line with a plane.
In the appendix (A-32 to A-36) it is shown that the subtended angle α * to the plane can be calculated using equation (1). In this case the constants A, B and C are given by the expressions (25), (26) and (27).
A = (m • t 1 )sin α
…(25)
B = (m • t 1 ) cosα − (m • t 2 )
…(26)
C =
m • p
− p1 x S 12
α sin α
x
m
t 2 p*
•
m • p
…(27)
m • p
− p1 ζ 1 − (1 − ζ ) (m • t 1 ) 2 x
…(28)
…(29)
Orientation of the Target Plane The normal vector m defining the orientation of the target plane can be constructed in terms of the dip angle Θ and dip azimuth Φ of the plane (30).
− sin Θ cos Φ m = − sin Θ sin Φ cos Θ
…(30)
North, East and Vertical Interpolation The interpolation on north, east or vertical ordinates are particular cases of the more general expression. The values of the vectors px and m to use in each of these are listed in Table 1. For example, if the well is building and approaching an eastern lease line use the vectors listed under the heading East, with E * set equal to the numeric value of the eastern boundary limit.
p x
m
+ (m • t 1 ) cosα − (m • t 2 )
≈
− p 1 [(m • t 2 ) − (m • t 1 )] 2 S 12 (m • t 1 )
North
East
Vertical
N * 0 0 1 0 0
0 E * 0 0 1 0
0 0 * V 0 0 1
Table 1. Choice of vectors east and vertical ordinates.
px and m for
interpolating on north,
SPE 84246
7
Turning Point In horizontal wells, the trajectory is steered to remain in the reservoir section, successively building and dropping inclination to avoid breaching the bottom or top of the pay zone. The calculation of the turning points at which the well becomes horizontal is of interest from a reservoir perspective Fig. 8. Since these points represent the maximum and minimum vertical depths, they are also needed to determine the numerical range of axes in well plots.
t 2
p2 t *
Their algorithms were iterative and depended on the explicit determination of the centre of curvature of the arc, which causes problems for small angle changes. In the appendix (A48 to A-65) it is shown how the minimum curvature and minimum distance can be determined explicitly, without reference to the centre of curvature. For a fixed radius of curvature there are two possible trajectories from p1 to p3 shown by the solid and dashed lines in Fig. 10. All the trajectories lie in a plane. The solid line shows the minimum distance to the target. The other trajectory shown by the dashed line requires an angle change greater than π . Note that there are four mathematical solutions but that two of these are physically unrealisable because they would require the well to turn back on itself.
n12
p3
p1 p* t 1
p1 v
Fig. 8
This case is equivalent to interpolating on an inclination. In this special case, direct determination of the azimuth φ * at the turning point is possible (31).
φ *
cosθ 1 sin θ 2 sin φ 2 − cosθ 2 sin θ 1 sin φ 1 = tan −1 θ θ φ θ θ φ − cos sin cos cos sin cos 1 2 2 2 1 1
…(31) In deriving (31), it is assumed that θ 1 < θ 2 and so represents a minimum inclination. Should θ 1 > θ 2 then the inclination is a maximum. In this case the direction of t * will change by π and the azimuth becomes (φ * + π). The course length S * to the turning point can then be calculated from the relationship S * = S 12α */α using (13). Finally, its position vector p* can be calculated using the t 1 and t * vectors with the minimum curvature equation (11). See the appendix (A-45 to A-47).
Fig. 10
The changes in angle can be categorised according to the position of the target point p3 relative to the starting point p1 and the direction t 1. These categories are represented by the regions A to F in Fig. 11. Targets falling in region A can be hit with angle changes less than π /2 (90 deg.), B with angle changes less than π (180 deg.), C with angle changes less than 3π /2 (270 deg.) and D with angles changes less than 2π (360 deg.). Areas E and F cannot be hit at the prescribed build up rate.
C
Position at Target Defined Hogg and Thorogood15 described expressions for the minimum curvature and minimum distance from a point p1 on a wellbore with direction t 1 to a target p3, Fig. 9.
D F
p1
t 1
E
p1 p2 R
B
t 1
t 2
A
Fig. 11
p3 Fig. 9
t 3
The appropriate region is determined by calculating the distance ψ between the points p1 and p3, the perpendicular distance η from p1 in the direction t 1 and the distance ξ normal to t 1, equations (32), (33) and (34) respectively. Refer to Fig. 12.
8
SPE 84246
critical radius of curvature Rc given by (38). The corresponding angle change α c is given by equation (39).
t 1 p3
Rc
ψ 2
=
(
2
2 ψ
1 2
− η ) 2
…(38)
2 2 −1 (ψ − η ) α c = 2 tan η
1 2
p1
…(39)
The course length can be calculated as S 12 = Rcα c. Finally, t 2 is calculated from (37) with the tangent section length β equal to zero.
Fig. 12
ψ 2
2
= p 3 − p1
η = p 3 − p 1
…(32)
• t 1
…(33)
1
ξ = (ψ 2 − η 2 )2
…(34)
From Fig. 11 it can be seen that if
η ≤ 0 and ξ ≤ 2 R then
α ≥ π which is outside the imposed limit on the angle change. With the restriction 0 ≤ α < π only targets in areas A, B and E apply. Minimum Distance To Target If R
2
< ψ 2 + η 2 then the
tangent section length β is greater than zero. Because of the imposed restriction 0 ≤ α < π , the target must be in areas A or B, Fig. 11. The values of the tangent section length and angle change corresponding to the minimum distance to the target are calculated using (35) and (36).
β = ψ 2 − 2 R(ψ 2 − η 2 ) 1 2
1 2
…(35)
η − β 2 R − (ψ 2 − η 2 )
α = 2 tan −1
1 2
…(36)
Position and Direction at Target Defined Advances in surveying and geosteering techniques have enabled multiple targets to be penetrated as a matter of course. The targets may be at different geological horizons or different fault blocks in the same horizon. In these cases, the well’s trajectory must be lined up and its direction on entry to or exit from the target must be defined as well as its position. The trajectory can no longer be achieved with a simple build and hold profile. An additional curved section must be added and in the general case the trajectory is three-dimensional, Fig. 13. However, since the target can be specified so that the trajectory lies completely in a plane this calculation can also be used to design nudge profiles to increase well separation directly beneath a well cluster. Liu and Shi16,17 described an iterative scheme using coordinate transforms for the solution of the equations. Different calculations were used for each of the two arcs.
We offer an alternative iterative solution based on the geometrical symmetry of the problem and the minimum distance to target scheme described earlier in this paper. The advantage of this scheme is that only the subtended angles α 1 and α 2 of the arcs need to be determined each iteration. Inspection of Fig. 13 shows that in the general case, each of the two sets of points p1, p2, p3 and p2, p3, p4 are geometrically similar. This suggests applying the minimum distance to target algorithm alternately for each of the sets of points. t 1 p1
The course length is then calculated as S 12 = Rα . In straight hole ψ equals η and as long as R is finite (35) reduces to the distance between the points p3 and p1. The straight-hole case therefore degenerates to a straight line rather than an arc with an infinitely large radius, conferring stability to the calculation. Finally, t 2 is calculated using (37).
t 2
=
p 3 − p1 −
S 12 f (α )
S 12 f (α ) 2
2
t 1
p2 R1
R2
p3
…(37)
+ β
Minimum Curvature To Target If R 2
t 2
t 3
t 4 p4
Fig. 13
≥ ψ 2 + η 2 the
target must lie in region E and the build rate is insufficient to hit it. To hit the target, the build rate must be increased to the
Let p1, j and p4, j be the targets and α 1, j and α 2, j be the th subtended angles at the j iteration. To start the scheme, the
SPE 84246
9
subtended angle α 1,0 is calculated using the minimum distance to target algorithm between the points p1 and p4. This corresponds to the trajectory 1 in Fig. 14.
p1,1 p1,2
t 2 po
p2
p1
t * 2
n12
px
p1
3 1
p*
p4,1
t 1 b*
b1
Fig. 15
p4
p4,0
The closest approach is reached when the vector ( px – p*) is normal to the curve. The subtended angle α * to this point is calculated using expressions (43), (44) and (45).
Fig. 14
A new target p1,1 is calculated using expression (40).
α 1, j p1, j +1 = p1 + R1 tan t 1 2
α = p 4 − R2 tan 2, j t 4 2
…(41)
The next iteration is started by calculating the subtended angle α 1,1 using the points p1 and p4,1 and so on. Iteration continues until the desired precision is achieved. A suitable criterion for convergence is given by (42).
[
ε > (α 1, j +1 − α 1, j )
2
= p x − p1 • t 1
…(43)
η 2
= p x − p1 • t 2
…(44)
α *
η cosα − η 2 S 12 = tan −1 η 1 1 + α sin α
…(45)
…(40)
The subtended angle α 2,0 is now calculated for the second arc using the points p4 and p1,1. When performing the calculation, the direction of the borehole at p4 must be set to – t4 in order to calculate the correct angle. This corresponds to the trajectory 2 in Fig. 14. A new target p4,1 is calculated using expression (41) which completes the first iteration.
p 4, j +1
η 1
+ (α 2, j +1 − α 2, j )
2
] 1 2
…(42)
Four iterations are usually sufficient to reduce the error ε -5 below 10 radians. Note that on completing the calculation, the direction vectors of the second build and hold trajectory must be reversed before use. From a safety critical systems perspective, neither the above scheme, nor that proposed by Liu and Shi are completely satisfactory. In neither case is convergence proved and no definitive statement is made regarding the conditions under which no solutions exist, for example the radii of curvature are too large to hit the target. Further work is required. Closest Approach The calculation of the closest distance of a point px to a circular arc, Fig. 15 is the key to the construction of the normal plane collision avoidance diagram described by 18 Thorogood and Sawaryn .
The corresponding position p* on the arc is determined by interpolating on the subtended angle α * and then using the minimum curvature equation. Lastly, the distance is calculated from the magnitude of the vector ( px – p*). Mathematical and practical difficulties arise when the distance to the point is of the same magnitude or exceeds the radius of curvature of the arc. Restrictions on space preclude further discussion. Straight Hole Conditions Expression (45) is indeterminate -4 in straight hole, when α equals zero. For α < 10 radians, small angle approximations are used and expression (46) is used to calculate the course length S * directly.
S *
≈
S 12 η 1
η 1 − η 2
+ S 12
…(46)
The point p* is calculated by interpolating on measured depth assuming straight hole conditions. For these very small subtended angles the value of the shape factor f (α *) in the minimum curvature equation equals unity. Again, the distance is calculated from the magnitude of the vector (px – p*). Toolface Angle 16 Liu and Shi provided useful expressions for the toolface at the start and end of the arc in terms of the inclination and azimuth values at its ends. In the appendix (A-66 to A-79) it is shown how the general vector equation provided by Thorogood and Sawaryn18 may be expanded in terms of the inclination and azimuth values at any intermediate point on th e arc as well as at its ends (47).
10
SPE 84246
=
tan τ *
[
− φ * ) + sin θ 1 cos(α − α * )sin (φ * − φ 1 )] cos(α − α * ) cos θ 1 − cos α * cos θ 2
sin θ * sin θ 2 cos α * sin (φ 2
…(47)
The expressions (51) and (52) can be used to derive the equivalents (54) and (55) of Liu and Shi’s expressions. The representation of the azimuthal component is simplified.
At the ends, this expression reduces to those presented by Liu and Shi. The forms (48) and (49) have a small advantage because they are not singular when the inclination at either of the ends is zero. Dividing both numerator and denominator by a factor that can be zero causes the singularity. τ 1
τ 2
= tan
−1
sin θ 2 sin(φ 2 − φ 1 ) sinθ cosθ cos(φ − φ ) − sinθ cosθ …(48) 2 1 2 1 1 2
…(49) sin θ 1 sin (φ 2 − φ 1 ) = tan −1 sin θ 2 cosθ 1 − sin θ 1 cosθ 2 cos(φ 2 − φ 1 )
Straight Hole Conditions For α < 10-4 radians, small angle approximation to (47) is used giving (50). The toolface is constant over the arc and is undefined when α equals zero. If the azimuth values in (50) straddle north, the minimum angular difference should be used.
φ − φ 1 tan τ ≈ sin θ 1 2 θ 2 − θ 1 *
…(50)
Curvature The total curvature κ is a constant 1/ R on a circular arc, Fig. 16. Liu and Shi presented expressions for the inclination and azimuthal components of the curvature at any point on an arc in terms of the toolface angle at its start.
p* S r *
b*
t *
S sin Fig. 16
Using vector methods the inclination and azimuthal curvatures can be shown to be (51) and (52). For further details see the appendix (A-80 to A-84). These two components satisfy Wilson’s19 expression (53) for the total curvature κ . κ θ * = κ cosτ *
…(51)
sin τ * κ = κ * sin θ
…(52)
* φ
κ = (κ θ *2
+ κ φ *2 sin 2 θ * )
1 2
κ φ *
= κ
= κ
(
)
cos α − α cos θ 1 *
− cos α * cos θ 2
…(54)
*
sin α sin θ
sin θ 2 cos α * sin (φ 2
− φ * ) + sin θ 1 cos(α − α * )sin (φ * − φ 1 ) sin α sin θ *
…(55) Straight Hole Conditions Again, for α < 10-4 radians, small angle approximations to (54) and (55) are used and expressions (56) and (57) are used to calculate the curvatures. If the azimuth values in (57) straddle north, the minimum angular difference should be used.
* θ
κ
* φ
κ
≈
≈
θ *
− θ 1 *
S φ *
− φ 1
S *
≈
≈
θ 2
− θ 1
S 12 φ 2
− φ 1
S 12
…(56)
…(57)
Implementation and Testing The algorithms in this paper are presented in logical order with the later, more complex cases using results of earlier ones. It is recommended that the routines are coded and tested in this order. Once coded, a good test procedure is to calculate values in two different ways. For example, a point can be interpolated using each of measured depth, inclination and azimuth determined from the previous calculation in cyclic order. The results should be identical. The trajectory presented in Table 2 is constructed with both the wellbore position and direction at the target defined. This example can be used to test most of the algorithms presented here. Station numbers with alphabetical 8 suffixes indicate interpolated points. Williamson presents other trajectories that may be used as tests in both oilfield and metric units.
h* R
κ θ *
…(53)
The simple representation afforded by the circular arc construct allows for a consistent treatment of all the other associated mathematical operations. Representation of the wellbore path in a different form, for example a spline or polynomial would necessitate rederivation of all the constructs in this paper. This may not be a simple task.
Conclusions 1. Previously unpublished algorithms have been presented for small angle approximations associated with a circular arc, the determination of a turning point, a general expression for toolface angle and minimum distance to a target with and without the d irection at the target defined. 2. Vector methods are a useful tool for 3D directional calculations and often result in simpler expressions compared with other means. Their use is recommended.
SPE 84246
3. Because of trigonometric identities, many forms exist for any one expression. This can cause confusion. Some forms are computationally more efficient and reliable than others. 4. A standard nomenclature is required for all directional work that is compatible with other, related subject areas. Consistency with accepted mathematical conventions should also be reviewed. 5. Further work is needed on the point to target algorithm with both position and direction at the target defined. In particular, the conditions under which no solutions are possible should be identified. 6. Multiple solutions to the trigonometric equations exist. The range of the variables must be carefully stated and tested for. 7. Increasingly complex trajectory plans emerging from work on designer wells means segments of the well plan may exceed 180 deg and the alternative solutions may need to be considered. 8. Care should be taken to design expressions so that both the numerator and denominator are not divided by a term that would make it indeterminate at some point. 9. Representation of the wellbore path in a form other than a circular arc would necessitate rederivation of all the constructs in this paper.
Nomenclature b = Negative unit normal vector c = Ratio of the course lengths S */S 12 on an arc f = Geometrical shape factor h = Unit highside vector m = Unit vector normal to a plane n = Unit binormal vector p = Position vector in N , E ,V coordinates, L, ft r = Unit rightside vector t = Unit direction vector v = Unit vertical vector in N , E ,V coordinates, L, ft A = Constant B = Constant C = Constant D = Measured depth, L, ft E = Easting, L, ft N = Northing, L, ft R = Radius of curvature, L, ft S = Course length, L, ft V = Vertical, L, ft Greek Symbols α = Subtended angle, radians β = Length of the tangent section, L, ft γ = Angle between the binormal and vertical vectors, radians ε = Angular error tolerance, radians ζ = Substituted variable η = Substitution for a dot product, L, ft θ = Inclination angle, radians κ = Curvature, radians/L, radians/ft ξ = Substitution for a dot product, L, ft τ = Toolface angle, radians
11
φ = Azimuth angle, radians ψ = Substitution for a dot product, L, ft ∆ = A difference in a parameter Θ = Dip angle, radians Φ = Dip azimuth, radians Subcripts and Superscripts c = Critical value j = Iteration counter o = Centre of curvature of the arc x = Position of a defining point 1,2,3,4 = First, second etc. point, arc or property θ = Inclination component φ = Azimuthal component * = Component to be determined
Acknowledgments The authors thank BP plc. for permission to publish this paper. References 1. 2. 3. 4.
5.
6.
7. 8.
9. 10. 11.
12. 13.
14. 15.
16.
Mason, C.M., Taylor, H.L.: “Systematic Approach to Well Surveying Calculations”, SPE3362, June 1971 Zaremba, W.A.: “Directional Survey by the Circular Arc Method, SPEJ February 1973;Trans., AIME, 255. API Bull. D20, Directional Drilling Survey Calculation Methods and Terminology, first edition, December 31, 1985 Walstrom, J.E., Harvey R.P., Eddy H.D.: A Comparison of Various Directional Survey Models and an Approach To Model Error Analysis”, JPT August 1972 Sawaryn, S.J., Sanstrom, W., McColpin, G.: “The Management of Drilling Engineering and Well Services Software as Safety Critical Systems”, SPE 73893 presented at the International Conference on Health, Safety and Environment in Oil and Gas Exploration and Production, Kuala Lumpur, Malaysia, 20-22 March 2002 Sawaryn, S.J. et. al.: ”Safety Critical Systems Principles Applied to Drilling Engineering and Well Services Software”, SPE84152 presented at the Annual Technical Conference and Exhibition, Denver, Colorado, 5-8 October 2003. “Delivery Update - Autotrack Milestones”, Baker Hughes INTEQ (2002) 18 Williamson, H.S.: “Accuracy Prediction for Directional Measurement While Drilling”, SPEDC December 2000, Vol 15 No. 4 SPE Letter and Computer Symbols Standard, 1993 Williamson, H.S., Wilson, H.F.: “Directional Drilling and Earth Curvature”, SPEDC March 2000, Vol 15 No. 1 Stockhausen, E.J., Lesso, W.G.: “Continuous Direction and Inclination Measurements Lead to an Improvement in Wellbore Positioning”, SPE/IADC 79917, Amsterdam February 2003. Thorogood, J.L.: “Well Surveying Data”, World Oil, April 1986, 100 Knuth, D.E.: “The Art of Computer Programming Volume 2: Seminumerical Algorithm”, Third Edition, Addison Wesley, 1997, 485-488 Wolfram Research, “Mathematica 4 Standard Add-On Packages”, Wolfram Media Inc.,1999, 34 Hogg T.W., Thorogood J.L.: “Performance Optimisation of Steerable Systems”, ASME Energy Resources Technology Conference, New Orleans, January 1990 Liu, X., Shi, Z.: “Improved Method Makes a Soft Landing of Well Path”, OGJ October 22nd 2001
12
SPE 84246
17. Xiushan, L., Jun, G., "Description and Calculation of the Well Path with Spatial Arc Model," Natural Gas Industry, Beijing, Vol. 20, No. 5, 2000, 44-47. 18. Thorogood J.L., Sawaryn S.J.: “The Travelling Cylinder: A Practical Tool for Collision Avoidance”, SPE 19989 19. Wilson, G.J.: "An Improved Method for Computing Directional Surveys," J. Pet. Tech., Vol. 20, Aug. 1968, 87176. 20. Weatherburn, C.E.: Elementary Vector Analysis, G. Bell and Sons, London (1967) 21. Abramowitz and Stegun: Handbook of Mathematical Functions, Dover, New York (1972), Chap. 4, page 75
SI Metric Conversion Factors ft x 3.048* E-01 = m deg/100 ft x 0.984252 E-00 = deg/30 m *Conversion factor is exact Appendix Summary of Vector Methods 20 Weatherburn presents details of all the vector methods used in this paper. The key constructions are the dot (or scalar) product • and vector cross product × . Let a1, a2, a3 and b1, b2, b3 and c1, c2, c3 be the N , E ,V components of the vectors a, b and c respectively. Let α be the angle between the a and b vectors and n be a unit vector normal to both a and b. The dot and cross products are then described by expressions A-1, A-2 and A-3.
a • b = a1b1 + a2 b2
=
+ a3b3
a b cosα
…(A-1) …(A-2)
a × b = a b sin α n
…(A-3)
a × (b × c ) = (a • c ) b − (a • b) c
…(A-4)
a • (b × c ) = a1 (b2 c3
− b3c2 ) + a2 (b3 c1 − b1c3 )
+ a3 (b1c2 − b2 c1 )
…(A-5)
t 2
sin θ 2 cos φ 2 = sin θ 2 sin φ 2 cosθ 2
= t 1 × n12
…(A-9)
b2
= t 2 × n12
…(A-10)
Substituting n12 from A-8 into A-9 and A-10 and using the expression A-4 for the vector triple products gives A-11 and A-12.
b1
=
b2
=
t 1 cosα − t 2 sin α t 1
…(A-11)
− t 2 cosα
…(A-12)
sin α
The position p2 is calculated from p1 using A-13. Substituting b1 and b2 from A-11 and A-12 into A-13 gives A-14.
p 2
= p1 + R(b 2 − b1 )
p 2
= p1 + R
…(A-13)
(1 − cosα ) sin α
(t 1 + t 2 )
…(A-14)
Finally, recalling the trigonometric identity tan(α /2) = (1cosα )/sinα and that S 12 = Rα we obtain the now familiar expression for the minimum curvature equation A-15.
p 2
S 12
= p 1 +
α
α (t + t ) 1 2 2
tan
…(A-15)
Straight Hole Conditions For small angles, tan( z ) can be 21 expanded in a Taylor series , A-16. Using the first five terms of the expansion, the factor tan(α /2)/(α /2) can be put in Horner 12,13 form, A-17.
tan
α
2 α 2
z 3 3
+
2 z 5 15
+
17 z 7 315
+
62 z 9 2835
+ ...
α 2 α 2 31α 2 1 + ≈ 1 + 1 + 1 + 12 10 168 18 α 2
…(A-16)
…(A-17)
For small angles the shape factor can be treated as unity and A-15 reduces to A-18, which is recognised as the balanced tangential3 survey calculation method.
p 2
The binormal vector n12 and vectors b1 and b2 at the arc’s ends are given by the expressions A-8, A-9 and A-10.
…(A-8)
b1
…(A-6)
…(A-7)
sin α
tan ( z ) = z +
Minimum Curvature Zaremba2 presented the following derivation. Refering to Fig. 4, the direction vectors t 1 and t 2 at the arc’s ends are given by A-6 and A-7.
sin θ 1 cosφ 1 t 1 = sin θ 1 sin φ 1 cosθ 1
t 1 × t 2
=
n12
≈ p1 +
S 12 2
(t 1 + t 2 )
…(A-18)
Interpolation Refering to Fig. 6, interpolation involves determining the position p* at some point on the arc given a criterion. The corresponding direction vector at the point is t *, A-19.
*
t
sin θ * cos φ * = sin θ * sin φ * cosθ *
…(A-19)
SPE 84246
13
The binormal vector n12 can also be written in terms of the direction vectors at the start point and point at which the interpolation is to take place, A-20.
t 1 × t
*
n12
=
sin α *
…(A-20)
Equating A-8 and A-20 and taking the cross product of both sides of the equality with t 1 gives A-21.
t 1 × t 1 × t
*
=
sin α *
t 1 × (t 1 × t 2 ) sin α
…(A-21)
Using A-4 to expand the triple cross products and rearranging for t * gives A-22. *
t
= t 1 cosα − *
sin α * (t 1 cosα − t 2 ) sin α
…(A-22)
Multiply the numerator and denominator of the first term in A22 by sinα and collect terms in t 1 and t 2. After simplification, this gives the important relationship A-23 that is the foundation for all the interpolation formulae. *
t
=
sin (α − α
*
sin α
)t + sin α t 1
sin α
2
sin θ 2 cos φ 2 sin θ sin φ 2 2 sin α cosθ 2
sin α *
α
…(A-24)
t
=
S α S 12 *
t
= sin θ 1 cos(φ * − φ 1 )sin θ * + cosθ 1 cosθ *
If the arc lies in the vertical plane then (v • n12) equals zero and the solution is single valued. The vector n12 is given by A-8. Expanding the scalar triple product using A-5 gives the expression A-30.
sin α
t 1
+
S α S 12 t
sin
sin α
2
S * S * ≈ 1 − t 1 + t 2 S 12 S 12
− φ 1 )
…(A-27)
…(A-30)
= (cosθ 2 − cosα cosθ 1 )sin α * + sin α cosθ 1 cosα *
…(A-31)
Interpolation at a Plane Refering to Fig. 7, A-14 can be used to calculate the point p* from p1. The radius R can be expressed as S 12/α to give A-32. *
= p1 +
S 12 (1 − cosα * )
(t + t ) *
…(A-32)
1
*
α sin α
Now use A-22 to substitute for t * in A-32 to give A-33.
p …(A-26)
= 0
Interpolation on Inclination Extracting and rearranging the vertical component of A-24 gives A-31 which is of the form C = A sin(α *) + B cos(α *). Once α * has been found, the azimuth component is determined from A-28.
*
For small values of α the expression A-26 can be expanded in a Taylor series. Evaluation of the terms is tedious and computer assistance14 was used to establish the expression (14) presented in the body of the paper. For small angles, the simple expression A-27 can be used. *
cos α *
p
Substituting α * from A-25 into A-23 gives the expression A26 for t * in terms of the course lengths.
*
…(A-28) Next, expand the terms sin(α -α *) in both numerator and denominator of A-26 using the trigonometric identity sin(α α *) = sinα cosα * - sinα * cosα . The terms involving sinα * are then collected on the left hand side of the equals sign and terms involving cosα * are collected on the right hand side. The azimuth terms are then combined using the same trigonometric identity to give the expression (15) for α * used to interpolate on azimuth in the body of the paper. The 3 traditional form of the dogleg severity equation A-29 can be used to determine θ *. This is of the form C = A sin(θ *)+ B cos(θ *).
sin α * cosθ *
S 12
sin θ 1 cos φ 1 sin (α − α * ) + sin θ 2 cos φ 2 sin α *
sin α
…(A-25)
sin 1 −
=
sin θ 1 sin θ 2 sin (φ 2
Interpolation on Measured Depth Since the radius of the arc is fixed, the ratio of the subtended angles is identical to the ratio of the course lengths to the same points, A-25.
=
cos φ *
sin θ 1 sin φ 1 sin (α − α * ) + sin θ 2 sin φ 2 sin α *
…(A-29) …(A-23)
sin θ * cos φ * sin θ 1 cos φ 1 * sin α α ( ) − * * sin θ sin φ = sin α sin θ 1 sin φ 1 + cosθ * cosθ 1
S *
sin φ *
*
Substituting A-6, A-7 and A-18 in A-23 gives A-24.
α *
Interpolation on Azimuth Dividing the easting and northing components of A-24 eliminates sinθ * giving the expression A-28.
*
= p 1 +
S 12 (1 − cos α * ) *
α sin α
(1 + cosα ) t − *
1
sin α * (t 1 cosα − t 2 ) sin α
…(A-33)
Now take the dot product of A-33 with the normal vector m of the plane and rearrange slightly.
14
SPE 84246
(
m • p
*
− p 1 ) =
S 12 (1 − cos α
*
α sin α *
) (1 + cosα ) (m • t ) − *
1
sin α * ((m • t 1 ) cosα − (m • t 2 )) sin α
ζ =
…(A-34)
*
S
The equation of the plane20 is given by A-35, showing (m • px) equals (m • p*).
(
m • p x
− p * ) = 0
…(A-35)
Equation A-35 can be used to eliminate p* from A-34. After rearranging, the resulting expression A-36 is of the form C = A sin(α *) + B cos(α *).
m • p x − p1 α sin α S 12
+ (m • t 1 )cosα − (m • t 2 ) =
(m • t 1 )sin α sin α + [(m • t 1 )cosα − (m • t 2 )]cosα *
*
…(A-36) Straight Hole Conditions For small angles, equation A-36 is badly behaved. Small angle approximations must be used and the interpolation must be conducted with respect to measured depth. The expression A-18 can be used to calculate the position of p* from p1, A-37.
p
*
= p1 +
S * 2
(t + t ) *
− (m • t 1 ) ± (m • t 1 )(1 + 2ς ) = (m • t 2 ) − (m • t 1 ) S 12
(1 + z ) = 1 + − 2
2S 12
(t 2 − t 1 ) + S t 1 − ( p − p1 ) = 0
…(A-38)
Taking the dot product of A-38 with the normal vector m of the plane and substituting (m • px) for (m • p*) as before, gives a quadratic equation A-39 in the course length S *. A-39 has the solution A-40.
(m • t 2 ) − (m • t 1 ) 2S 12
S *2
+ (m • t 1 )S * − m • ( p x − p1 ) = 0 …(A-39)
2m • (t 2 − t 1 ) m • ( p x − p1 ) − (m • t 1 ) ± (m • t 1 )2 + S 12 * S = m • (t 2 − t 1 ) S 12
1 2
…(A-40) To simplify the manipulation, define a variable ζ according to A-41. Expression A-40 can then be written as A-42.
…(A-42)
8
+
z 3 16
−
z 4 128
+L
…(A-43)
…(A-44)
From A-44, it can be seen that the positive root must be chosen so that the expression degenerates to the straight-line case when ζ equals zero. Factoring A-44 gives the expression (29) in the body of the paper. Turning Point From Fig. 8 the vector t * can be written as A-45. The angle between the n12 and v vectors is γ. *
*
ς 2 ς 3 − (m • t 1 ) ± (m • t 1 )1 + ς − + 2 2 S * ≈ (m • t 2 ) − (m • t 1 ) S 12
t
*
z z 2
1 2
=
Sustituting A-27 into A-37 gives A-38.
S *2
1 2
…(A-41)
Since 2ζ is much smaller than unity, the square root may be expanded in a series using the first four terms of the binomial expansion21, A-43 to give the expression A-44.
…(A-37)
1
− p1 [(m • t 2 ) − (m • t 1 )] 2 S 12 (m • t 1 )
m • p x
n12 × v sin γ
…(A-45)
Using A-8 for n12 and expanding the resulting vector triple product gives A-46. *
t
=
− (v • t 2 ) t 1 + (v • t 1 ) t 2 sin γ sin α
…(A-46)
Note that (v • t 2) equals cosθ 2 and (v • t 1) equals cosθ 1 and that θ * equals π /2 at the turning point. Using these values and expressions A-19, A-6 and A-7 for t *, t 1 and t 2 in A-46 gives A-47.
cos φ * sin θ 2 cos φ 2 cos θ * 1 sin φ = sin γ sin α sin θ 2 sin φ 2 − 0 cosθ 2
sin θ 1 cos φ 1 sin θ sin φ 1 1 sin γ sin α cosθ 1 cosθ 2
…(A-47)
Dividing the easting by the northing components gives expression (31) in the body of the paper for the azimuth φ * of the turning point.
SPE 84246
15
Position at Target Defined In the most general case, the target p3 can be hit with one curved section of radius R and one straight section of length β , Fig. 9. Using A-15 the position p3 of the target can be written as A-48.
α = p1 + R tan (t 1 + t 2 ) + β t 2 2
p 3
…(A-48)
=
2
β
η 2
p 3 − p 1
• t 1 = R sin α + β cosα
2
…(A-49)
= 2 R 2 (1 − cosα ) + 2 β R sin α + β 2 …(A-50)
2
The substitution of β from A-53 and β from A-60 into A-54 results in A-61.
+
…(A-51)
• t 1
…(A-52)
Making the substitution A-51 in A-50 and A-52 in A-49 gives A-53 and A-54. The variable ξ is determined using Pythagoras’s theorem, Fig. 12.
η = R sin α + β cos α
…(A-53)
= 2 R 2 (1 − cosα ) + 2 β R sin α + β 2
…(A-54)
Multiplying A-53 by 2 R gives A-55. Rearranging A-54 gives A-56.
2 R sin α + 2 β R cosα = 2 Rη
…(A-55)
2 β R sin α − 2 R 2 cosα = ψ 2
…(A-56)
2
− β 2 − 2 R 2
Squaring both A-55 and A-56 and adding the results 2 eliminates α and results in a quadratic in β , A-57. Solving the quadratic gives an explicit expression for β , A-58
− 2 β 2ψ 2 + (ψ 4 − 4 R 2ψ 2 + 4 R 2η 2 ) = 0
β = ψ 2 ± 2 R(ψ 2 − η 2 ) 1 2
1 2
…(A-57)
…(A-58)
Inspecting A-58, the minimum curvature to target will occur when β equals zero and the radius equals Rc, A-59. Since the radius cannot be negative, the negative root of A-58 must be the correct one.
Rc
=
(η − R sin α ) +
− 2 Rη sin α + R 2 sin 2 α 2
cos α
ψ 2 2(ψ
2
1 2
− η ) 2
…(A-59)
To find the angle α , rearrange A-53 for β and square the result to give A-60.
…(A-61)
+
R 2 (1 + cos 2 α ) − 2 R 2 cosα
2
= p 3 − p 1
η = p 3 − p1
β 4
η 2
cos α
Multiplying through by cos2α and rearranging terms gives A62. Completing the square on the right hand side of A-62 gives A-63.
(ψ
2
ψ 2
2 R sin α
= 2 R 2 (1 − cos α ) +
ψ 2
ψ 2 cos 2 α = 2 Rη sin α (cosα − 1) + η 2
Let
ψ 2
…(A-60)
cos 2 α
Taking the dot product of A-48 with t 1 and also with itself gives A-49 and A-50 respectively.
p 3 − p 1
− 2 Rη sin α + R 2 sin 2 α
− η 2 )cos 2 α = [η sin α − R(1 − cosα )]2
…(A-62) …(A-63)
Taking the square root of A-63 and rearranging for R results in equation A-64 of the form C = A sin(α ) + B cos(α ). The solution of such an equation is given by (1).
R = η sin α + R ± (ψ
2
− η 2 ) cosα 1 2
…(A-64)
After substituting the values of the constants A, B and C and comparing the result with A-58 we obtain the important result 2 2 2 2 A + B – C = β . Setting β equal to zero provides a useful test to determine if the target can be hit at the specified radius of curvature. Finally, the expression for α /2 is given by A-65.
α η − β tan = 2 2 R − (ψ 2 − η 2 )
1 2
…(A-65)
Setting β equal to zero and R equal to Rc in A-65 gives the expression (39) in the body of the paper for the subtended angle α c for the minimum curvature to target. Toolface Angle The toolface angle is determined by the dot products between the -b, h and r vectors, Fig. 2.
(− b
*
• h * ) = cosτ *
…(A-66)
(− b
*
• r * ) = sinτ *
…(A-67)
Expressions A-66 and A-67 lead directly to the vector equation A-68 for toolface angle presented by Thorogood 18 and Sawaryn .
− b * • r * tan τ = * * − • b h *
…(A-68)
The rightside unit vector lies in the horizontal plane normal to both the v and t * vectors, A-69. Evaluating the expression gives (7) in the body of the paper.
16
SPE 84246
v × t
Inserting the expressions for the scalar triple products A-77 and A-78 into A-76 gives the final form for (-b* • r *), A-79.
*
=
*
r
*
…(A.69)
sinθ
The highside unit vector lies in the vertical plane normal to both the r * and t * vectors, A-70. Combining A-69 with A-70 and expanding the triple vector product gives A-71. Evaluating the expression gives (6) in the body of the paper.
h
*
h
*
= r * × t *
(− b
*
=
sin θ *
…(A-71)
b
= t × n12 *
sin θ 2 cos α sin (φ 2 *
sin α
*
d t
*
dS
…(A-72)
= −κ b *
…(A-80) *
cos α − α ) t 1 *
b
*
=
− cos α
*
t 2
sin α
The definitions for κ θ and
…(A-73)
*
d t
(− b
*
•h
*
2
*
1
1
2
sin α
− v sin θ * *
• r ) =
κ θ *
• (v × t * ) − cos(α − α * ) t 1 • (v × t * )
(
t 1 • v × t
) = − sinθ sinθ sin(φ − φ )
…(A-77)
t 2 • v × t
= sinθ 2 sinθ * sin(φ 2 − φ * )
…(A-78)
*
No.
*
1
= κ cosτ *
…(A-83)
the dot products of both sides with r * gives A-84 for
1
κ φ * .
Using A-67 and the expression A-79 for (-b* • r *) gives expression (55) in the body of the paper.
Using A-5, the scalar triple products become A-77 and A-78. *
…(A-82)
Eliminating the differential between A-80 and A-82 and taking
sin α sin θ * …(A-76)
*
= κ φ * r *
Using A-66 and the expression A-75 for (-b* • h*) gives the expression (54) in the body of the paper.
Taking the dot product of A-73 and A-69 gives an expression for (-b* • r *), A-76.
(− b
*
*
sin α sin θ
cos α * t 2
d t
the dot product of both sides with h* gives A-83 for κ θ .
*
*
1
…(A-81)
Eliminating the differential between A-80 and A-81 and taking
− cos α * cosθ 2 + cos(α − α * )cosθ 1 …(A-75) (− b • h ) = *
*
= κ θ * h *
sin θ * dS *
…(A-74) *
*
dS
) = − (t • t )t − (t • t )t • t cosθ *
κ φ * as the inclination and azimuthal
components of the curvature give expressions A-81 and A-82.
Taking the dot product of A-73 and A-71 gives an expression for (-b* • h*), A-74. Taking the dot products gives A-75. *
− φ * ) + sin θ 1 cos(α − α * )sin (φ * − φ 1 )
Curvature Refering to Fig. 16, Frenets20 formula for total curvature κ gives A-80.
The negative normal unit vector b* is normal to both the t * and n12 vectors. Combining A-8 and A-19 with A-72 gives A73. *
• r * ) =
…(A-79) The general expression (47) given in the body of the paper for toolface angle τ * is obtained by substituting A-75 and A-79 into A-68. For straight hole conditions, small angle approximations must be used.
…(A-70)
t cosθ * − v
*
κ φ *
= κ
sin τ *
sin θ *
…(A-84)
The signed values of A-83 and A-84 are commonly referred to as the build rate and walk rate respectively.
MD
Inc.
Azi.
N
E
V
DLS
(ft)
(deg)
(deg)
(ft)
(ft)
(ft)
(deg/100ft)
R
(deg)
(deg/100ft)
(deg/100ft)
(deg/100ft)
(ft)
1
702.55
5.50
45.00
40.00
40.00
700.00
-
38.14
1.57
12.89
2.00
2864.78
2
1964.57
29.75
77.05
154.78
393.64
1895.35
2.00
6.85
1.99
0.48
2.00
2864.78
2a
4250.00
29.75
77.05
408.84
1498.82
3879.60
0.00
-
-
-
-
-
3
5086.35
29.75
77.05
501.82
1903.25
4605.73
0.00
-135.72
2.15
4.22
3.00
1909.85
3a
8504.11
80.89
300.71
1967.04
1033.30
7050.00
3.00
-20.54
2.81
1.07
3.00
1909.85
3b
8828.04
90.00
297.31
2123.40
751.22
7075.71
3.00
-20.27
2.81
1.04
3.00
1909.85
3c
9151.97
99.11
293.92
2262.88
460.41
7050.00
3.00
-20.54
2.81
1.07
3.00
1909.85
4 9901.68 120.00 285.00 2500.00 -200.00 6800.00 3.00 -23.58 2.75 Table 2. Test case based on a trajectory with both its position and direction defined at the target
1.39
3.00
1909.85
