IADC/SPE 112623

Drillstring Solutions Improve the Torque-Drag Model

Robert F Mitchell, Halliburton

Copyright 2008, IADC/SPE Drilling Conference

This paper was prepared for presentation at the 2008 IADC/SPE Drilling Conference held in Orlando, Florida, U.S.A., 4–6 March 2008.

This paper was selected for presentation by an IADC/SPE program committee following review of information contained in an abstract submitted by the author(s). Contents of the paper have not
been reviewed by the International Association of Drilling Contractors or the Society of Petroleum Engineers and are subject to correction by the author(s). The material does not necessarily
reflect any position of the International Association of Drilling Contractors or the Society of Petroleum Engineers, its officers, or members. Electronic reproduction, distribution, or storage of any
part of this paper without the written consent of the International Association of Drilling Contractors or 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 IADC/SPE copyright.

Abstract

The only standard drillstring model in use today is the torque-drag model, and because of the simplicity and general
availability of this model, it has been used extensively for planning and in the field. Field experience indicates that this model
generally gives good results for many wells, but sometimes performs poorly.

In the standard torque-drag model, the drillstring shape is taken as the wellbore shape. Considering that surveys are taken
within the drillstring, this is an excellent assumption. However, given that the most common method for determining the
wellbore shape is the minimum curvature method, the wellbore shape is less than ideal, because the bending moment is not
smooth at survey points. This defect is dealt with by neglecting bending moment, but as a result of this assumption some of
the contact force will also be neglected.

If the wellbore shape is the drillstring shape, why not use analytic drillstring solutions to model the wellbore trajectory?
Before these solutions can properly interpolate the survey data, they require some modification. This paper shows how this
can be done in a simple and efficient way, so that little of the calculation speed of the conventional model is lost. With this
new wellbore model, the simple torque-drag model becomes a full stiff-string formulation.

The paper gives a complete description of the wellbore trajectory calculation, including tests of circular, helical, and catenary
trajectories to verify the accuracy of the formulation. Several example problems with different inclinations and wellbore
curvatures are studied to compare the two torque-drag formulations. These studies give comparisons in execution time, drag
forces, and torques for the two models, and for the new formulation, the magnitude of the bending moments.

Introduction

Generally speaking, what do we mean by “torque-drag modeling?” Drag is the excess load compared to rotating drillstring
weight, which may be either positive when pulling the drillstring or negative while sliding into the well. This drag force is
attributed to friction generated by drillstring contact with the wellbore. When rotating, this same friction will reduce the
surface torque transmitted to the bit. It is useful to be able to estimate the friction forces when planning a well or doing post-
mortem analysis.

Analysis of these drillstring loads is done with drillstring computer models, and there have been many drillstring models
developed over the last 30 years (e.g. Walker, 1977). By far the most common method for drillstring analysis is the “torque-
drag” model originally developed by Dawson and Morehead (Johancsik, 1973) and put into differential equation form by
Sheppard (1987). Because of the simplicity and general availability of this model, it has been used extensively for planning
and in the field (e.g. Lesso, 1989). If any drillstring model could be called “standard,” this would be the one.

The torque-drag model formulation simplifies the general drillstring problem by assuming that the drillstring trajectory is the
same as the wellbore trajectory. Contact is assumed to be continuous. This assumed trajectory is certainly within inches of
the actual drillstring trajectory, and in fact, may be the same trajectory. Unfortunately, the most commonly used wellbore
trajectory model, the minimum curvature method, implies that bending moments are not continuous at survey points. The
2 IADC/SPE 112623

consequences of this assumption were probed by Mitchell and Samuel (Mitchell, 2007). The solution to this problem
proposed in this paper is to develop a new trajectory model that has sufficient smoothness to model the drillstring trajectory.

Wellbore Trajectory Models

Wellbore trajectory models are used for two distinct purposes. The first use is planning the well location, and consists of
determining kick-off points, build and drop rates, and straight sections needed to reach a specified target. The second use is to
integrate measured inclination and azimuth angles to determine a well’s location. Our interest here is the second purpose,
though angles determined for the first purpose can also be used in this model. Various trajectory models have been proposed,
with varying degrees of smoothness. The simplest model, the tangential model, consists of straight line sections. Thus, the
slope of this model is discontinuous at survey points. The most commonly used model is the minimum curvature model
(Taylor, 1972), which consists of circular arc sections. This model has continuous slope, but discontinuous curvature. In fact,
Taylor argues that a wellbore would not necessarily have continuous curvature. We argue, however, that since survey
measurements are made within the drillstring, the trajectory given by these measurements is the trajectory of the drillstring,
which must have continuity of bending moment (proportional to curvature).

If one wanted higher continuity in his trajectory model, Taylor proposed the use of cubic splines as an example. Before we go
further, a brief description of the use of splines is in order. Suppose we have a table of {x i , y i } , and we want to determine
intermediate values of y as a function of x. One simple method would be linear interpolation:

⎛ x −x ⎞ ⎛ x − xj ⎞
y( x ) = y j ⎜ j+1 ⎟ + y j+1 ⎜ ⎟ . . . . . .(1)
⎜x −x ⎟ ⎜x −x ⎟
⎝ j+1 j ⎠ ⎝ j+1 j ⎠
where we are interpolating between xj and xj+1. Now suppose that we want our interpolation to have smooth first and second
derivatives at the xj points. We might propose this interpolation:

y( x ) = y jf1 ( x ) + y′j′f 2 ( x ) + y j+1f 3 ( x ) + y′j′+1f 4 ( x ) . . . . . .(2)

where the functions fj have been devised so that:
y( x j ) = y j
y( x j+1 ) = y j+1
. . . . . .(3)
y′′( x j ) = y′j′
y′′( x j+1 ) = y′j′+1
In the classic cubic spline formulation (Press, 1992), the fj are cubic functions of x and we determine the unknown
coefficients y′j′ by requiring continuity of the first derivatives of y(x) at each xj. Note that the functions in equation (2) need
not be cubic functions. All they must do is satisfy equation (3).

The idea of using spline formulations to model wellbore trajectories is clearly not a new idea. For example, cubic splines and
tension-spines have been used to plan trajectories (Sampaio, 2006a and 2006b). Our goal is somewhat different. We wish to
find the wellbore trajectory from the survey data.
v
Given survey data, we can calculate the tangent vector t j at each survey point sj. Our formula for interpolating the tangent
vectors is:
r
v Tj (s)
t j (s) = v v
Tj (s) • Tj (s)
v v
Tj (s) = t jf1 j (s) + κ jn jf 2 j (s) . . . . . .(4)
v v
+ t j+1f 3 j (s) + κ j+1n j+1f 4 j (s)
v
where s is measured depth, κj is the curvature at sj, and n j is the normal vector at sj. This formulation has two purposes. The
first purpose is to satisfy the Frenet equation (Zwillinger, 1996) for a curve (by suitable choice of functions fij):
IADC/SPE 112623 3

v
d t (s) v
= κ(s) n (s) . . . . . .(5)
ds
The second reason is to insure that s is indeed measured depth. This requirement means:

du12 + du 22 + du 32 = ds 2
(an incremental change of position equals the incremental arc length) or, in terms of the tangent vectors:

2 2 2 v v
⎛ du1 ⎞ ⎛ du 2 ⎞ ⎛ du 3 ⎞ du du v v
⎜ ⎟ ⎜+ ⎟ ⎜+ ⎟ = • = t • t =1 . . . . . .(6)
⎝ ds ⎠ ⎝ ds ⎠ ⎝ ds ⎠ ds ds
It is easy to show that equation (4) satisfies this condition (see Appendix A). The unknowns in equation (4) are the normal
vectors and the curvatures. Appendix A gives the details for determining these unknowns.

The final missing piece of this problem is the choice of the interpolating functions fij. Since we wish to model drillstrings, the
best choice for interpolating functions should be solutions to actual drillstring problems.

The equation for the mechanical equilibrium of a weightless elastic rod with large displacement is (Nordgren, 1974, p. 778):
v r v
EIu iv − [(F − EIκ 2 )u′]′ = 0 . . . . . .(7)

where EI is the bending stiffness, F is the axial force (tension positive), and κ is the curvature of the rod. If we look at a small
interval of s, F and κ are roughly constant, so the solution to equation (7) is:

u (s) = c 0 + c1s + c 2 sinh(λs) + c3 cosh(λs)

when : EIλ2 = F − EIκ 2 > 0
u (s) = c 0 + c1s + c 2 sin(λs) + c3 cos(λs)
. . . . . .(8)
when : EIλ2 = EIκ 2 − F > 0
u (s) = c 0 + c1s + c 2s 2 + c3s 3
when : EIκ 2 − F = 0
where the c0-c3 are four constants to be determined. The third equation is simply a cubic equation, so cubic splines are a
candidate solution, even though they represent a special case of zero axial loads. Equation (8a) can be used to define what are
known as tension-splines, and we can use equation (8b) to define “compression” splines. Appendix B explains in detail how
this can be done.

One problem is that the λ coefficients are functions of the axial force, which we do not know until we solve the torque-drag
equations. In practice, λ tends to be small, so that the solution approximates a cubic equation. We can use the cubic
interpolation to approximate the trajectory, solve the toque-drag problem, then use the torque-drag solution to refine the
trajectory, iterating if necessary.

How does this trajectory model compare to the standard minimum curvature model? The two models were compared with
three analytic wellbore trajectories: circular arc, catenary, and helix. The comparisons of the displacements are shown in
Figures 1-3. Minimum curvature matches the circular arc exactly, since the method is a circular arc. Only one displacement is
shown for the helix, but is representative of the other displacements. Both trajectory models match the analytic solutions with
excellent results. The spline model was also used to calculate the rate of change of curvature for the catenary, Figure 4, and
the geometric torsion for the helix, Figure 5. The minimum curvature models predicts zero for both quantities. The spline
model determines both quantities accurately, though there is some end effect apparent in the geometric torsion calculation.
4 IADC/SPE 112623

Torque-Drag Calculations

Torque-drag calculations were made using the comprehensive torque-drag model developed in (Mitchell, 2007). The
equilibrium equations are integrated using the Bulirsh-Stoer method (Press, 1992). Otherwise, the only difference in the
solutions is the choice of the trajectory model.

For the first test we will consider the drag and torque properties of an idealized well plan based on Well 3 described in
Sheppard (Sheppard, 1986), see Figure 6. The fixed points on the model trajectory are as follows. The well is considered to
be drilled vertically to a KOP at a depth of 2,400 ft. The inclination angle then builds at a rate of 5°/100 ft. The target location
is considered to be at a vertical depth of 9,000 ft and displaced horizontally from the rig location by 6,000 ft. Drilled as a
conventional build-tangent well, this would correspond to a 44.5° well deviation. The model drillstring was configured with
372 feet of 6-1/2 inch drill collar (99.55 lbf/ft.) and 840 ft of 5 inch heavyweight pipe (50.53 lbf/ft.) with 5 inch drillpipe
(20.5 lbf/ft.) to the surface. A mud weight of 9.8 lbm/gal was used. In this example, a value of 0.4 was chosen for the
coefficient of friction to simulate severe conditions. Torque-loss calculations were made with an assumed WOB of 38,000
lbf. and with an assumed surface torque of 24,500 ft.-lbf.

Hook load calculated for zero friction was 192202 lbf for the circular arc calculation, and 192164 lbf for the spline model,
which compare to a spreadsheet calculation of 192203 lbf. The slight difference, 38 lbf, is due to the spline taking on a
slightly different shape (due to smoothness requirements) from the straight-line/circular arc shapes specified, which the
minimum curvature model exactly duplicated. Other than the slight difference in the spline trajectory, all other aspects of the
axial force calculations are identical between the two models. Tripping out, with friction coefficient 0.4, the hook load was
313474 lbf for the circular arc model and 319633 lbf for the spline model, for a difference of 6159 lbf. If we calculate from
the zero friction base line, this represents a difference of 5% in the axial force loading. Figure 7 shows the contact force
through the build section of the model. In overall magnitude, they are similar, although the spline formulation forces
smoothness that is not present in the minimum curvature model. The bending moments for the drillstring through the build
section are given in Figure 8. Note that minimum curvature does give a lower bending moment than the spline, but that the
spline results are smoother. With a surface torque of 24,500 ft-lbs, the torque at the bit was 15387 ft-lbs for the minimum
curvature model and 14961 ft-lbs for the spline model. This represents a 3% difference in the distributed torque between the
two models. Since this case has a relatively mild build rate, and since the build section was only about 8% of the total well
depth, the expectation was to see a relatively small effect from the spline formulation. Because the classic torque-drag
analysis has historically given good results, the agreement of the two models for this case verifies that the overall formulation
is correct.

For a more demanding example, the short-radius wellbore described in SPE 85328 (Grinde, 2003) was used. Figure 9 shows
the vertical and horizontal views of the end of the wellpath. The fixed points on the model trajectory are as follows. The well
is considered to be drilled vertically to a KOP at a depth of 2,800 ft. The inclination angle then builds at a rate of 5°/100 ft to
an angle of 25 degrees. The remainder of the trajectory is given in Figure 9. The model drillstring was configured with a 2-
7/8 inch motor, two 3-1/4 inch MWD tools, and 2-7/8 inch drillpipe to the surface. A mud weight of 10.1 lbm/gal was used.
In this example, a value of 0.4 was chosen for the coefficient of friction to simulate severe conditions. The build rate for this
example was 42°/30m, roughly ten times the build rate of the first case. In the abstract it was asserted that by neglecting
bending moment, some of the contact force would also be neglected. Figure 10 shows that this is indeed the case, for the
contact force for the spline model at the end of the build is four times that of the minimum curvature model. Figure 11 shows
the bending moment for this case, with the minimum curvature model again giving a lower bending moment than the spline,
but again the spline results are smoother.

Execution Times
v
The position vector u for the spline model has to be determined by numerical integration. This position vector is not needed
for the torque-drag calculations, but otherwise would be needed for the trajectory definition. If we do not include the position
vector determination, the execution timing is comparable between the two methods. Calculation of the position vector takes
about 15 times as much time as the minimum curvature method, but since this calculation only needs to be made once, it has
little effect on the overall execution time. All test cases executed in less than .10 second on a modern personal computer.

Conclusions

The classic torque-drag model, because of its simplicity and general availability, has been used extensively for planning and
in the field. Experience indicates that this model generally gives good results for many wells, but sometimes performs poorly.

In the standard torque-drag model, the drillstring shape is taken as the wellbore shape. However, use of the minimum
IADC/SPE 112623 5

curvature method the wellbore shape is less than ideal, because the bending moment is not smooth at survey points.

A more robust model for the wellbore shape has been developed using spline functions derived from drillstring solutions.
These solutions are not only better for torque-drag analysis, but are more characteristic of the actual trajectory measured by
surveys.

Sample calculations show that the spline model matches well with the minimum curvature model for analytic solutions of
wellbore trajectories. Comparison of the two models for wellbore trajectories with conventional build rates shows good
agreement. Since the classic model is known to be accurate for this type of well, this verifies that the spline model is
comparable to the classic model. For a more demanding case with higher build rates, the spline model shows superior
performance and predicts higher contact loads.

References

Grinde, Jan, and Haugland, Torstein 2003. Short Radius TTRD Well with Rig Assisted Snubbing on the Veslefrikk Field,
SPE 85328 presented at the SPE/IADC Middle East Drilling Technology Conference and Exhibition, Abu Dhabi, UAE,
(20-22 October).

Johancsik, C.A., Dawson, R. and Friesen, D.B. 1983. Torque and Drag in Directional Wells - Prediction and Measurement,
paper SPE 11380 presented at the IADC/SPE Drilling conf., New Orleans.

Lesso, W.G., Mullens, E. and Daudey, J. 1989. Developing a Platform Strategy and Predicting Torque Losses for Modeled
Directional Wells in the Amaulijak Field of the Beaufort Sea, Canada, SPE 19550 presented at the 64th Annual Technical
Conference and Exhibition, San Antonio, (Oct 8-11).

Mitchell, Robert F. and Samuel, Robello 2007. How Good is the Torque-Drag Model? paper SPE 105068 presented at the
IADC/SPE Drilling conf., Amsterdam (February 20-22).

Nordgren, R. P. 1974, On Computation of the Motion of Elastic Rods, Trans. ASME, Journal of Applied Mechanics,
(September), pp. 777-780.

Press, William H. et.al. 1992, Numerical Recipes in Fortran 77, Second Edition, Cambridge, England: Cambridge University
Press., pp. 107-110.

Sampaio, Jorge H. B. 2006a, Planning 3D Well Trajectories Using Cubic Functions, Trans. ASME, Journal of Energy
Resources Technology, Vol. 128, (December), pp. 257-267.

Sampaio, Jorge H. B. 2006b, Planning 3DWell Trajectories Using Spline–in–Tension Functions, JERT–05–1043, pp. 1-22.

Sheppard, M.C., Wick, C. and Burgess, T.M. 1986. Designing Well Paths to Reduce Drag and Torque, paper SPE 115463,
presented at the 61st Annual Technical Conference and Exhibition of the SPE, New Orleans, (Oct).

Taylor, Howard L. and Mason, C. Mack 1972. A systematic Approach to Well Surveying Calculations, SPE Journal,
(December), pp. 474-488.

Walker, B.R. and Friedman, M.B. 1977. Three-Dimensional Force and Deflection Analysis of a Variable Cross-Section
Drillstring, Trans ASME, Journal of Pressure Vessel Tech., (May), pp. 367-373.

Zwillinger, Daniel (editor). 1996. CRC Standard Mathematical Tables and Formulae, 30th Edition, Boca Raton, Florida:
CRC Press, pp 321-322.

Nomenclature
r
b = binormal vector
~
b = special binormal vector
E = Young’s elastic modulus (psf)
f ij (s) = interpolation function
6 IADC/SPE 112623

F = the effective axial force (lbf)

~
F = F − EIκ 2
moment of inertia (ft4)
vI =
i = unit vector in east direction
vE
i = unit vector in north direction
vN
i = unit vector in downward direction
rZ
n = normal vector
~
n = special normal vector
rs = measured depth (ft)
t = tangent vector
v
Tv = spline tangent vector function
u = position vector, (ft)
v
u 0j = initial position vector, increment j (ft)
{x j , y j} = set of numbers to be interpolated
y′j′ = the second derivative of y(x) at xj

αj = component of curvature in ~
n direction (ft-1)
~
βj = component of curvature in b direction (ft-1)
Δsj = s j+1 − s j (ft)
λj = coefficient in spline functions
r
εj = angle between n and ~
n
κ = wellbore curvature (ft-1)
ϕ = wellbore trajectory inclination angle
ϑ = wellbore trajectory azimuth angle
τ = geometric torsion of a curve (ft-1)
ξj = (s − s )/(s
j j+1 − sj)
f′ = df/ds
f iv = d4f/ds4

subscripts
j = the interval from sj to sj+1

Appendix A: The Spline Wellbore Trajectory

The normal method for determining the well path is to use some type of surveying instrument to measure the inclination and
azimuth at various depths and then to calculate the trajectory.

At each station i, inclination angle ϕi and azimuth angle ϑi are measured, as well as the course length Δsi =s i+1-si between
stations. These angles have been corrected to true north, if a magnetic survey, or for drift, if a gyroscopic survey. The survey
v
angles define the tangent t i to the trajectory at each station i, where the tangent vector is defined in terms of inclination ϕi
and azimuth ϑi in the following formulas:
v r
ti • i N = cos(ϑi ) sin(ϕ i )
v v
ti • i E = sin(ϑi ) sin(ϕ i ) . . . . . .(A-1)
v v
ti • i z = cos(ϕ i )
IADC/SPE 112623 7

If we knew how the angles ϕ and ϑ varied between stations, or equivalently, if we knew how the tangent vectors varied
between stations, then we could determine the trajectory by integrating the tangent vector:
v
v du i
ti = , so
ds . . . . . .(A-2)
v v sv
u i ( s ) = u i0 + ∫ t i ds
si

v v v v
Given tangent vectors t j and t j+1 and associated normal vectors n j and n j+1 , we can create a tangent vector interpolation
function connecting these vectors. First, we will need a set of interpolation functions fij(s), s in [sj, sj+1], with the following
properties:

df1 j (s j ) df1 j (s j+1 )

f1 j (s j ) = 1, = 0, f1 j (s j+1 ) = 0, =0
ds ds
df 2 j (s j ) df 2 j (s j+1 )
f 2 j (s j ) = 0, = 1, f 2 j (s j+1 ) = 0, =0
ds ds
. . . . . .(A-3)
df 3 j (s j ) df 3 j (s j+1 )
f 3 j (s j ) = 0, = 0, f 3 j (s j+1 ) = 1, =0
ds ds
df 4 j (s j ) df 4 j (s j+1 )
f 4 j (s j ) = 0, = 0, f 4 j (s j+1 ) = 0, =1
ds ds
There are a variety of functions that satisfy equations (A-3). If we define the following spline function Tj(ξ):
v v
Tj (ξ) = t jf1 j (s) + κ jn jf 2 j (s)
v v . . . . . .(A-4)
+ t j+1f 3 j (s) + κ j+1n j+1f 4 j (s)

we see that:
v r
Tj (s j ) = t j
v v
Tj (s j+1 ) = t j+1
v
dTj v . . . . . .(A-5)
(s j ) = κ jn j
ds
v
dTj v
(s j+1 ) = κ j+1n j+1
ds
The function Tj satisfies the Frenet equation:
v
d t (s) v
= κ(s) n (s) . . . . . .(A-6)
ds
for a tangent vector at s = sj and sj+1. However, Tj is not a tangent vector because it is not a unit vector. This can be corrected
by normalizing Tj:
r
v Tj (s)
t j (s) = v v . . . . . .(A-7)
Tj (s) • Tj (s)
8 IADC/SPE 112623

where it is easy to show that equation (A-6) is still satisfied. In order to evaluate the curvatures κj, we differentiate equation
(A-7) twice and evaluate at s = sj and sj+1:
v
d 2 t (s j ) v
• t j = − κ 2j
ds 2
v
d 2 t (s j ) v d 2f 2 j (s j ) v v d 2f 3 j (s j ) 2
v d f 4 j (s j )
• n j = κj + n j • t j+1 + κ j+1n j+1 • n j
ds 2 ds 2 ds 2 ds 2
v
d 2 t (s j ) v v v d 2 f 3 j (s j ) v v d 2f 4 j (s j )
• b j = t j+1 • b j + κ n
j+1 j+1 • b j
ds 2 ds 2 ds 2
v . . . . . .(A-8)
d 2 t (s j+1 ) v
2
• t j+1 = − κ 2j+1
ds
2
v
d t (s j+1 ) v d 2f 4 j (s j+1 ) v v d 2f1 j (s j+1 ) v d 2f 2 j (s j+1 )
• n j+1 = κ j+1 + n j+1 • t j + κ jn j • n j+1
ds 2 ds 2 ds 2 ds 2
v
d 2 t (s j+1 ) v v v d 2f1 j (s j+1 ) v v d f 2 j (s j+1 )
2

• b j+1 = t j • b j+1 + κ jn j • b j+1

ds 2 ds 2 ds 2
Using the Frenet formulae (A-6) and
v v v
dn (s)
= − κ(s) t (s) + τ(s)b(s)
ds
v . . . . . .(A-9)
2
d t (s) v v v
= − κ 2
(s ) t (s ) + κ′( s ) n (s ) + κ (s ) τ(s ) b (s)
ds 2
we can see that:
v
d 2 t (s j ) v
2
• t j = − κ 2j
ds
v
2
d t (s j ) v dκ
2
•nj = j
ds ds
v
d t (s j ) v
2

• b j = κ jτ j
ds 2
v . . . . . .(A-10)
d 2 t (s j+1 ) v
• t j+1 = − κ 2j+1
ds 2
v
d 2 t (s j+1 ) v dκ j+1
• n j+1 =
ds 2 ds
v
d t (s j+1 ) v
2

• b j+1 = κ j+1τ j+1

ds 2
The Frenet formulae are identically satisfied for (A-9)a and (A-9)d.
v
Before we can solve this set of equations for curvatures κj, we need a representation for the normal vector, n j , and the
v
binormal vector, b j . If the tangent vector is defined by the inclination angle, ϕ j , and the azimuth angle, ϑ j , in the following
way:
IADC/SPE 112623 9

⎡sin ϕ j cos ϑ j ⎤
v ⎢ ⎥
t j = ⎢ sin ϕ j sin ϑ j ⎥ . . . . . .(A-11)
⎢ cos ϕ j ⎥
⎣ ⎦

Then the Frenet equation (7) requires

⎡cos ϕ j cos ϑ j ⎤ ⎡− sin ϑ j ⎤

d v ⎢ ⎥ d
t j = ⎢ cos ϕ j sin ϑ j ⎥ ϕ j + ⎢⎢ cos ϑ j ⎥⎥ sin ϕ j ϑ j
d
ds ds ds . . . . . .(A-12)
⎢ − sin ϕ j ⎥ ⎢⎣ 0 ⎥⎦
⎣ ⎦
v
= κ jn j

From equation (15), we can determine the well known equation for the curvature κj :

2 2
⎛d ⎞ ⎛d ⎞
κ j = ⎜ ϕ j ⎟ + sin 2 ϕ j ⎜ ϑ j ⎟ . . . . . .(A-13)
⎝ ds ⎠ ⎝ ds ⎠

If we define:

⎡cos ϕ j cos ϑ j ⎤
~ ⎢ ⎥
n j = ⎢ cos ϕ j sin ϑ j ⎥
⎢ − sin ϕ j ⎥
⎣ ⎦ . . . . . .(A-14)
⎡− sin ϑ j ⎤
~ ⎢
b j = ⎢ cos ϑ j ⎥⎥
⎢⎣ 0 ⎥⎦
v ~ v
Note that t j ,~
n j ,and b j form a right-handed coordinate system at sj. We can define the normal vector, n j , and the binormal
v
vector, b j , by rotation through the angle ε j around the tangent vector:

v ~
nj = ~ n j cos ε j + b j sin ε j
v ~ . . . . . .(A-15)
b j = −~
n j sin ε j + b j cos ε j
v
Then n j is a unit vector consistent with Frenet equation (5), given:

1 d sin ϕ j d
cos ε j = ϕ j and sin ε j = ϑj . . . . . .(A-16)
κ j ds κ j ds

The variables κj and εj are not the most convenient choices because of the nonlinearity introduced by the sine and cosine
functions. An alternate choice would be:
10 IADC/SPE 112623

v ~
κ jn j = α j~
n j + β jbj
α j = κ j cos ε j
β j = κ j sin ε j . . . . . .(A-17)

κ j = α 2j + β 2j
βj
ε j = tan −1
αj

~
Equations (A-10) can be rewritten in terms of the vectors ~
n and b to give:
v
d 2 t (s j ) ~ d 2 f 2 j (s j ) ~ v d 2 f 3 j (s j ) ~
[ ]
2
~ ~ d f 4 j (s j )
• n j = α j + n j • t j+1 + n j • α n
j+1 j+1 + β b
j+1 j+1
ds 2 ds 2 ds 2 ds 2
v
d 2 t (s j ) ~ d 2f 2 j (s j ) ~ v d 2f 3 j (s j ) ~
[ ]
2
~ ~ d f 4 j (s j )
• bj = β j + b j • t j+1 + b j • α j+1n j+1 + β j+1b j+1
ds 2 ds 2 ds 2 ds 2
v . . . . . .(A-18)
d 2 t (s j+1 ) ~ d 2 f 4 j (s j+1 ) ~ v d 2f1 j (s j+1 ) ~ ~
2
~ d f 2 j (s j+1 )
• n j+1 = α j+1 + n j+1 • t j + n j+1 • (α j n j + β j b j )
ds 2 ds 2 ds 2 ds 2
v
d 2 t (s j+1 ) ~ d 2f 4 j (s j+1 ) ~ v d 2f1 j (s j+1 ) ~ ~
2
~ d f 2 j (s j+1 )
• b j+1 = β j+1 + b j+1 • t j + b j+1 • ( α n
j j + β b
j j )
ds 2 ds 2 ds 2 ds 2
2
v 2
Continuity of d t / ds at survey points requires for j=2, N-1:

[α ] ⎛ d 2f 2 j d 2 f 4 j−1 ⎞
[ ]
2 2
~ ~ ~ ~ d f 2 j−1 ⎜ ⎟ ~ ~ ~ ~ d f4 j
n • n + β n • b + α j⎜ + + α n • n + β n • b
ds 2 ⎟⎠
j−1 j j−1 j−1 j j−1 j+1 j j+1 j+1 j j+1
ds 2 ⎝ ds
2
ds 2
⎛ v d 2f v d 2f ⎞
n j • ⎜⎜ t j+1 23 j − t j−1 12j−1 ⎟⎟
=~
⎝ ds ds ⎠
. . . . . .(A-19)

[α ] d 2f 2 j−1 ⎛ d 2f 2 j d 2f 4 j−1 ⎞
[ ]
2
~ ~ ~ ~ ⎜ ⎟ ~ ~ ~ ~ d f4 j
b j • n j−1 + β j−1b j • b j−1 + βj⎜ 2 + + α j+1b j • n j+1 + β j+1b j • b j+1
ds 2 ⎟⎠
j−1
ds 2 ⎝ ds ds 2

~ ⎛⎜ v d f 3 j v d f1 j−1 ⎞⎟
2 2

= b j • ⎜ t j+1 2 − t j−1
⎝ ds ds 2 ⎟⎠

The set of equations (A-19) together with boundary conditions defined at the initial and end points forms a diagonally
dominant block tridiagonal set of equations that are relatively easy to solve. Note also that by solving for αj and βj, we have
also solved for dϕj/ds and dϑj/ds through equation (A-17), and further, there is no ambiguity about the magnitude of ϑj (±nπ)
in the definition of these derivatives.

We need expressions for the parameters κ, τ, and κ′ that appear in the torque-drag equilibrium equations. Recalling the Frenet
formulae (equations A-6 and A-9):
IADC/SPE 112623 11

v
d t (s) r
= κ(s)n (s)
ds
v v v
dn (s)
= − κ(s) t (s) + τ(s)b(s)
ds
v . . . . . .(A-20)
2
d t (s) v v v
2
= − κ 2 (s) t (s) + κ ′(s)n (s) + κ(s)τ(s)b(s)
ds
v
v d t (s) v r v
t (s) × = t (s) × κ(s)n (s) = κ(s)b(s)
ds
We find that:

d v d v
κ(s) = t j (s) • t j (s)
ds ds
d d v d2 v
κ(s) κ(s) = t j (s) • 2 t j (s) . . . . . .(A-21)
ds ds ds
d v ⎡v d v ⎤
2
κ(s) 2 τ(s) = 2 t j (s) • ⎢ t j (s) × t j (ξ)⎥
ds ⎣ ds ⎦

If κ is non-zero at a given point, then:

d v d v
κ(s) = t j (s) • t j (s)
ds ds
d v d2 v
t j (s) • 2 t j (s)
d
κ(s) = ds ds . . . . . .(A-22)
ds d v d v
t j (s) • t j (s)
ds ds
d v ⎡v d v ⎤
2
t (s) • ⎢ t j (s) × t j (s)⎥
2 j
τ(s) =
ds ⎣ ds ⎦
d v d v
t j (s) • t j (s)
ds ds

Appendix B: Drillstring Solutions As Interpolation functions

As stated in Appendix A, we will need a set of interpolation functions fij(s), s in [sj, sj+1], with the following properties:

df1 j (s j ) df1 j (s j+1 )

f1 j (s j ) = 1, = 0, f1 j (s j+1 ) = 0, =0
ds ds
df 2 j (s j ) df 2 j (s j+1 )
f 2 j (s j ) = 0, = 1, f 2 j (s j+1 ) = 0, =0
ds ds
. . . . . .(B-1)
df 3 j (s j ) df 3 j (s j+1 )
f 3 j (s j ) = 0, = 0, f 3 j (s j+1 ) = 1, =0
ds ds
df 4 j (s j ) df 4 j (s j+1 )
f 4 j (s j ) = 0, = 0, f 4 j (s j+1 ) = 0, =1
ds ds
12 IADC/SPE 112623

For example, the following cubic functions satisfy the requirements of equation (B-1):

f1 j (s) = 1 + (2ξ − 3)ξ 2

f 2 j (s) = ξ(ξ − 1) 2 (s j+1 − s j )
f 3 j (s) = (3 − 2ξ)ξ 2 . . . . . .(B-2)

f 4 j (s) = ξ 2 (ξ − 1)(s j+1 − s j )

s − sj
ξ=
s j+1 − s j
The cubic spline functions defined in equation (B-2) are not the only possible choices. An alternate formulation that has
direct connection to drillstring solutions is the tension spline:

[cosh(λ) − 1)][1 − cosh(λξ)] sinh(λ)[λξ − sinh(λξ)]

f1 j (ξ) = 1 + −
λ sinh(λ) + 2[1 − cosh(λ)] λ sinh(λ) + 2[1 − cosh(λ)]
⎧ [sinh(λ) − λ cosh(λ))][1 − cosh(λξ)] [λ sinh(λ ) + 1 − cosh(λ)][λξ − sinh(λξ)] ⎫
f 2 j (ξ) = ⎨ξ − − ⎬(s j+1 − s j )
⎩ λ2 sinh(λ) + 2λ[1 − cosh(λ)] λ2 sinh(λ ) + 2λ[1 − cosh(λ)] ⎭
[1 − cosh(λ))][1 − cosh(λξ)] sinh(λ)[λξ − sinh(λξ)]
f 3 j ( ξ) = +
λ sinh(λ) + 2[1 − cosh(λ)] λ sinh(λ) + 2[1 − cosh(λ)]
⎧[sinh(λ ) − λ)][1 − cosh(λξ)] [1 − cosh(λ)][λξ − sinh(λξ)] ⎫
f 4 j ( ξ) = ⎨ 2 + 2 ⎬(s j+1 − s j )
⎩ λ sinh(λ) + 2λ[1 − cosh(λ)] λ sinh(λ) + 2λ[1 − cosh(λ)] ⎭
s − sj
ξ=
s j+1 − s j
. . . . . .(B-3)
where λ is a parameter to be determined. In the context of beam-column solutions,

~
F
λ = Δs
EI . . . . . .(B-4)
~
F = F − EIκ 2 > 0
A similar solution for strings in compression is:

[cos(λ) − 1)][1 − cos(λξ)] sin(λ )[λξ − sin(λξ)]

f1 j (ξ) = 1 − −
λ sin(λ) − 2[1 − cos(λ)] λ sin(λ ) − 2[1 − cos(λ)]
⎧ [sin(λ) − λ cos(λ))][1 − cos(λξ)] [λ sin(λ) − 1 + cos(λ)][λξ − sin(λξ)] ⎫
f 2 j (ξ) = ⎨ξ + − ⎬(s j+1 − s j )
⎩ λ2 sin(λ) − 2λ[1 − cos(λ)] λ2 sin(λ) − 2λ[1 − cos(λ)] ⎭
[cos(λ) − 1)][1 − cos(λξ)] sin(λ)[λξ − sin(λξ)]
f 3 j ( ξ) = +
λ sin(λ) − 2[1 − cos(λ)] λ sin(λ) − 2[1 − cos(λ)]
⎧ [sin(λ) − λ )][1 − cos(λξ)] λ[cos(λ) − 1][λξ − sin(λξ)] ⎫
f 4 j ( ξ) = ⎨ − 2 + 2 ⎬(s j+1 − s j )
⎩ λ sin(λ) − 2λ[1 − cos(λ)] λ sin(λ) − 2λ[1 − cos(λ)] ⎭
s − sj
ξ=
s j+1 − s j
. . . . . .(B-5)
IADC/SPE 112623 13

where λ is a parameter to be determined. In the context of beam-column solutions,

~
−F
λ = Δs
EI . . . . . .(B-6)
~
F = F − EIκ 2 < 0

8000

7000
lateral displacement - analytic
vertical displacement - analytic
6000 lateral displacement - minimum curvature
vertical displacement - minimum curvature
lateral displacment - spline
5000 vertical displacment - spline
displacement (ft)

4000

3000

2000

1000

0
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
measure depth (ft)

Figure 1: Circular Arc Comparison Test

14 IADC/SPE 112623

9000

8000 lateral displacement - analytic

vertical displacement - analytic
lateral displacement - minimum curvature
7000
vertical displacement - minimum curvature
lateral displacement - spline
6000 vertical displacement - spline
displacement (ft)

5000

4000

3000

2000

1000

0
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
measured depth (ft)

Figure 2: Catenary Comparison Test

15

10

north - analytic
north - minimum curvature
5 north - spline
displacement (ft)

-5

-10

-15
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
measured depth (ft)

Figure 3: Helix Comparison Test

IADC/SPE 112623 15

0.000E+00

-2.000E-09

analytic
change of curvature (1/ft2)

spline
-4.000E-09

-6.000E-09

-8.000E-09

-1.000E-08
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
measured depth (ft)

Figure 4: Catenary Rate-of-Change of Curvature

1.40E-03

1.20E-03

analytic
1.00E-03
spline
torsion (1/ft)

8.00E-04

6.00E-04

4.00E-04

2.00E-04

0.00E+00
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
measured depth (ft)

Figure 5: Helix Torsion Test

16 IADC/SPE 112623

2400 ft
9000 ft

5o per 100 ft

6000 ft

Figure 6: Test Case 1 Wellbore

300

250

200 minimum curvature

contact force (lbf/ft)

spline

150

100

50

0
0 1000 2000 3000 4000 5000 6000
measured depth (ft)

Figure 7: Test Case 1 - Contact Forces

IADC/SPE 112623 17

3500

3000

2500

Minimum Curvature
bending moment (ft-lbf)

spline
2000

1500

1000

500

0
0 1000 2000 3000 4000 5000 6000
measured depth (ft)

Figure 8 - Test Case 1 Bending Moment

Figure 9: Test Case 2 - Short Radius Wellpath SPE85328

18 IADC/SPE 112623

50

45

40

constant curvature
35
spline
contact force (lbf/ft)

30

25

20

15

10

0
9700 9900 10100 10300
measured depth (ft)

Figure 10: Test Case 2 Contact Force

6000

constant curvature
5000
spline

4000
bending moment(lbf-ft)

3000

2000

1000

0
9700 9900 10100 10300
measured depth (ft)

Figure 11: Test Case 2 Bending Moment