Modeling of drill string dynamics in directional wells for real-time simulation

Abstract

Simple and computationally efficient drill string models running real-time describing motion in all axes in directional wells are important for the implementation of closed-loop control and assisted monitoring during drilling operations. This paper proposes a new simplified three-dimensional model based on a parametric curve and lumped-parameter modeling, where Kane’s method is used to establish the equations of motion. Validation of the steady-state motion and convergence for the lumped model in vertical and horizontal alignment was compared with a finite-element model. The configuration and restoring forces show good results compared with finite-element analysis. Hence, the model demonstrate the axial contraction as a function of the body restoring forces being oriented to the inertial frame, inherently producing nonlinear coupled axial tension forces. The qualitative response of the model is confirmed in simulation case studies, being showcased by a deviated J-well configuration. Traveling block velocity and top drive torque are included as actuated inputs to analyze off-bottom friction and contact along the wellbore. The model is proposed to act as a virtual sensor for drilling directional wells.

Keywords

Offshore-drilling real-time Kane’s method coupled vibrations Runge–Kutta methods

1. Introduction

Thin hydrocarbon reservoir zones require drilling of complex directional wells, including long horizontal lengths for the well to be economically feasible. On the Norwegian Continental Shelf (NCS), several of the oil & gas zones are classified as thin (22–25 m¹). Advanced directional drilling methods have been applied to increase oil and gas recovery from these fields on the NCS, such as the Troll and Oseberg field.²

Oil & gas wells drilled at an angle from the straight vertical geometry are commonly denoted deviated wells. These wells can be preceded with multiple directional changes, including horizontal sections to reach challenging reservoir targets for maximizing hydrocarbon fluid entry into the production well. However, challenges with regards to drilling directional wells are factors such as substantial friction from drill string interacting with the wellbore wall during translatory motion and rotation, challenging hole-cleaning of cuttings, and drill string sections being stuck in the formation (known as key seating) in inclined sections.³

Predicting drill string motion in three-dimensional (3D) directional wells has been researched to some extent during the past years, typically using distributed models or finite element methods (FEM).^4,5 In Aarsnes and Shor,⁶ a distributed model of the torsional drill string dynamics for inclined wellbores was presented. Distributed friction was included, and the model was validated toward field data. A continuous Cosserat rod model was proposed in Goicoechea et al.,⁷ for arbitrary wellbores. The model equations were solved by implementation in a FEM environment. Valuable insight into the dynamics of drilling can also include heat-transfer and fluid dynamical models to predict wellbore instabilities (see, e.g., Wang et al.⁸). This is, however, outside the present paper scope of work. The FEM models are comprised of elements typically with two nodes having six degrees of freedom in each node. The complexity increases for a higher number of elements to approximate the distributed effects of a structure, hence, the model becomes computationally inefficient. The implication of this is necessary post-data analysis to match a model with field data.

The lumped-parameter method for modeling physical systems comprises of collecting distributed properties of, such as, mechanical and fluid dynamic systems into discrete lumps.⁹ These lumps are then describing the inertia and compliance effects of the model. Moreover, for mechanical vibrations, the natural frequencies tend toward the continuum properties when the number of lumped elements increases. Thus, when physical parameters in the model are lumped together, the dynamic model is commonly tuned to the application. Physical parameters such as friction coefficients and fluid parameters being obtained empirically can be calibrated online in real-time drilling performance optimization applications.

Simplified lumped-parameter models for drilling applications become relevant for control systems design and online calibration toward downhole bit motion and interactions with the wellbore. Real-time sensor data of downhole conditions can be obtained by using the wired-drill pipe technology.¹⁰ Moreover, combining this with a simulation model (a virtual copy of the physical process) updated with the measurements can aid the driller with closed-loop control, monitoring, and safe decision-making during drilling. Critical aspects to achieve the best possible performance of an online system using measurement data are such as the accuracy of the sensor measurements, signal transmission rate, dynamic model formulation, and the applied numerical solver.¹¹

For deviated wellbore formulations, a two-dimensional (2D) lumped model (LM) was proposed in Zhao et al.,¹² with estimated geometrical axial stiffness based on the change in wellbore inclination. Friction due to sliding was included assuming full contact with the wellbore wall, to analyze surge and swab pressures without rotation. The drill string model in Cayeux¹³ was also confined to motion in two axes, and included contact forces and friction at the tool-joint. The finite-difference method was applied and the model was argued to be suitable for real-time transient torque and drag analysis. A comprehensive study of a two-element 3D LM for a straight horizontal well was given in Xie et al.¹⁴ The equations of motion were formulated using Lagrange’s equations. A penalty-based contact law was assumed, and the interaction between bit and formation was included to model drilling conditions. A similar contact law was used in Qin et al.,¹⁵ where the Newmark iteration method used for time integration of a FEM model derived from spatial beam elements. The latter work also considered the segment of drill string confined inside the drilling riser, where its motion is influenced by the sea currents.

In this work, we develop a lumped multielement formulation of a drill string for arbitrary 3D wellbore configurations. The wellbore is simply a parametric curve in space, and the drill string motion is assumed to be a perturbation from the nominal wellbore curve.

The main contributions from this work are as follows:

A real-time oriented derivation of a drill string model for directional wells, where arbitrary wellbores can be simulated.

Application of Kane’s method for deriving a lumped parameter drill string model, in coordinate form. This method is, to the authors’ knowledge, novel in the genre of modeling drilling systems and contributes with its minimal set of equations compared with the traditional Lagrangian formulation. One effectively avoids lengthy partial differentiations as required in the more well-known Lagrange’s equations, critical for increasing number of bodies in the system (see, e.g., Kane and Levinson¹⁶ for a comparison of this method to well-established formulations).

Case studies for comparing the lumped-model to nonlinear FEM in horizontal configuration, and in case studies applying a wellbore contact model including a realistic drive configuration.

Introduction of Generalized- $α$ for numerical analysis of drill strings in contact with wellbore. In addition, this paper provides a introductory comparison of the well-known Runge–Kutta 4 to the Generalized- $α$ .

This work also includes a field-validated model from Gjerstad et al.¹⁷ for predicting fluid frictional shear stresses, to estimate surge and swab pressure. Extension and change of the wellbore during time is not considered in this paper.

Kane’s method introduces a minimal set of ordinary differential equations, suitable for real-time simulation.¹⁸ The correspondence to robotic manipulators is evident, and as such the drill string segment can further be developed as a link, where the link Jacobian can be given accordingly. Control synthesis and stability analysis in the robotics framework can then be applied to the model. The proposed model is demonstrated in an off-bottom scenario for a deviated well, and simulation performance is discussed.

The rest of the paper is structured as follows. Section 2 gives a preliminary introduction to the material applied in this paper, sections 3 and 4 describe the kinematics and dynamics of the proposed model, sections 5 and 6 give a study on dynamic convergence of the model and the transient response to common drilling scenarios, while section 7 includes a discussion of the real-time applicability. The paper is concluded in section 8.

2. Preliminaries

This section includes a brief definition of the curve parameterization representing a generic wellbore, and formulation of the equations of motion by using Kane’s method for a multibody system. The latter is used to develop a dynamic model of multiple lumped elements connected by elastic mass-less springs, in space.

The position of a point $i$ along a curve $C$ (see, e.g., Figure 1) is expressed relative to an inertial frame 0, given by its orthogonal unit vectors ${{\vec{e}}_{i}}$ , in coordinate free form, yielding

{\vec{p}}_{0 i} = c {\vec{e}}_{1} + f {\vec{e}}_{2} + h {\vec{e}}_{3}

(1)

and its coordinate form is given as $p_{0 i}^{0} = {[\begin{matrix} c & f & h \end{matrix}]}^{⊤}$ , where ${c (s), f (s), h (s)} \in C^{1}$ are parameterized by $s$ . In general, ${\vec{p}}_{ki}^{c}$ and $p_{ki}^{i / c}$ describe the position of $i$ relative to a frame $k$ , in the coordinates of $i$ being referenced to $c$ (superscripts), in coordinate free and coordinate form, respectively.

Figure 1.

A parametric curve in space.

Notice that $s$ resembles the natural paramterization used for a parametric curve. We assume that the curve does not change with time. Furthermore, the length of $C$ can be computed as

L_{C} = \int_{a}^{b} ‖ \frac{d {\vec{p}}_{0 i} (τ)}{d τ} ‖ d τ

(2)

where $‖ x ‖ = \sqrt{x_{1}^{2} + \dots + x_{n}^{2}}$ is the Euclidean norm.

A frame local to the origin of ${\vec{p}}_{0 i}$ is derived from the unit tangent, normal and bi-normal vectors for $C$ at the origin of ${\vec{p}}_{0 i} (s)$ . The unit tangent vector is given as

{\vec{c}}_{1} = \frac{{\vec{p}}_{0 i}^{'}}{| | {\vec{p}}_{0 i}^{'} | |}

(3)

where $'$ denotes the derivative with respect to $s$ . Furthermore, the unit normal vector is computed as the rate of change of the unit tangent with respect to the arc length, given as

κ {\vec{c}}_{3} = d {\vec{c}}_{1} / d s

(4)

where ${\vec{c}}_{1} \cdot d {\vec{c}}_{1} / d s = 0$ , where we have used the definition of curvature at $i$ given as $κ = ‖ d {\vec{c}}_{1} / d s ‖,$ being the rate at which the unit tangent vector is turning. In the case of zero curvature (completely straight sections) the unit normal vector are undefined due to $κ = 0$ . The unit normal and bi-normal vector given as

{\vec{c}}_{3} = {\vec{c}}_{1}^{'} / ‖ {\vec{c}}_{1}^{'} ‖

(5)

{\vec{c}}_{2} = {\vec{c}}_{1} \times {\vec{c}}_{3}

(6)

Finally, the representation of a skew-symmetric matrix is denoted by a hat notation yielding

\hat{x} = [\begin{matrix} 0 & - x_{3} & x_{2} \\ x_{3} & 0 & - x_{2} \\ - x_{2} & x_{1} & 0 \end{matrix}]

(7)

with $x \times y = \hat{x} y$ , and $\hat{x} = - {\hat{x}}^{⊤}$ , and Equation (7) will be used in the following sections.

2.1. Formulation of equations of motion

In this section, we introduce Kane’s formulation for establishing the equations of motion for a multi-body system. The equations of motion are then expressed in coordinate form.

The force and torque equilibrium of a rigid body $i$ , where the line of action is through the mass center, yields

{\vec{f}}_{i} + {\vec{f}}_{i}^{٭} + {\vec{f}}_{i}^{c} = 0, {\vec{T}}_{i} + {\vec{T}}_{i}^{٭} + {\vec{T}}_{i}^{c} = 0

(8)

where ${\vec{f}}_{i}$ , ${\vec{T}}_{i}$ are the external acting forces and torques, ${\vec{f}}_{i}^{٭}$ , ${\vec{T}}_{i}^{٭}$ are the inertia forces and torques, and ${\vec{f}}_{i}^{c}$ , ${\vec{T}}_{i}^{c}$ are vectors of constraint forces and torques.

To arrive at the equations of motion, the principle of virtual work and virtual displacements are used. The virtual displacement of a particle $k$ , which position is described by ${\vec{p}}_{k}$ , on body $i$ can be expressed as the sum of the partial contribution from each generalized coordinate $q_{r}$ , $r = 1, \dots, n_{q}$ , in the system. Hence, the virtual displacement is defined as

δ {\vec{p}}_{k} = \sum_{r}^{n_{q}} \frac{\partial {\vec{p}}_{k}}{\partial q_{r}} δ q_{r}

(9)

where $δ q_{r}$ is the virtual displacement of the generalized coordinate $q_{r}$ . Equation (9) is then a vector tangent to the position ${\vec{p}}_{k} (q, t)$ with $q$ being a vector of generalized coordinates.¹⁸ Moreover, assuming that each time derivative of $q_{r}$ , being ${\overset{\cdot}{q}}_{r}$ , is independent, the system will have $n_{q}$ degrees of freedom. In addition, to satisfy the constraint ${\vec{p}}_{k} = {\vec{p}}_{k} (q, t)$ , a set of constraint forces are utilized which will, according to the principle of virtual work, perform no work on the system, i.e.

\sum_{k}^{n_{p}} {\vec{f}}_{k}^{c} δ {\vec{p}}_{k} = 0

(10)

where $n_{p}$ is the number of particles in $i$ .

Multiplying the virtual displacements with Equation (8) yields

\begin{matrix} ({\vec{f}}_{i} + {\vec{f}}_{i}^{٭} + {\vec{f}}_{i}^{c}) δ {\vec{p}}_{i} + \\ ({\vec{T}}_{i} + {\vec{T}}_{i}^{٭} + {\vec{T}}_{i}^{c}) δ {\vec{θ}}_{i} = 0 \end{matrix}

(11)

where

δ {\vec{p}}_{i} = \sum_{r}^{n_{q}} \frac{\partial {\vec{v}}_{0 i}^{c}}{\partial {\overset{\cdot}{q}}_{r}} δ q_{r}, δ {\vec{θ}}_{i} = \sum_{r}^{n_{q}} \frac{\partial {\vec{ω}}_{0 i}}{\partial {\overset{\cdot}{q}}_{r}} δ q_{r}

(12)

where ${\vec{v}}_{0 i}^{c}$ is the linear velocity of the mass center of $i$ relative to the inertial frame 0, ${\vec{ω}}_{0 i}$ is the angular velocity of $i$ relative to 0, and $\partial {\vec{v}}_{0 i}^{c} / \partial {\overset{\cdot}{q}}_{r} = \partial {\vec{p}}_{0 i}^{c} / \partial q_{r}$ is obtained by first taking the time derivative of ${\vec{p}}_{0 i}^{c} (q, t)$ , being

{\vec{v}}_{0 i}^{c} = \frac{d {\vec{p}}_{0 i}^{c} (q, t)}{d t} = \sum_{r}^{n_{q}} \frac{\partial {\vec{p}}_{0 i}^{c}}{\partial q_{r}} {\overset{\cdot}{q}}_{r} + \frac{\partial {\vec{p}}_{0 i}^{c}}{\partial t}

(13)

and furthermore, the partial derivative of Equation (13) with respect to ${\overset{\cdot}{q}}_{r}$ .

Summing over the number of bodies $i$ in Equation (11) yields the d’Alembert principle of equation of motion, expressed as

\sum_{i}^{n} ({\vec{f}}_{i} + {\vec{f}}_{i}^{٭}) δ {\vec{p}}_{i} + \sum_{i}^{n} ({\vec{T}}_{i} + {\vec{T}}_{i}^{٭}) δ {\vec{θ}}_{i} = 0

(14)

where the principle of virtual work has been used to eliminate the constraint forces and torques. Inserting for Equation (12) for body $i$ in Equation (14) while keeping the index $r$ fixed, we get

\sum_{i}^{n} \frac{\partial {\vec{v}}_{0 i}^{c}}{\partial {\overset{\cdot}{q}}_{r}} ({\vec{f}}_{i} + {\vec{f}}_{i}^{٭}) + \sum_{i}^{n} \frac{\partial {\vec{ω}}_{0 i}}{\partial {\overset{\cdot}{q}}_{r}} ({\vec{T}}_{i} + {\vec{T}}_{i}^{٭}) = 0

(15)

by considering that the variations in $δ q_{r}$ are arbitrary.

In Kane and Levinson,¹⁹ the linear and angular velocities of $i$ are expressed according to quantities known as partial velocities. Furthermore, these acts as projection vectors for the generalized speeds used in Kane and Levinson,¹⁹ being a linear combination of the time derivative of the generalized coordinates. In this work we assume that the generalized speeds are equal to the generalized velocities. However, utilizing the generalized speeds may result in simpler kinematic and, hence, dynamic equations if formulated efficiently. The linear and angular velocities are given as

{\vec{v}}_{0 i}^{c} = \sum_{r = 1}^{n_{q}} {\vec{\tilde{v}}}_{ir}^{c} {\overset{\cdot}{q}}_{r} + {\vec{v}}_{t}, {\vec{ω}}_{0 i} = \sum_{r = 1}^{n_{q}} {\vec{\tilde{ω}}}_{ir} {\overset{\cdot}{q}}_{r} + {\vec{ω}}_{t}

(16)

where ${\vec{\tilde{v}}}_{ir}^{c} = \partial {\vec{v}}_{0 i}^{c} / \partial {\overset{\cdot}{q}}_{r}$ is the partial linear velocity, ${\vec{\tilde{ω}}}_{ir} = \partial {\vec{ω}}_{0 i} / \partial {\overset{\cdot}{q}}_{r}$ is the partial angular velocity. Note that the partial velocities are equal to the fractions in Equation (12), and act as projection vectors for the generalized velocities.

The general form of the equations of motion by Kane and Levinson¹⁹ denoted Kane’s equations, can be obtained from Equation (15) describing the dynamic equilibrium for body $i$

F_{ir} + F_{ir}^{٭} = 0

(17)

where $F_{ir}$ is all of the generalized active forces and torques acting on body $i$ , and $F_{ir}^{٭}$ is the generalized inertia forces and torques for body $i$ , in the direction of each generalized velocity ${\overset{\cdot}{q}}_{r}$ .

The projection of the inertia forces and torques on ${\overset{\cdot}{q}}_{r}$ can be expressed by the terms in Equation (15) yielding¹⁶

F_{ir}^{٭} = {\vec{f}}_{i}^{٭} \frac{\partial {\vec{v}}_{0 i}^{c}}{\partial {\overset{\cdot}{q}}_{r}} + {\vec{T}}_{i}^{٭} \frac{\partial {\vec{ω}}_{0 i}}{\partial {\overset{\cdot}{q}}_{r}}

(18)

{\vec{f}}_{i}^{٭} = - m_{i} {\vec{a}}_{0 i}^{c}

(19)

{\vec{T}}_{i}^{٭} = - {\vec{I}}_{i}^{c} \cdot {\vec{α}}_{0 i} - {\vec{ω}}_{0 i} \times ({\vec{I}}_{i}^{c} \cdot {\vec{ω}}_{0 i})

(20)

where $m_{i}$ is the mass of $i$ , ${\vec{a}}_{0 i}^{c}$ is the linear acceleration of the mass center, ${\vec{α}}_{0 i}$ is the angular acceleration of $i$ , and ${\vec{I}}_{i}^{c}$ is the central inertia dyadic of $i$ .

Performing a summation over $n$ for Equations (17)–(20) we get

\begin{matrix} \sum_{i}^{n} ({\vec{\tilde{v}}}_{ir}^{c} \cdot {\vec{f}}_{i} + {\vec{\tilde{ω}}}_{ir} \cdot {\vec{T}}_{i} - m_{i} {\vec{\tilde{v}}}_{ir}^{c} \cdot {\vec{a}}_{0 i}^{c} - \\ {\tilde{ω}}_{ir} \cdot ({\vec{I}}_{i}^{c} \cdot {\vec{α}}_{0 i} + {\vec{ω}}_{i} \times ({\vec{I}}_{i}^{c} \cdot {\vec{ω}}_{0 i}))) = 0 . \end{matrix}

(21)

A coordinate form representation of Equation (21) is given next. A projection matrix is formulated for body $i$ from the partial velocities in Equation (16), given as

J_{i} = [\begin{matrix} {\tilde{v}}_{i 1}^{i / c} & \dots & {\tilde{v}}_{i n_{q}}^{i / c} \\ {\tilde{ω}}_{i 1}^{i} & \dots & {\tilde{ω}}_{i n_{q}}^{i} \end{matrix}]

(22)

where ${\tilde{v}}_{ir}^{c}, {\tilde{ω}}_{ir} \in R^{3 \times 1}$ , and the dimensions of $J_{i}$ is $6 \times n_{q}$ . Hence, if the equation of motion are for a link, $J_{i}$ is the link Jacobian as derived in Egeland and Sagli.²⁰ The active and inertia forces and torques are then expressed as a column vector yielding

ℱ_{i}^{⋆} = [\begin{matrix} m_{i} a_{0 i}^{i / c} \\ ℐ_{i}^{i / c} α_{0 i}^{i} + {\hat{ω}}_{0 i}^{i} (ℐ_{i}^{i / c} ω_{0 i}^{i}) \end{matrix}]

(23)

ℱ_{i} = {[\begin{matrix} {f_{i}^{i}}^{⊤} & {t_{i}^{i}}^{⊤} \end{matrix}]}^{⊤}

(24)

where $a_{i}^{i / c}$ , $α_{0 i}^{i} \in R^{3 \times 1}$ are the coordinate form linear and angular accelerations of the mass center of $i$ , respectively, $ω_{0 i}^{i} \in R^{3 \times 1}$ is the angular velocity of $i$ and ${\hat{ω}}_{0 i}^{i} \in R^{3 \times 1}$ is the skew symmetric form as given in Equation (7), $I_{i}^{i / c} R^{3 \times 3}$ is the inertia dyadic, $f_{i}^{i}, t_{i}^{i} \in R^{3 \times 1}$ are the active forces and torques, all expressed in the coordinates of $i$ .

Premultiplying Equations (23) and (24) with the transpose of Equation (22) yields

\sum_{i}^{n} J_{i}^{⊤} F_{i}^{٭} = τ

(25)

where $τ = \sum_{i}^{n} J_{i}^{⊤} F_{i}$ is the generalized active forces and torques, and each row of $J_{i}^{⊤}$ is multiplied with the inertia vector yields the specific row associated with ${\overset{\cdot}{q}}_{r}$ , equal to Equation (21).

For systems exhibiting compliance and damping effects, such as a drill string, Equation (25) can be extended to

\sum_{i}^{n} J_{i}^{⊤} F_{i}^{٭} + D \overset{\cdot}{q} + K q = τ

(26)

where $K, D \in R^{n_{q} \times n_{q}}$ are the stiffness and damping matrix included by superposition, respectively.

3. Wellbore and drilling system representation

The wellbore and drill string kinematics are presented in this section in order to build a simulation model for experimental studies in drilling conditions. The drill string is in this paper assumed to be residing on the center of the wellbore curve $C$ , and can be illustrated as shown in Figure 2.

Figure 2.

Drill string system confined in a wellbore (wb). The nomenclature is DP—drill pipe, DC—drill collar. At the drill rig floor, a top drive rotates the drill string with to achieve pipe rotational velocity ${\overset{\cdot}{ϕ}}_{s}$ , and a hoisting system controls the position in the well by $u_{h}$ . The spiral connection between the black dots, that is, the discretized model elements, resembles massless springs.

The representation is based upon the work in Tengesdal et al.²¹ for drill string models in two dimensional space being discretized into lumped masses. The model in this paper is developed for 3D space. Each of the model elements $(i = 1, \dots, n)$ can either represent drill pipe (DP) or the drill string bottom-hole assembly (BHA). The BHA includes drill collars (DC) — heavier, thick-walled pipes with larger stiffness, the drill bit and tools to shape the wellbore. The element $n$ , i.e., the last discretized element is then the lowermost part of the BHA.

A set of assumptions is applied to the derivations followed in this paper, and are outlined below:

Assumption (1). We assume that the curve $C$ is a static parametric curve independent of time $t$ . Hence, frame $i$ with $^{i} {\vec{c}}_{i}$ , denotes a fixed frame in space.

Assumption (2). The drill string element has uniform density, and a symmetric cross-section.

Assumption (3). The wellbore is circular, defined by the diameter $d_{w}$ , with its center on the curve $C$ .

Assumption (4). The rotation about the normal and bi-normal axes of $b_{i}$ is neglected, since the angular displacement about the tangent axes are larger, simplifying the modeling.

Assumption (5). The distributed properties of the drill string assembly can be described by discrete masses with inertia, connected by massless springs, reducing the system into finite dimension.

The curve $C$ defines the nominal wellbore configuration, and at $t = t_{0}$ the drill string is initially located on $C$ . The system reference frames are shown in Figure 3. Frame 0 is further on defined with orthogonal unit vectors ${\vec{a}}_{i}$ and is denoted as the fixed base frame, with positive ${\vec{a}}_{3}$ axis point downward $({\vec{a}}_{3} = - {\vec{e}}_{3})$ , frame $i$ located on $C$ is defined by the unit tangent $^{i} {\vec{c}}_{1}$ , normal $i_{3}$ and bi-normal $i_{2}$ vectors described in section 2, and a local frame $b_{i}$ describes the dynamics of each drill string element. Vertical reach is defined along the ${\vec{a}}_{3}$ axis, and horizontal reach is defined along the ${\vec{a}}_{1}$ axis. In addition, the horizontal reach in the direction of ${\vec{a}}_{2}$ can be specified for the model.

Figure 3.

Wellbore profile and curvature. The parameters $L_{vr}$ and $L_{hr}$ denote the vertical and horizontal reach, respectively.

The drill string model can ultimately be described by the equations of motion similar to what was introduced in section 2 as

M (q) \overset{\cdot\cdot}{q} + C (q, \overset{\cdot}{q}) \overset{\cdot}{q} + D (q, \overset{\cdot}{q}) \overset{\cdot}{q} + K q = τ (q, \overset{\cdot}{q}, t)

(27)

where $M$ is the system inertia matrix, $C$ is the Coriolis matrix, $D$ is the system damping matrix from viscous friction of the fluid in the well, $K$ is the stiffness matrix, and $τ$ are the forces and torques acting on the particular element. The forces and torques will in this paper comprise of the applied forces from the topdrive, gravitational forces of the weight of the drill string (assuming no added mass effects), nonlinear internal restoring forces and torques to due to elasticity of the long drill string, and contact forces arising from friction and impact with the wellbore wall. Note that Equation (27) will in this paper be formulated on the basis of lumped parameter modeling.

Hence, the objective is to generalize a model with the above formulation, comprising of four degree-of-freedoms (three translatory and one rotational), coupling of rotational and lateral motion due to eccentricity in the elements, external forces due to contact and input from the rig, and internal forces describing the elasticity of the structure. Moreover, the lumped-element configuration presented in this paper is developed with non-rigidly attached elements as illustrated in Figure 2. Furthermore, the work in this paper is limited to off-bottom analysis.

3.1. Kinematics

The static wellbore curve $C$ is defined by a vector in the coordinates of $0$ relative to an inertial frame and the arc-length parameter $s$ as

p_{0 i}^{0} (s) = {[\begin{matrix} c (s) & f (s) & h (s) \end{matrix}]}^{⊤}

(28)

where the bold notation denotes a column vector, $i = 1, \dots, n$ , and the three functions $c (s), f (s)$ , and $h (s)$ describe the wellbore profile in space as a function of the arc-length $s$ . Furthermore, $C$ is defined as $s \in [0, 1]$ .

Suppose that a survey point, $s_{i}$ , is located at discrete points along the curve $C$ . In a survey point the dynamics of a drill string segment can be observed, where $s_{i} = i l_{i} / L_{w}$ , $l_{i}$ is the distance between each $s_{i}$ , and $L_{w}$ is the total wellbore length. The total length of the well is given by Equation (2).

From Figure 4, the body frame $b_{i}$ describes the perturbed drill string motion at the survey-point from the nominal point $p_{0 i}^{0} (s_{i})$ in the directions of $^{i} {\vec{c}}_{1}$ , $^{i} {\vec{c}}_{2}$ , $^{i} {\vec{c}}_{3}$ . In addition, each body frame is subject to rotation of about $^{i} {\vec{b}}_{1}$ by $ϕ_{i}$ .

Figure 4.

Illustration of the fixed frame for a point along $C$ and perturbed motion of local frame $b_{i}$ .

The center of mass can be shifted from the geometrical center $b_{i}$ due to a mass out-of-balance in the normal plane. The out-of-balance center is given by the constant eccentricities as ${ϵ_{1}}^{i} {\vec{b}}_{2}, {ϵ_{2}}^{i} {\vec{b}}_{3},$ where ${\vec{b}}_{j}$ with $j = 1, 2, 3$ are the orthogonal unit vectors of frame $b_{i}$ coinciding with frame $i$ in the undeformed case.

The origin of the element frame $b_{i}$ (the geometric center) taking into account the eccentricity, in coordinates of $i$ is written in coordinate free form as

\begin{matrix} {\vec{p}}_{i b_{i}}^{c} (t) = ξ_{i} (t)^{i} {\vec{c}}_{1} + ζ_{i} (t)^{i} {\vec{c}}_{2} \\ + χ_{i} (t)^{i} {\vec{c}}_{3} + {ϵ_{1}}^{i} {\vec{b}}_{2} + {ϵ_{2}}^{i} {\vec{b}}_{3} \end{matrix}

(29)

where superscript $c$ denotes the geometric center, $^{i} {\vec{c}}_{1}$ is the unit tangent vector of frame $i$ computed from Equation (3), $^{i} {\vec{c}}_{3}$ is the unit normal vector computed according to Equation (5), $^{i} {\vec{c}}_{2} =^{i} {\vec{c}}_{1} \times^{i} {\vec{c}}_{3}$ is the bi-normal, and $ξ_{i}, ζ_{i}$ , and $χ_{i}$ are the magnitudes in each of the directions. The curvature $κ$ at $s_{i}$ is then computed according to Equation (4). Noting that $^{i} {\vec{b}}_{2} = \cos (ϕ_{i})^{i} {\vec{c}}_{2} - \sin (ϕ_{i})^{i} {\vec{c}}_{3}$ , and $^{i} {\vec{b}}_{3} = \sin (ϕ_{i})^{i} {\vec{c}}_{2} + \cos (ϕ_{i})^{i} {\vec{c}}_{3}$ , we express the position in the coordinates of $i$ as

\begin{matrix} {\vec{p}}_{i b_{i}}^{c} (t) = ξ_{i}^{i} {\vec{c}}_{1} + (ζ_{i} + ϵ_{1} c_{ϕ_{i}} + ϵ_{2} s_{ϕ_{i}})^{i} {\vec{c}}_{2} + \\ (χ_{i} - ϵ_{1} s_{ϕ_{i}} + ϵ_{2} c_{ϕ_{i}})^{i} {\vec{c}}_{3} \end{matrix}

(30)

where $s_{ϕ_{i}} = \sin (ϕ_{i})$ , and $c_{ϕ_{i}} = \cos (ϕ_{i})$ . The position of $b_{i}$ relative to frame 0, in the coordinates of frame 0, given in coordinate form is given as

p_{0 b_{i}}^{0 / c} (s, t) = p_{0 i}^{0} (s) + R_{i}^{0} p_{i b_{i}}^{i / c} (t)

(31)

where $R_{i}^{0} = {r_{jk}}, R_{i}^{0} \in SO (3)$ is the rotation matrix from $i$ to $0$ , given by the dot product of the orthogonal unit vectors of frame $i$ and frame 0, where each matrix element is given as

r_{jk} = {\vec{a}}_{j} \cdot {\vec{c}}_{k}

(32)

where $j, k \in [1, 2, 3]$ . $SO (3) = {R | R^{3 \times 3}, R R^{⊤} = I, \det (R) = 1}$ denoted special orthogonal group.

3.1.1. Remark

For completely straight wells, we have ${\vec{c}}_{3}$ and ${\vec{c}}_{2}$ undefined since $κ = 0$ . Hence, for a strictly vertical section $\vec{p} = L_{vr} {\vec{a}}_{3}$ , $^{i} {\vec{c}}_{3}$ are aligned with the ${\vec{a}}_{1}$ axis, and for strictly horizontal wells $\vec{p} = s L_{hr} {\vec{a}}_{1}$ or $\vec{p} = s L_{hr} {\vec{a}}_{2}$ , $^{i} {\vec{c}}_{3} = - {\vec{a}}_{3}$ .

Since the curve $C$ is static, we assume that the angular velocity ${\vec{ω}}_{0 i}$ of $i$ relative to $0$ is zero, ${\vec{ω}}_{i b_{i}} = {\overset{\cdot}{ϕ}}_{i}^{i} {\vec{b}}_{1} = {\overset{\cdot}{ϕ}}_{i}^{i} {\vec{c}}_{1}$ , and we neglect the rotation about the normal and bi-normal axes. The linear and angular velocity of $b_{i}$ in $i$ , in the coordinates of $i$ in coordinate form are expressed as

v_{i b_{i}}^{i / c} (t) = [\begin{matrix} {\overset{\cdot}{ξ}}_{i} \\ {\overset{\cdot}{ζ}}_{i} - (ϵ_{1} s_{ϕ_{i}} - ϵ_{2} c_{ϕ_{i}}) {\overset{\cdot}{ϕ}}_{i} \\ {\overset{\cdot}{χ}}_{i} - (ϵ_{1} c_{ϕ_{i}} - ϵ_{2} s_{ϕ_{i}}) {\overset{\cdot}{ϕ}}_{i} \end{matrix}], ω_{i b_{i}}^{i} (t) = [\begin{matrix} {\overset{\cdot}{ϕ}}_{i} \\ 0 \\ 0 \end{matrix}]

(33)

and the linear velocity of the origin of $b_{i}$ relative to the base frame 0, is derived as

v_{0 b_{i}}^{0 / c} (t) = \frac{^{0} d p_{0 b_{i}}^{0 / c}}{d t} = R_{i}^{0} \frac{^{i} d p_{i b_{i}}^{i / c}}{d t}

(34)

since ${\vec{v}}_{0 i} = 0$ due to $(\partial {\vec{p}}_{0 i} / \partial s) \partial s / \partial t = 0$ , and ${\vec{ω}}_{0 i} = 0$ .

The generalized coordinates define the configuration of the system, comprising of

q = {[\begin{matrix} ξ_{1} \dots ξ_{n} & ζ_{1} \dots ζ_{n} & χ_{1} \dots χ_{n} & ϕ_{1} \dots ϕ_{n} \end{matrix}]}^{⊤}

(35)

where $q \in R^{n_{q} \times 1}$ , $n_{q}$ is the number of generalized coordinates. Time differentiation of $q_{r}$ in Equation (35) yields the generalized velocities ${\overset{\cdot}{q}}_{r}$ . Furthermore, partial differentiation of the linear and angular velocities in Equation (33) with respect to the generalized velocities ${\overset{\cdot}{q}}_{r}$ yields the partial velocities.¹⁹ These velocities define the projection of the external and inertial forces onto each ${\overset{\cdot}{q}}_{r}$ since the latter has no directional quantity.

A velocity vector can be expressed as

ν_{i b_{i}}^{i / c} (t) = [\begin{matrix} v_{i b_{i}}^{i / c} \\ ω_{i b_{i}}^{i} \end{matrix}]

(36)

and the partial linear and angular velocity associated with $b_{i}$ is derived as

v_{b_{i} r}^{i / c} = \frac{\partial v_{i b_{i}}^{i / c}}{\partial {\overset{\cdot}{q}}_{r}}, {\tilde{ω}}_{b_{i} r}^{i} = \frac{\partial ω_{i b_{i}}}{\partial {\overset{\cdot}{q}}_{r}}, r = 1, \dots, n_{q} .

(37)

from which we can define the projection Jacobian as

J_{i} = [\begin{matrix} {\tilde{υ}}_{b_{i} 1}^{i / c}, \dots, {\tilde{υ}}_{b_{i}, n_{q} .}^{i / c} \\ {\tilde{ω}}_{b_{i} 1}^{i}, \dots, {\tilde{ω}}_{b_{i}, n_{q} .}^{i} \end{matrix}] .

(38)

The vector in Equation (36) and its time-derivative are expressed in terms of the generalized velocities using Equation (38) given by Egeland and Sagli²⁰ as

ν_{i b_{i}}^{i / c} = J_{i} \overset{\cdot}{q}, {\overset{\cdot}{ν}}_{i b_{i}}^{i / c} = J_{i} \overset{\cdot\cdot}{q} + {\overset{\cdot}{J}}_{i} \overset{\cdot}{q} .

(39)

where $\overset{\cdot}{J} = (\partial J_{i} / \partial q) \overset{\cdot}{q}$ is describing the Coriolis effects if drill string eccentricity is taken into account.

4. Drill string dynamics

The distributed motion of the system is discretized into $n$ mass points, comprising of the lumped system inertia. Longitudinal mass-less springs representing the axial and torsional compliance, and rotational springs located at $i$ yield the lumped bending moment.

The configuration is shown for $i - 1$ to $i + 1$ in Figure 5. Furthermore, the equations of motion as presented in Equation (26) are utilized in this section to obtain a dynamic drill string model.

Figure 5.

The mass points located at $l_{i}$ on $C$ .

The drill string element inertia forces of the center of mass yield

{\vec{f}}_{i}^{٭} = - {\tilde{m}}_{i} {\vec{a}}_{i b_{i}}^{c}

(40)

where ${\vec{a}}_{i b_{i}}^{c} = {\overset{\cdot}{\vec{v}}}_{i b_{i}}^{c}$ , and ${\tilde{m}}_{i} = m_{i} + m_{f}$ is the mass of the element given by its dry weight and added mass, respectively, due to acceleration of mud around and inside the element. These are given as²²

m_{i} = ρ_{s} l_{i} π (d_{i, dp / dc}^{o^{2}} - d_{i, dp / dc}^{i^{2}}) / 4

(41)

m_{f} = ρ_{s} l_{i} π (d_{i, dp / dc}^{o^{2}} - d_{i, dp / dc}^{i^{2}}) / 4

(42)

where $d_{i, dp / dc}^{o^{2}}$ is the outer diameter for drill pipe or collar (notation dp/dc) for element $i$ , $d_{i, dp / dc}^{i}$ is the inner diameter, $ρ_{s}$ is the density of the pipe material, $ρ_{m}$ is the mud density, and $C_{a}$ is the added mass coefficient.

The inertia torques are given as

{\vec{T}}_{i}^{٭} = - {\vec{I}}_{b_{i}}^{i / c} {\vec{α}}_{i b_{i}}^{c}

(43)

where ${\vec{α}}_{i b_{i}}^{c} = {\overset{\cdot}{\vec{ω}}}_{i b_{i}}$ , and noting that due to Assumption (2) we have ${\vec{ω}}_{i, b_{i}} \times ({\vec{I}}_{b_{i}}^{i / c} {\vec{ω}}_{i b_{i}}) = 0$ . Furthermore, the central inertia dyadic of the element is given as

\begin{matrix} {\vec{I}}_{b_{i}}^{i / c} = I_{xx}^{ci} {\vec{b}}_{1}^{i} {\vec{b}}_{1} + I_{yy}^{ci} {\vec{b}}_{2}^{i} {\vec{b}}_{2} + I_{zz}^{ci} {\vec{b}}_{3}^{i} {\vec{b}}_{3} \\ I_{yy}^{c} = I_{zz}^{c} = \frac{m_{i}}{12} (3 (d_{i, dp / dc}^{o^{2}} + d_{i, dp / dc}^{i^{2}}) + l_{i}^{2}) \\ I_{xx}^{c} = \frac{m_{i}}{2} (d_{i, dp / dc}^{o^{2}} + d_{i, dp / dc}^{i^{2}}) \end{matrix}

(44)

where $I_{yy}^{c}, I_{xx}^{c}, I_{zz}^{c}$ are the moments of inertia about each principal axis.

The internal forces and torques contribute to restoring the structure to its original configuration and resist elastic deformation. Assuming that a mass-less spring element connecting $i - 1$ to $i$ and $i + 1$ represents the elasticity of the pipe section, we can express the tension-induced forces on element $i$ in vector form in frame $i$ , given as

\begin{matrix} f_{e, i}^{i} = k_{i}^{a} (1 - l_{i} / | | p_{b_{i} b_{i - 1}}^{i / c} | |) p_{b_{i} b_{i - 1}}^{i / c} \\ + k_{i + 1}^{a} (1 - l_{i + 1} / | | p_{b_{i} b_{i + 1}}^{i / c} | |) p_{b_{i} b_{i + 1}}^{i / c} \end{matrix}

(45)

where $p_{b_{i} b_{i - 1}}^{i / c}$ is the coordinate vector from $b_{i}$ to $b_{i - 1}$ , and $p_{b_{i} b_{i + 1}}^{i / c}$ is the coordinate vector from $b_{i}$ to $b_{i + 1}$ , both expressed in the coordinates of $i$ , and $k_{i}^{a} = E A_{i} / l_{i}$ the axial stiffness coefficient for the lumped element.

The restoring torques about the centroid axis are given as

{\vec{T}}_{e, i} = - \frac{G J_{i}}{l_{i}} (ϕ_{i} - ϕ_{i - 1})^{i} {\vec{c}}_{1} + \frac{G J_{i + 1}}{l_{i + 1}} (ϕ_{i + 1} - ϕ_{i})^{i} {\vec{c}}_{1}

(46)

where $G$ is the shear modulus, and $J_{i}$ is the $i th$ torsion constant (here equal to the polar moment of inertia, due to the symmetric cross-section).

The vector of tension forces and restoring torques imposed on element $i$ are then expressed as

F_{i}^{e} = {[\begin{matrix} f_{e, i}^{i} & t_{e, i}^{i} \end{matrix}]}^{⊤} .

(47)

The tension-induced forces comprise the axial stiffness of the drill string, and the coupling to the transverse motion along the assembly. To include the effect of bending moment, a linear approach by superposition on the generalized coordinates is taken. The forces on each element due to bending can be derived from the potential energy, yielding

τ_{b} = - \frac{\partial U_{b}}{\partial q}

(48)

where the term $U_{b}$ comprises the bending strain component. Assuming that the drill string cross-section remains plane in bending, we neglect the rotary inertia about $^{i} {\vec{b}}_{2}$ and $^{i} {\vec{b}}_{3}$ for the drill string element. Consequently, the element falls under the category of Euler–Bernoulli beam theory, with infinite shear stiffness.

From Figure 5, the angles formed from a change in lateral perturbation from nominal configuration at $i$ along $C$ can be expressed as

θ_{1, i} = \arcsin (\frac{ζ_{i + 1} - ζ_{i}}{l_{i}}), θ_{2, i} = \arcsin (\frac{χ_{i + 1} - χ_{i}}{l_{i}})

(49)

where $ζ_{i}$ and $χ_{i}$ are the local perturbations in the direction of ${\vec{c}}_{2}$ and ${\vec{c}}_{3}$ , respectively.

The lumped bending moment can be derived considering the potential energy of a rotational spring, located in each frame $i$ , given as

U_{b} = \sum_{i = 1}^{n - 1} \frac{k_{i}^{b}}{2} ({(θ_{1, i} - θ_{1, i - 1})}^{2} + {(θ_{2, i} - θ_{2, i - 1})}^{2})

(50)

where $k_{i}^{b} (E, I_{i}^{x / y}, l_{i})$ is the bending stiffness coefficient. Using the small angle approximation, i.e., $\sin θ \approx θ$ , we assume that the perturbed lateral motion from the initial wellbore configuration is small. Note that $θ_{1, 0} = ζ_{1} / l_{1}$ , and $θ_{2, 0} = χ_{1} / l_{1}$ . Using this in Equation (50), and performing the partial derivative in Equation (48), the generalized lateral restoring forces due to bending are derived as

\begin{matrix} τ_{b} = K q (t) \\ K = - Blockdiag (0, \bar{K} (k_{i}^{b}), \bar{K} (k_{i}^{b}), 0) \end{matrix}

(51)

where $0 \in^{n_{q} \times n_{q}}$ , and $\bar{K} \in^{n_{q} \times n_{q}}$ is the linear stiffness matrix in bending. Moreover, Equation (51) constitutes generalized forces due to linearized bending moments.

The gravitational forces are included assuming equal fluid densities inside the drill string and annulus. Hence, the buoyancy factor is expressed according to Aadnøy and Andersen²³ as

β = (1 - ρ_{m} / ρ_{s})

(52)

where the gravitational forces acting on $i$ in inertial vertical direction can be expressed as

f_{g, i}^{0} = m_{i} β g a_{3}

(53)

where $a_{3}$ is the fixed base frame unit vector pointing strictly downward, hence, the external forces from the acceleration of gravity are projected on $i$ by $R_{i}^{0}$ computed at $s_{i}$ .

The traveling block adjusts the vertical position of the drill string in the drilling rig. The control input is the block velocity, used to manipulate tripping speed. A change in block position results in a reaction force at $s = 0$ on $i = 1$ , given as

f_{h, 1}^{1} = {[\begin{matrix} \frac{E A_{1}}{l_{1}} u_{h} & 0 & 0 \end{matrix}]}^{⊤}

(54)

where $u_{h}$ is the block position.

The applied torque from the top drive connection shaft yields the torque input to the first element of the drill string. From the shaft rotation angle, given by the drive system, the input torque is given as

t_{td, 1}^{1} = {[\begin{matrix} \frac{G J_{1}}{l_{1}} ϕ_{s} (t) & 0 & 0 \end{matrix}]}^{⊤}

(55)

which is applied at $i = 1$ , and $ϕ_{s} (t)$ is the rotation angle of the top drive connection shaft. The top-of-string torque at $s = 0$ can then be written as

T_{0} = \frac{G J_{1}}{l_{1}} (ϕ_{s} - ϕ_{1}) .

(56)

The applied forces from the block and drive are then

ℱ_{i}^{a} = (\begin{array}{l} {[f_{h, i}^{i^{⊤}} t_{td, i}^{i^{⊤}}]}^{⊤} & i = 1 \\ 0_{6 \times 1} & i \neq 1 \end{array} .

(57)

Drag forces due to pipe motion in the surrounding fluid are assumed to be proportional to the linear and angular velocities. The effects of damping are included by considering possible turbulent and laminar flow from viscous drag forces. We can express the forces and torques from fluid drag as ^24,25

\begin{matrix} f_{d, i}^{i} = - c_{1, i}^{f} | | v_{i b_{i}}^{i / c} | |^{2} \frac{v_{i b_{i}}^{i / c}}{| | v_{i b_{i}}^{i / c} | |} - c_{2, i}^{f} v_{i b_{i}}^{i / c} \\ = - c_{1, i}^{f} | | v_{i b_{i}}^{i / c} | | v_{i b_{i}}^{i / c} - c_{2, i}^{f} v_{i b_{i}}^{i / c} \end{matrix}

(58)

where $c_{1, i}^{f} = ρ_{m} C_{d} (A_{w} - A_{i}) / 2$ , $C_{d}$ is an estimated drag force coefficient, $A_{w} = π d_{w}^{2} / 4$ , $A_{i} = π d_{i, dp / dc}^{o^{2}} / 4$ , and $c_{2, i}^{f}$ is the coefficient related to the viscous part of the drag forces. Similarly, the estimated damping due to rotational motion is assumed to be dominated by the viscous part proportional to the angular velocity

t_{d, i}^{i} = - c_{i}^{r} ω_{i b_{i}}^{i}

(59)

and $c_{i}^{r}$ is an estimated coefficient. The viscous damping forces and torques can be rewritten using Equations (58), (59) and (36), yielding

\begin{matrix} F_{i}^{d} = - D_{i} ν_{i b_{i}}^{i / c} \\ D_{i} = [\begin{matrix} I (c_{1, i}^{f} | | ν_{i b_{i}}^{i / c} | | + c_{2, i}^{f}) & 0 \\ 0 & I c_{i}^{r} \end{matrix}] \end{matrix}

(60)

where $D_{i} \in^{6 \times 6}$ .

4.1. Wellbore contact and friction

In this section, we define the interaction forces due to the transient motion of the drill string in the wellbore. A pipe cross-section is presented in Figure 6 in the normal and bi-normal plane, i.e., the local $^{i} {\vec{c}}_{2} -^{i} {\vec{c}}_{3}$ plane.

Figure 6.

Drill pipe or collar in contact with wellbore wall. Note that all the forces and torques are lumped to the $i th$ element. The green dot–dashed curve represents the wellbore curve $C$ .

The impulse response from drill string contact with the wellbore formation can be expressed by a penalty formulation,^14,24 being derived next.

The element $b_{i}$ centroid position is given by the normal and binormal component of ${\vec{p}}_{i b_{i}}^{c}$ in Equation (30). The angle $ψ_{i} = \arctan ({\vec{p}}_{i b_{i}, [2]}^{c} / {\vec{p}}_{i b_{i}, [1]}^{c}) \neq ϕ_{i}$ describes the planar whirl of the element, where subscripts 1 and 2 refer to the first and second vector components. At the moment of contact, $| | {\vec{p}}_{i b_{i}, [1, 2]}^{c} | | - r_{0} = 0$ where $r_{0} = (d_{w} - d_{i, dp / dc}^{o}) / 2$ is the nominal clearance for each element on the curve $C$ . The contact normal force arising from impact can be expressed as

{\vec{f}}_{n, i} = - k_{co} (1 - \frac{r_{0}}{δ_{c, i}}) {\vec{p}}_{i b_{i}, [1, 2]}^{c} - c_{co} | {\overset{\cdot}{δ}}_{c, i} | \frac{{\vec{v}}_{i b_{i}, [1, 2]}^{c}}{| | {\vec{v}}_{i b_{i}, [1, 2]}^{c} | |}

(61)

where $δ_{c, i} = | |_{i b_{i}, [1, 2]}^{c} | |$ , with $δ_{c, i}$ being the displaced length of the $i th$ center. The second term on the left-hand side of Equation (61) is the formation damping, assumed proportional to the magnitude of ${\overset{\cdot}{δ}}_{c, i}$ in the direction of the pipe segment velocity. Hence, the formation stiffness $k_{co}$ and damping $c_{co}$ are chosen by trial-and-error during the initial simulation phase, such that each element rests on the wellbore in steady state.

A friction force perpendicular to the contact normal force and in the opposite direction of the segment velocity is generated at the point of contact. The slip velocity between the two surfaces (pipe and wellbore) in contact (as sketched in Figure 6) can be expressed as²⁶

{\vec{v}}_{i b_{i}}^{s} = {\vec{v}}_{i b_{i}}^{c} + {\vec{ω}}_{i b_{i}} \times (\frac{{\vec{p}}_{i b_{i}, [1, 2]}^{c}}{| | {\vec{p}}_{i b_{i}, [1, 2]}^{c} | |} \frac{d_{i, dp / dc}^{o}}{2})

(62)

where ${\vec{v}}_{i b_{i}}^{s}$ is the slip velocity. The friction force arising from wellbore contact is dependent on the unit vector of the slip velocity, given as

{\vec{f}}_{f, i} = - μ_{k} | | {\vec{f}}_{n, i} | | \frac{{\vec{v}}_{i b_{i}}^{s}}{| | {\vec{v}}_{i b_{i}}^{s} | |}

(63)

where $μ_{k}$ is the coefficient of kinetic friction (also referred to as the Coulomb friction coefficient), and $| | {\vec{f}}_{f, i} | | \leq μ_{s} | | {\vec{f}}_{n, i} | |$ , where $μ_{s} > μ_{k} > 0$ . Hence, when the slip velocity is zero, the friction force balances the externally applied forces of $i$ .

The frictional torque due to contact is given by the magnitude of the friction force, expressed as

t_{f, i}^{i} = - \frac{d_{i, dp / dc}^{o}}{2} | | f_{f, i}^{i} | | \frac{ω_{i b_{i}}^{i}}{| | ω_{i b_{i}}^{i} | |} - c_{f, i}^{c} ω_{i b_{i}}^{i}

(64)

where the rightmost term denotes the viscous torque proportional to element angular velocity.

The contact forces and torques are summarized as

F_{i}^{c} = {[{(f_{n, i}^{i} + f_{f, i}^{i})}^{⊤} {t_{f, i}^{i}}^{⊤}]}^{⊤}

(65)

4.2. Equations of motion

The equations of motion are formed using Equation (26) with Equations (40) and (43), and together with the defined generalized active and restoring forces and torques, given as

\sum_{i}^{n} J_{i}^{⊤} [\begin{matrix} m_{i} a_{i b_{i}}^{i / c} \\ I_{b_{i}}^{i / c} α_{i b_{i}}^{i} \end{matrix}] = \sum_{i}^{n} J_{i}^{⊤} (F_{i}^{g} + F_{i}^{e} + F_{i}^{d} + F_{i}^{a} + F_{i}^{c})

which can be rewritten using Equations (39) and (60) as

\sum_{i}^{n} J_{i}^{⊤} (M_{i} {\overset{\cdot}{ν}}_{i b_{i}}^{i / c} + D_{i} ν_{i b_{i}}^{i / c}) = \sum_{i}^{n} J_{i}^{⊤} (F_{i}^{g} + F_{i}^{e} + F_{i}^{a} + F_{i}^{c})

where

M_{i} (q) = [\begin{matrix} (m_{i} + m_{f}) I, & 0 \\ 0 & I_{i}^{i / c} \end{matrix}] \in R^{6 \times 6} .

Adding in the forces resulting from linear bending in Equation (51), and using the generalized velocities in Equation (39), we get

\begin{matrix} \sum_{i}^{n} (J_{i}^{⊤} M_{i} J_{i} \overset{\cdot\cdot}{q} + J_{i}^{⊤} M_{i} {\overset{\cdot}{J}}_{i} \overset{\cdot}{q} + J_{i}^{⊤} D_{i} J_{i} \overset{\cdot}{q}) \\ + K q = \sum_{i}^{n} J_{i}^{⊤} (F_{i}^{g} + F_{i}^{e} + F_{i}^{a} + F_{i}^{c}) \end{matrix}

(66)

and hence, the equations of motion of the drill string system can be written as second-order ordinary differential equations including the linear bending stiffness, yielding

M (q) \overset{\cdot\cdot}{q} + C (q, \overset{\cdot}{q}) \overset{\cdot}{q} + D (q, \overset{\cdot}{q}) \overset{\cdot}{q} + K q = τ (q, \overset{\cdot}{q}, t)

(67)

where $M = \sum_{i}^{n} J_{i}^{⊤} M_{i} J_{i} \in R^{n_{q} \times n_{q}}$ is the system inertia matrix, $C (q, \overset{\cdot}{q}) = \sum_{i}^{n} J_{i}^{⊤} M_{i} {\overset{\cdot}{J}}_{i} \in R^{n_{q} \times n_{q}}$ is the Coriolis matrix, $D (\overset{\cdot}{q}) = \sum_{i}^{n} J_{i}^{⊤} D (\overset{\cdot}{q})_{i} J_{i} \in R^{n_{q} \times n_{q}}$ is the system damping matrix, and $τ = τ_{a} + τ_{g} + τ_{e} + τ_{c}$ are the generalized forces and torques defined as

τ_{a} = \sum_{i}^{n} J_{i}^{⊤} F_{i}^{a}, τ_{g} = \sum_{i}^{n} J_{i}^{⊤} F_{i}^{g},

(68)

τ_{c} = \sum_{i}^{n} J_{i}^{⊤} F_{i}^{c}, τ_{e} = \sum_{i}^{n} J_{i}^{⊤} F_{i}^{e},

(69)

where $τ_{a}$ is the applied forces and torques, $τ_{g}$ is the external gravitational forces, $τ_{c}$ is the generalized contact forces and torques, and $τ_{e}$ represents the generalized tension forces.

5. Steady-state configuration

The LM configuration imposed by the elastic deformation described by Equation (67) is in this section compared to a similar FEM model in Abaqus, an acknowledged commercial software for finite element analysis (FEA). For linear analysis, the B21 element is used, which is a 2-node linear beam in a plane. For nonlinear analysis, the B22 element is used, which is a 3-node quadratic beam in a plane. The analysis is performed to validate the elastic force contributions in the LM, which represents bending, given by the lateral displacement of each element.

For a deviated wellbore configuration the drill string model in Equation (67) is deformed by gravity in both the lateral and longitudinal directions, where both the axial stiffness and bending stiffness of the pipe affect the configuration. To verify the steady-state solution of the linear bending forces, the two effects are isolated in separate tests: a vertical and a horizontal configuration.

The drill string is stretched due to the structural weight, and under the influence of gravity and buoyancy with mud present in the annulus. We assume that the wellbore is arbitrarily large, i.e., $F_{i}^{c} = 0$ in this section. When the drill string motion reaches steady state, we have $\overset{\cdot\cdot}{q} (t) = \overset{\cdot}{q} (t) = 0$ and Equation (67) reduces to

K q_{eq} + τ_{a} (q_{eq}) + τ_{e} (q_{eq}) + τ_{g} = 0

(70)

where $q_{eq}$ is the equilibrium position when the drill string is subject to externally applied forces and torques, along with the tension forces and gravitational forces.

5.1. Free-hanging vertical solution

The vertical hanging static case is compared with the FEA results, and we assume that the drill string only subject to its dry-weight, hence $β = 0$ . The Case 1 parameters are presented in Table 1.

Table 1.

Case 1 parameters.

	Description	Value	Unit
$E$	Young’s modulus	200 $\cdot$ 10 $^{9}$	Pa
$G$	Shear modulus	79 $\cdot$ 10 $^{9}$	Pa
$ρ_{s}$	Density, steel	7850	kg $\cdot$ m $^{- 3}$
$L_{dc}$	Length BHA	120	m
$L_{vr}$	Vertical reach	1000	m
$L_{hr}$	Horizontal reach	0	m
$d_{i, dp}^{o}$	Drill pipe OD	0.127	m
$d_{i, dp}^{i}$	Drill pipe ID	0.1086	m
$d_{i, dc}^{o}$	Drill collar OD	0.165	m
$d_{i, dc}^{i}$	Drill collar ID	0.0572	m
$I_{dp}^{x}$	$2^{nd}$ moment of inertia	0.0594 $\cdot$ 10 $^{- 4}$	m $^{4}$
$I_{dc}^{x}$	$2^{nd}$ moment of inertia	0.3595 $\cdot$ 10 $^{- 4}$	m $^{4}$

OD: Outer diameter; ID: inner diameter.

The displacements at six locations from discretization of the vertical drill string configuration into 10, 20, and 100 lumped elements are compared with the converged FEA solution. The displacements are shown in Figure 7.

Figure 7.

Static configuration for Case 1, with converged FEA solutions for comparison.

The general observation in Figure 7 is that by increasing $n$ a better convergence to the distributed properties of the structured is obtained for the vertical configuration.

5.2. Convergence to a horizontal solution

A simply supported beam in horizontal orientation with a pinned-slider boundary condition is used to verify the forces due to bending and tension. The model representation is depicted in Figure 8.

Figure 8.

Static configuration for Case 2. Dotted line represents a sketch of the deformed configuration.

A static force, $f_{e}$ , applied at the bit $(s_{n})$ in the direction of $^{i} {\vec{c}}_{1}$ , is included to keep tension in the drill string to avoid a steep U-shape, keeping the drill bit fixed. For a linear analysis, the horizontal force does not affect the bending deformation. This means that a linear model configured with sliders yields the same shape as a cable pinned at both ends.

The LM is moment free at $i = 1$ , i.e., the first element agrees with the pinned boundary condition of the simply supported beam. For the slider boundary condition in the beam, springs are used in the LM to fix the origin position of element $i = n$ in the $^{n} {\vec{c}}_{2} -^{n} {\vec{c}}_{3}$ plane.

The supporting forces acting on element $b_{n}$ in Equation (67) are given as

f_{n}^{s} = {[\begin{matrix} f_{e} & - k_{1} ζ_{n} & - k_{2} χ_{n} \end{matrix}]}^{⊤}

(71)

where $k_{1}, k_{2}$ are lateral spring coefficients for the LM.

The Case 2 parameters different from Case 1 are given in Table 2. Three values of the acceleration of gravity constant $(g_{(1 - 3)})$ are tested to investigate the effects of pure bending and the nonlinear effects of bending and axial tension for large deformations. Note that $f_{e} = 0$ only for the pure bending analysis. The LM in this paper are discretized into 10 and 20 elements to provide insight into the lumped-parameter approximation of the distributed pipe properties.

Table 2.

Case 2 parameters for $g_{(1 - 3)}$ , and the structural parameters are given in Table 1.

	Description	Value (1–3)	Unit
$L_{vr}$	V. reach	$0$ , $0$ , $0$	m
$L_{hr}$	H. reach	$10^{3}$ , $10^{3}$ , $10^{3}$	m
$k_{n 1}$	Support	$10^{9}$ , $10^{9}$ , $10^{9}$	$N \cdot m^{- 1}$
$k_{n 2}$	Support	$10^{9}$ , $10^{9}$ , $10^{9}$	$N \cdot m^{- 1}$
$g$	Acc. of gravity	$9.81 \cdot {10^{- 6}, 10^{- 4}, 1}$	$m \cdot s^{- 2}$
$f_{e}$	Static force	$0$ , $4$ *, $1500 \cdot 10^{3}$	N
$p_{\infty}$	$G_{α}$	$0.05$	–
$α_{m}$	$G_{α}$	$(2 p_{\infty} - 1) / (p_{\infty} + 1)$	–
$α_{f}$	$G_{α}$	$p_{\infty} / (p_{\infty} + 1)$	–

LM: lumped model; FEA: finite element analysis.

Applies to LM only, to replicate the nonlinear FEA model. Common $G_{α}$ parameters (shown at the bottom in this table) for all configurations.

The Generalized- $α$ $(G_{α})$ scheme (see, e.g., Chung and Hulbert²⁷) is implemented in the LM for nonlinear analysis. The latter method is an extension of the classical method by Newmark²⁸ for second-order systems, including an optimal combination of high accuracy at lower frequencies and damping for high-frequency modes. This method enforces dynamic equilibrium at every time step. The algorithm used in this work is drawn from Arnold and Brüls.²⁹ Furthermore, the parameter $p_{\infty} \in [0, 1]$ in Table 2 determines the amount of high-frequency numerical damping present in the integration, where a low value yields higher dissipation of the high-frequency response.

For pipe weights with $g_{1}$ and $g_{2}$ , the static configurations are shown in Figure 9. In these configurations, the Newton–Raphson method is used to solve Equation (70).

Figure 9.

Drill string curvature due to internal weight with $g_{1} = g \cdot 10^{- 6}$ (left plot) and $g_{2} = g \cdot 10^{- 4}$ . The black dots mark the nodes in the LM. FEA-L: linear; FEA-NL: nonlinear.

In the leftmost plot, for the loading conditions with the lowest acceleration of gravity and zero horizontal force, the linear and nonlinear converged FEA together with the LM shows similar results. The FEA solution is discretized into 20 m long elements, leading to a total of 50 elements. For this configuration the bending stiffness governs the shape of the drill string.

When the gravity is increased, as seen in the rightmost plot in Figure 9, the linear analysis will start to deviate from the nonlinear analysis. The nonlinear analysis for both FEA and LM also shows the contraction of the pipe due to coupling to the axial tension, with a negative horizontal displacement.

When $g_{3} = g$ , the FEA and LM are subject to a static horizontal force $f_{e} = 1500 kN$ , which is the same as a pinned boundary condition at the bit end. Furthermore, $LM$ is implemented with 10 and 20 uniformly distributed elements to compare the effect of increasing the number of lumped elements. The results are shown in Figure 10.

Figure 10.

Drill string curvature due to internal weight with $g_{3} = 9.81 m s^{- 2}$ . The black dots mark the nodes in the LM.

A close match with the FEA is shown both for the pure bending case and for the more involved case of large deformations. Since many of the significant effects such as complex geometrical stiffening effects are neglected in the LM, deviations are evident.

The elasticity of the LM increases the time for its structural configuration to reach steady state. In addition, the displacement is only an effect of the imposed external and restoring forces. Hence, our model converges to the free-hanging physical configuration of the drill string when the number of elements increase. Similar conclusions were drawn for a simplified discrete rod analysis in Dreyer and Vuuren,³⁰ however, without axial contraction between the model elements.

6. Simulation and results

In this section, a simulation study is performed, and the qualitative properties of the model are presented. The numerical solution forward in time $t$ is obtained using the $G_{α}$ integrator from Section 5. Note that the tangent damping and stiffness matrices are derived from Equation (67), preliminary to the simulation procedure from the model configuration. A J-well configuration with ${\vec{p}}_{0 i} (s_{i}) = L_{hr}$ $s_{i} {\vec{a}}_{3}$ is used in this section, and the wellbore and drill string parameters are given in Table 3.

Table 3.

Model parameters with $n = 12$ for the J-well studies.

	Description	Value	Unit
$ρ_{m}$	Density, mud	1840	kg $\cdot$ $m^{- 3}$
$L_{dc}$	Length BHA	240	m
$L_{vr}$	Vertical reach	3200	m
$L_{hr}$	Horizontal reach	1000	m
$d_{w}$	Borehole diameter	0.2032	m
$ϵ_{1}, ϵ_{2}$	BHA eccentricity*	0, $10^{- 4}$	m
$μ_{k}$	Kinetic friction factor	0.2	–
$μ_{s}$	Static friction factor	0.3	–
$k_{co}$	Formation stiffness	$10^{9}$ $(d_{i, \cdot}^{o} / d_{i, dc}^{o})$	N $\cdot$ $m^{- 1}$
$c_{co}$	Formation damping	$5 \cdot 10^{6}$ $(d_{i, \cdot}^{o} / d_{i, dc}^{o})$	N $\cdot$ s $\cdot$ $m^{- 1}$
$v_{tr}$	Slip velocity threshold	$10^{- 5}$	m $\cdot$ $s^{- 1}$
$C_{a}$	Added mass coeff.	1	–
$C_{d}$	Drag force coeff.	1	–
$c_{i}^{r}$	Damping coeff.	20	N $\cdot$ m $\cdot$ s $\cdot$ ${rad}^{- 1}$
$c_{2, i}^{f}$	Damping coeff.	10	N $\cdot$ s $\cdot$ $m^{- 1}$
$c_{f, i}^{c}$	Viscous friction coeff.	1.5	N $\cdot$ m $\cdot$ $s$ ${rad}^{- 1}$
$p_{\infty}$	$G_{α}$	0.1	–
$h$	Simulation time-step	0.001	s

BHA: bottom-hole assembly.

The structural parameters are given in Table 1.

Applies for the model in rotation and lifting.

The friction forces are modeled by Coloumb friction in Equation (63). These forces balance the inertia and external forces on $i$ when the slip velocity is zero. From Equation (66), the equations of motion for equilibrium of mass $i$ are expressed as

K q + J_{i}^{⊤} (ℳ_{i} J_{i} \overset{..}{q} + ℳ_{i} {\dot{J}}_{i} \dot{q} + D_{i} J_{i} \dot{q} - ℱ_{i}) = 0

where $F_{i}$ is the sum of all forces on the right-hand side of Equation (66). Then, the equilibrium forces when the slip velocity is zero is expressed as

ℱ_{i}^{s} = ℳ_{i} (J_{i} \overset{..}{q} + {\dot{J}}_{i} \dot{q}) + D_{i} J_{i} \dot{q} -

(72)

ℱ_{i}^{g} - ℱ^{e} - ℱ_{i}^{a} - {[\begin{matrix} {ℱ_{n, i}^{i}}^{⊤} & 0 \end{matrix}]}^{⊤}

(73)

where $F_{i}^{s}$ is the equilibrium force required to keep the mass $i$ at rest during zero slip, and we have assumed negligible bending stiffness contribution due to large lengths, and that $J_{i}^{⊤} \neq 0$ . Hence, for the $i th$ element to break away from stiction, the sum of forces need to be larger than zero for the element to start accelerating. For computational purposes, we implement the approach by Karnopp³¹ with a slip velocity dead zone given as $v_{tr}$ .

6.1. Block and torque control

We assume that we can manipulate the traveling block velocity in terms of controlling drill string position in the wellbore. The block is subject to heave and velocity control from the driller, and its position is defined by the block position input $u_{h}$ in Equation (54). When the drill string is stationary in the rig floor slips, $u_{h}$ follows the rig floor elevation imposed by sea heave motion.

A variable speed drive (VSD) controls the torque supplied to the drill string. For the purposes of including a dynamic top drive model, the motor is given as

J_{m} {\overset{\cdot\cdot}{ϕ}}_{m} = u - b_{m} {\overset{\cdot}{ϕ}}_{m} - T_{0} n_{g}

(74)

where $J_{m}$ is the motor plus gear inertia, $u$ is the top drive input, i.e., the drive shaft torque from the VSD, $b_{m}$ represents the motor losses being proportional to motor speed, and $n_{g} < 1$ is the reduction gear ratio. Note that $ϕ_{s} = \int^{t} n_{g} {\overset{\cdot}{ϕ}}_{m} d τ$ used in Equation (55), as the gear reduction reduces speed, and increases the pipe torque.

Using a proportional-integral (PI) controller, the input torque supplied by the VSD is computed from the desired top-of-string angular velocity, yielding

u = k_{p} ({\overset{\cdot}{ϕ}}_{d} / n_{g} - {\overset{\cdot}{ϕ}}_{m}) + k_{i} \int^{t} ({\overset{\cdot}{ϕ}}_{d} / n_{g} - {\overset{\cdot}{ϕ}}_{m}) d τ

(75)

where ${\overset{\cdot}{ϕ}}_{d}$ is the desired angular velocity at the top-of-string, and $k_{p}, k_{i}$ are the controller parameters. Note that standard anti-windup of the integrator terms in Equation (75) by back-calculation is implemented. A high-gain controller with parameters drawn from Aarsnes et al.³² is implemented, and the drive parameters are shown in Table 4.

Table 4.

Top drive parameters.

	Description	Value	Unit
$n_{g}$	Gearbox ratio	1/7	–
$J_{m}$	Top drive inertia	900	$kg \cdot m^{2}$
$b_{m}$	Viscous friction in rotor	1.0	$N \cdot m \cdot s \cdot ra d^{- 1}$
$k_{p}$	Proportional gain, VSD	$100 J_{1} \sqrt{G ρ_{s}}$	$N \cdot m \cdot s \cdot ra d^{- 1}$
$k_{i}$	Integral gain, VSD	$5 J_{m}$	$N \cdot m \cdot s \cdot ra d^{- 1}$

VSD: variable speed drive.

6.2. Surge and swab pressures

The pressure profile in the wellbore is an important parameter for monitoring and maintaining well-stability, by safeguarding the drill string velocity. In this section, we include a static model for computing the wellbore annular pressures. Hence, no coupling is included such as the effect of damping of drill string vibrations due to mud flow (see, e.g., Wilson and Noynaert⁴).

During drill string movement, the pressure profile in the wellbore is subject to surge (running the drill string down in the well) and swab (retrieving the drill string) pressures. To avoid large surge and swab pressures causing instabilities in the well, the tripping velocity of the block is manipulated.

A symmetrical cross-section of the element was assumed, and we neglect the mud flow variations due to tool joints. The average mud flow velocity can be calculated assuming closed pipe ends, when all the fluid is displaced upward in the annular section.³³ The mud velocity is given by the pipe velocity of the $i th$ element as

{\bar{v}}_{m, i} = - ({\vec{v}}_{i b_{i}}^{c} \cdot^{i} {\vec{c}}_{1}) \frac{d_{i, dp / dc}^{o^{2}}}{d_{w}^{2} - d_{i, dp / dc}^{o^{2}}}

(76)

where ${\bar{v}}_{m, i}$ is the average annular mud flow velocity assumed to be defined by the axial velocity component of the drill string, and the equation yields the case of no mud being pumped into the drill string from top-side.

The pressure loss over the drill string segment is a function of mud flow velocity and its flow regime. The Darcy–Weisbach equation can be used as an approximation to obtain the pressure loss due to pipe movement.³⁴ It is defined as

Δ P_{f, i} = f_{D} \frac{l_{i}}{D_{h, i}} ρ_{m} \frac{{\bar{v}}_{m, i}^{2}}{2}

(77)

where $Δ P_{f, i}$ is the pressure loss, $f_{D}$ is the Darcy friction factor, $l_{i}$ is the section length, and $D_{h, i}$ is the hydraulic diameter.

In annular flow, where the inner section moves and the outer section is stationary, it can be useful to define the effective velocity,^17,33 given as

v_{e, i} = {\bar{v}}_{m, i} + γ_{c} ({\vec{v}}_{i b_{i}}^{c} \cdot {\vec{c}}_{1})

(78)

= - (\frac{d_{i, dp / dc}^{o^{2}}}{d_{w}^{2} - d_{i, dp / dc}^{o^{2}}} + γ_{c}) ({\vec{v}}_{i b_{i}}^{c} \cdot {\vec{c}}_{1})

(79)

where $γ_{c}$ is denoted as the clinging factor which characterizes the portion of fluid following the pipe wall, being calculated on the basis of either laminar or turbulent flow.

The laminar and turbulent flow regimes are used to characterize the Darcy friction factor $f_{D}$ , in order to compute Equation (77). Instead, we express the pressure loss in terms of the wall shear stresses. Hence, the steady-state pressure loss due to pipe movement in the annulus for element $i$ is obtained by using $f_{D, i} = 8 τ_{w, i} / (ρ_{m} {\bar{v}}_{m, i}^{2})$ , yielding

Δ P_{f, i} (v_{e, i}) = \frac{4 l_{i}}{d_{w, i} - d_{i, dp / dc}^{o}} τ_{w, i} (v_{e, i})

(80)

where $d_{w} - d_{o, i} = D_{h, i}$ , and $τ_{w, i}$ is the wall shear stresses, which in this paper is computed from the work in Gjerstad et al.¹⁷ (section 4.1.1; Appendix I).

An estimate of the bottom hole pressure (BHP) is then obtained as

P_{bhp} = ρ_{m} g (TVD + \cos (η_{n} (s_{n})) l_{n}) + Δ P_{f, n}

(81)

where $η_{n} (s_{n}) = \arccos (c_{1, 1} / | | {\vec{c}}_{1} | |)$ is the inclination at $i$ along $C$ , $c_{1, 1}$ is the first component of ${\vec{c}}_{1}$ , the first term on the left-hand side is the hydrostatic pressure measured at the true vertical depth (TVD) measured along ${\vec{a}}_{3}$ to the position of the element, and the second term is the pressure loss due to inclination.

6.3. Heave-induced pressure oscillations

In this section, the drill string is fixed in the drill floor slips during pipe connection. Hence, the drill string is subject to heave motion from the rig.

A simple wave characteristic is assumed, where the heave amplitude and wave period are the only parameters considered.³⁵ The rig floor position due to heave is given as

u_{heave} = Z \sin (2 π t / T)

(82)

where $Z$ is the wave amplitude, $T$ is the wave period, and $t$ is the time.

At time $t = 0$ , the drill string initial configuration is subject to gravity and contact forces, resting on the wellbore wall. The wellbore configuration and element position are seen in Figure 11.

Figure 11.

Wellbore configuration with model element position indicated by circles.

The drill string is being suspended in the rig-floor slips at $t = 10$ being subject to the heave motion given by Equation (82). The heave trajectory is shown in the topmost plot in Figure 12. The BHP (seen in the lower-most plot) fluctuates due to transient surge and swab pressure at the BHA given by $Δ P_{f, n}$ . The maximum and minimum BHP are generated for the 3 m amplitude waves with approximately 622 and 592 bar, respectively. The difference between the maximum and minimum BHP values is close to what is reported in Zhao et al.¹² for axial velocities in the vicinity of 2 m/s.

Figure 12.

The effect of heave due to fixation in rig floor slips. Solid line represent wave height of 2 m, and dashed line for wave height of 6 m.

From plots two and three (counting top-down) in Figure 12, the clearance variable $δ_{c, n}$ and the corresponding norm of the contact force are shown. During heave, the contact forces in Equation (61) tend to vary with the minor changes in lateral position and velocity.

As observed from Figure 12, in plot four and five the bit sticks for a moment during the wave top. This correlates with plot five showing the slip velocity, where axial stiction is obviously occurring for lower heave amplitudes.

6.4. Rotation while lifting

To evaluate the consistency of the model toward coupled rotation and translatory motion, we utilize the same wellbore and drill string configuration as in Section 6.3. A pipe-stand trip sequence is performed with and without rotation, to illustrate some transient characteristics related to axial friction while rotating.

Tripping out one stand and running it back in is seen in the first 90 s of Figure 13(a) and (b) for $μ \in {0.2, 0.3}$ . Note that BHA eccentricity of $ϵ_{2} = 0.1 mm$ is included. The corresponding drive outputs are presented in Figure 14(a) and (b).

Figure 13.

Rotating while running out and back in one stand. For solid lines the BHA elements has $ϵ_{1} = ϵ_{2} = 0$ , and for dashed lines the BHA has $ϵ_{1} = 0, ϵ_{2} = 10^{- 4}$ ; the rest of the drill string elements always have $ϵ_{1} = ϵ_{2} = 0$ : (a) $μ_{k} = 0.2$ and (b) $μ_{k} = 0.3$ .

Figure 14.

Down-hole bit angular velocity and motor and top-of-string torque while running out and in the drill string one stand. For solid lines the BHA elements has $ϵ_{1} = ϵ_{2} = 0$ , and for dashed lines the BHA has $ϵ_{1} = 0, ϵ_{2} = 10^{- 4}$ ; the rest of the drill string elements always have $ϵ_{1} = ϵ_{2} = 0$ : (a) $μ_{k} = 0.2$ and (b) $μ_{k} = 0.3$ .

The tripping velocity is $0.3 m s^{- 1}$ during the sequences, and we note that the bit experiences minor oscillations during axial movement (topmost plots of Figure 13(a) and (b)). Hence, when starting the top-drive and applying torque to the top of the drill string, a smoother trajectory for the bit axial motion is seen.

The intention of rotating while pulling is to reduce axial drag and friction between the drill string and the wellbore.³⁶ We observe from the second plot in Figure 13(a) and (b), that the tension difference in the uppermost segment of the drill string reduces by around 5 tons when rotating and tripping. The lowermost plot shows the slip velocity at the bit, seen to be dampened out due to increase in the angular velocity of the bit.

Increasing the friction coefficients in the model tends to lead to a slight increase in the uneven bit axial motion (topmost plot of Figure 13(b)) during tripping without rotation, and a general amplification of the upper segment tension and slip velocity. The bit angular velocity also experiences an increased oscillatory pattern, as seen in the topmost plot in Figure 14(b).

7. Discussion

In this work, no model for axial bit–rock interaction and coupled torque-on-bit has been included. Improvements in future work involve extending the model applicability to downhole conditions. Furthermore, comparison toward similar models and field data are steps to further establish the validity of the model. The authors do not presently have access to field data for comparison. As any well described by a parametric curve can be implemented, a realistic well from work such as Aadnøy and Andersen²³ could facilitate further model validation. The rest of this section includes a discussion on model effects and a study of two numerical integration methods.

7.1. Geometrical stiffening effects

The drill string is confined in the borehole such that limited lateral movement is allowed, and the curvature at $s_{i}$ can be assumed to be low. Hence, the ratio of the length-to-diameter of the element is large. Geometrical stiffening effects can increase the ability of the model to represent the actual stress and strain (by means of using Green’s strain) of a drill string segment. However, this will lead to coupling terms in Equation (67), and limit the efficiency of the implementation presented in this paper. Moreover, existing FEM could be better suited if exact replication of the structural properties is desired.

7.2. Selection of numerical solver

The numerical stability dictates the performance during the simulation and accuracy of the results obtained. For second-order models used to represent drill string dynamics, the method by Newmark²⁸ (see, e.g., Wilson and Noynaert,⁴ Butlin and Langley³⁷), Runge–Kutta, and the forward Euler methods^12,38 are among some of the more common integration schemes found in the literature. Nonlinear models comprising contact dynamics and stick-slip behavior may require a relatively short time step $(h \approx 10^{- 4})$ for convergence. As such, efficient numerical solvers that do not introduce higher frequencies over time is important, since increasing the physical damping in the model to overcome solver deficiencies is generally not a good approach.

The stability region for a Runge–Kutta method of order 4 in terms of an undamped oscillator can be drawn from Hairer and Wanner,³⁹ Figure 2.1, yielding

h \leq 2.8 / ω_{\max}

(83)

where $ω_{\max}$ is the largest eigenfrequency present in the system. The latter then dictates the critical time-step size for a stable RK4 solver.

In our model, the contact forces generated from drill string element impact with the formation is given by the penalty model in Equation (61). In free hanging form, boundary conditions alter the eigenfrequencies of the drill string in each of the degrees-of-freedom. However, when contact occurs, the formation stiffness is typically larger than the structural stiffness. For the Generalized- $α$ method, the tangent stiffness and damping matrix are formed from the partial derivative of $τ (q, \overset{\cdot}{q}, t)$ , with respect to $q$ and $\overset{\cdot}{q}$ , respectively. Hence, these matrices can in the case of a low number of lumped elements become ill-conditioned due to the large contact stiffness.

The time step for a simulation sequence has to be set according to the relevant frequencies in the system. Typically, when including a penalty-based model for the wall impact, the formation stiffness dictates the magnitude of the system eigenfrequencies. The eigenfrequency of the drill string when in contact with the wellbore can be given as

ω_{co, \max} = \sqrt{\frac{k_{co}}{min ([m_{1}, \dots, m_{n}])}}

(84)

where $\min (\cdot)$ denotes the minimum of all the individual drill string element masses.

For an RK4 method the maximum time-step size for a stable RK4 solver can be predicted for the drill string. For the simulated cases with $n = 12$ the mass of a drill pipe segment (from the upper block matrix in $M_{i}$ ) yields 19,540 kg. With a stiffness coefficient of $10^{9}$ , we get $ω_{co, \max} \approx 226 rad \cdot s^{- 1}$ , and $h \approx 0.0124$ . The time step defining RK4 numerical stability in equation (83) is then dependable for keeping a lower discretization of the drill string. Increasing $n$ requires a lower value of $h$ . As $m_{i}$ becomes sufficiently small, $h \approx 2.8 / \sqrt{k_{co}}$ . Hence, the Generalized- $α$ algorithm is better suited for shorter drill string element configurations for Equation (67), while simulating with larger time steps.

7.2.1. Brief comparison of solvers

This section gives a comparison of applying an RK4 and a $G_{α}$ for a periodic block position trajectory. As RK4 is among one of the most well-known numerical integration methods in the literature, we apply RK4 with $h = 10^{- 4}$ as the base case. Then, a measure of comparison can be given by the square of the error, computed according to

e_{n} = (| | {\vec{p}}_{n, b_{n}}^{RK 4_{1}} | | - | | {\vec{p}}_{n, b_{n}}^{S} | |)

(85)

where $S \in {G_{α 1}, \dots, G_{α 5}, RK 4_{1}, \dots, RK 4_{5}}$ , and the superscripts denote the solver configurations for the time steps $h \in {10^{- 4}, 10^{- 3}, 10^{- 2}, 2 \cdot 10^{- 2}, 10^{- 1}}$ .

The block position is given by a periodic function $u_{h} = \sin (2 π (t - 5) / 20 + (t - 5)^{2} / 400)$ to evaluate the error propagation and convergence. The bit position error between the RK4 and $G_{α}$ with Equation (85) using RK4 with $h = 10^{- 4}$ as the base case for comparison, is shown in Figure 15.

Figure 15.

Comparing $e_{n}^{2}$ for $h \in {10^{- 4}, 10^{- 3}, 10^{- 2}, 10^{- 1}}$ with $G_{α}$ and for $RK 4_{1}$ as base case in heave motion. For the lowermost plot the y-axis is in log scale, and solid line denotes $h = 10^{- 4}$ , dashed line $h = 10^{- 3}$ , dotted line $h = 10^{- 2}$ , and dash–dotted line $h = 10^{- 1}$ , respectively.

The square of the input to the block, along with the norm of the bit velocity for $RK 4_{1}$ and $G_{α 1}$ are shown in the first plot of Figure 15. It is noted a difference between the solvers between $t = 10$ and $20 s$ , where in general, the stiction phase is more stable for the $G_{α 1}$ below the threshold velocity. A general observation from the second plot in Figure 15, and the data in Table 5, is that we obtain root mean square error (RMSE) of the solvers in the range of centimeters for time steps beyond $10^{- 3}$ . For longer wellbores and drill string configurations, it could be expected that these errors increase, such that displacement errors lead to less physically accurate behavior downhole.

Table 5.

Comparing RTFs for simulation studies in this paper.

$h$	Solver	RTF	RMSE
0.0001	RK4	23.8	0
	$G_{α}$	36	$3.65 \cdot 10^{- 2}$
0.001	RK4	2.29	$5.18 \cdot 10^{- 4}$
	$G_{α}$	19.9	$6.1 \cdot 10^{- 3}$
0.01	RK4	0.226	$3.43 \cdot 10^{- 3}$
	$G_{α}$	3.06	$9.03 \cdot 10^{- 3}$
0.02	$G_{α}$	2.06	$1.51 \cdot 10^{- 2}$
0.1	$G_{α}$	0.475	$6.45 \cdot 10^{- 2}$

Note that RK4 for $h = 0.02$ and $0.1 s$ is numerically unstable, and did not converge.

RTF: real-time factor; RMSE: root-mean-square of the error.

The RMSE is computed from

RMSE = \sqrt{(1 / N) \sum_{k = 1}^{N} e_{k}^{2}}

(86)

where $e_{k}$ is the position error at each numerical time step. The RMSE, and a real-time factor (RTF) defined here from the actual code runtime to the fixed simulation time of one case, given as $RTF = t_{actual} / t_{final}$ , are presented in Table 5. The $G_{α}$ solver uses $p_{\infty} = 0.1$ .

The $G_{α}$ includes an internal iteration to force dynamic equilibrium in each time step depending on a tolerance criteria and a set maximum number of iterations. Hence, from the lower time-step values in Table 5 the solver is limited in achieving RTFs in real-time range. However, both implementation and the use of adaptive time-step algorithms need to be evaluated further before flagging this solver as inefficient for real-time. On the contrary, by increasing the time step the strength of $G_{α}$ appears as this method is unconditionally stable for larger time steps than RK4.

8. Conclusion

From the work presented in this paper, we have proposed an LM for analyzing transient drill string motion in arbitrary wellbore configurations. The wellbore is represented as a parametric curve in space, and the equations of motion are derived using Kane’s method.

A static configuration test toward FEA acting as a benchmark was performed. The LM, subject to axial tension and linear bending stiffness, conforms with the FEA for completely straight configurations. Evidently, increasing the number of elements in the discrete model yields better agreement. The model axial friction forces and contact forces were analyzed in two simulation studies for a deviated J-well drill string. The first presented the behavior during pipe connection, and the second showed the effect of rotating while tripping the drill string. Well-known characteristics such as reduced axial tension due to rotation were confirmed in simulation. A PI torque controller was implemented, along with manipulation of the block velocity.

The model computational efficiency in terms of error propagation and simulation time is compared with a Generalized- $α$ solver and the well-known Runge–Kutta order 4 method. The comparison shows that the model is suitable for real-time implementation with RTFs lower than 1 for this prototype implementation with an RK4, and a Generalized- $α$ for larger time steps. Combining this with the real-time-oriented drill string model formulation using Kane’s method, it is concluded that the work showcases novel solver and model configurations for drilling analysis. As such, the model is suitable as a virtual sensor for downhole monitoring.

Footnotes

Appendix 1 Funding

The research presented in this paper has received funding from the Norwegian Research Council, SFI Offshore Mechatronics, Project No. 237896.

ORCID iD

Njål Kjærnes Tengesdal

Author biographies

Njål Kjærnes Tengesdal received the PhD degree from the Department of Mechanical and Industrial Engineering, Norwegian University of Science and Technology, in 2021. He is currently positioned as an automation engineer, working with control systems in the oil and gas industry.

Gaute Fotland received the PhD degree from the Department of Mechanical and Industrial Engineering, Norwegian University of Science and Technology, in 2022. He is currently positioned as a mechanical design engineer, working with well intervention in the oil and gas industry.

Christian Holden is an Associate Professor at the Norwegian University of Science and Technology NTNU, Dept. of Mechanical and Industrial Engineering, Trondheim, Norway.

Bjørn Haugen is an Associate Professor at the Norwegian University of Science and Technology NTNU, Dept. of Mechanical and Industrial Engineering, Trondheim, Norway.

References

Henriksen

Storegjerde

. Cost effective horizontal drilling in the troll field through use of state of the art technology and optimal operations. In: Proceedings of the SPE/IADC drilling conference, 1997, pp. 57–68, https://onepetro.org/SPEDC/proceedings-abstract/97DC/All-97DC/SPE-37577-MS/189505

Andresen

Hovda

Olsen

. Experience with drilling C-26A, a world record extended reach horizontal well in the Oseberg Field, North Sea. In: SPE annual technical conference and exhibition: SPE-30463, 1995, https://onepetro.org/SPEATCE/proceedings-abstract/95SPE/All-95SPE/SPE-30463-MS/57638

Aarrestad

Blikra

Torque and drag-two factors in extended-reach drilling. J Pet Technol 1994; 46: 800–803.

Wilson

Noynaert

. A new damping model for nonlinear drillstring dynamics: understanding the effects of rotation, eccentricity, and confined fluid flow and their impact on unconventional drillstring design. In: IADC/SPE drilling conference and exhibition, 2016, https://onepetro.org/SPEDC/proceedings-abstract/16DC/2-16DC/D021S017R006/207821

Belyaev

Eliseev

VV.

Flexible rod model for the rotation of a drill string in an arbitrary borehole. Acta Mech 2018; 229: 841–848.

Aarsnes

UJF

Shor

. Torsional vibrations with bit off bottom: modeling, characterization and field data validation. J Pet Sci Eng 2018; 163: 712–721.

Goicoechea

Buezas

Rosales

MB.

A nonlinear Cosserat rod model for drill-string dynamics in arbitrary borehole geometries with contact and friction. Int J Mech Sci 2019; 157–158: 98–110.

Wang

Hasan

Zhang

, et al. Study on the hydrate formation region in a wellbore during Subsea Mudlift Drilling. Simulation 2021; 97(2): 69–82.

Karnopp

Margolis

Rosenberg

RC.

System dynamics: modeling, simulation, and control of mechatronic systems. Chichester: John Wiley & Sons, 2012.

10.

Foster

Macmillan

. High speed telemetry on wired drill pipe, history, and impact on drilling process: OTC-28735-MS. In: Offshore technology conference, 2018, https://onepetro.org/OTCONF/proceedings-abstract/18OTC/3-18OTC/D031S032R006/179743

11.

Florence

Iversen

. Real-time models for drilling process automation: equations and applications. In: IADC/SPE drilling conference and exhibition, 2010, pp. 1218–1233, https://onepetro.org/SPEDC/proceedings-abstract/10DC/All-10DC/SPE-128958-MS/106004

12.

Zhao

Hovda

Sangesland

Abnormal down hole pressure variation by axial stick-slip of drillstring. J Pet Sci Eng 2016; 145: 194–204.

13.

Cayeux

. On the importance of boundary conditions for real-time transient drill-string mechanical estimations. In: IADC/SPE drilling conference and exhibition, 2018, https://onepetro.org/SPEDC/proceedings-abstract/18DC/2-18DC/D021S010R005/214167

14.

Xie

Huang

, et al. Nonlinear dynamics of lump mass model of drill-string in horizontal well. Int J Mech Sci 2020; 174: 105450.

15.

Qin

Zhou

, et al. Nonlinear dynamic characteristics of the drill-string for deep-water and ultra-deep water drilling. J Pet Sci Eng 2022; 209: 109905.

16.

Kane

Levinson

DA.

Formulation of equations of motion for complex spacecraft. J Guid Control Dyn 1980; 3(2): 99–112.

17.

Gjerstad

Time

Bjorkevoll

. A medium-order flow model for dynamic pressure surges in tripping operations. In: SPE/IADC drilling conference, 2013, https://onepetro.org/SPEDC/proceedings-abstract/13DC/All-13DC/SPE-163465-MS/176802

18.

Egeland

Gravdahl

JT.

Modeling and simulation for automatic control, vol. 76. Trondheim: Marine Cybernetics, 2002.

19.

Kane

Levinson

DA.

Dynamics, theory and applications. New York: McGraw-Hill, 1985.

20.

Egeland

Sagli

JR.

Coordination of motion in a spacecraft / manipulator system. Int J Robot Res 1993; 12: 366–379.

21.

Tengesdal

Hovda

Holden

A discussion on the decoupling assumption of axial and torsional dynamics in bit-rock models. J Pet Sci Eng 2021; 202: 108070.

22.

Zhao

Hovda

Sangesland

Whirl simulation of drill collar and estimation of cumulative fatigue damage on drill-collar connection. SPE J 2018; 23(2): 286–300.

23.

Aadnøy

Andersen

. Friction analysis for long-reach wells. In: IADC/SPE drilling conference, 1998, pp. 819–834, https://onepetro.org/SPEDC/proceedings-abstract/98DC/All-98DC/SPE-39391-MS/190854

24.

Jansen

DJ.

Nonlinear dynamics of oilwell drillstrings. PhD Thesis, Delft University of Technology, Delft, 1993.

25.

Fossen

TI.

Handbook of marine craft hydrodynamics and motion control. Chichester: John Wiley & Sons, 2011.

26.

Gulyayev

Shevchuk

Drill string bit whirl simulation with the use of frictional and nonholonomic models. J Vib Acoust 2016; 138: 011021.

27.

Chung

Hulbert

GM.

A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method. J Appl Mech 1993; 60: 371–375.

28.

Newmark

MN.

A method of computation for structural dynamics. J Eng Mech Div 1959; 85: 64–94.

29.

Arnold

Brüls

Convergence of the generalized-α scheme for constrained mechanical systems. Multibody Syst Dyn 2007; 18: 185–202.

30.

Dreyer

Vuuren

JHV

. A comparison between continuous and discrete modelling of cables with bending stiffness. Appl Math Model 1999; 23: 527–541.

31.

Karnopp

DC.

Computer simulation of stick-slip friction in mechanical dynamic systems. J Dyn Syst Meas Control 1985; 107: 100–103.

32.

Aarsnes

UJF

Meglio

Shor

. Benchmarking of industrial stick-slip mitigation controllers. In: 3rd IFAC workshop on automatic control in offshore oil and gas production, 2018, pp. 233–238, https://folk.ntnu.no/skoge/prost/proceedings/ifac-oogp18/OOGP18/0040.pdf

33.

Whittaker

Staff

Theory and application of drilling fluid hydraulics. Boston, MA: D. Reidel Publishing Company, 1985.

34.

White

FM.

Fluid mechanics (McGraw-Hill Series in Mechanical Engineering). 4th ed. New York: McGraw-Hill, 2011.

35.

Haaø

Vangen

Tyapin

, et al. The effect of friction in passive and active heave compensators of crown block mounted compensators. IFAC-PapersOnLine 2012; 45(8): 316–320.

36.

Aadnøy

Fazaelizadeh

Hareland

A 3D analytical model for wellbore friction. J Can Pet Technol 2010; 49: 25–36.

37.

Butlin

Langley

RS.

An efficient model of drillstring dynamics. J Sound Vib 2015; 356: 100–123.

38.

Germay

Denoël

Detournay

Multiple mode analysis of the self-excited vibrations of rotary drilling systems. J Sound Vib 2009; 325: 362–381.

39.

Hairer

Wanner

Solving ordinary differential equations II: stiff and differential-algebraic problems. 2nd ed. New York: Springer, 2002.