Component-Based Modeling and Simulation of Nonlinear Drill-String Dynamics

Tengesdal, Njål Kjærnes; Holden, Christian; Pedersen, Eilif

doi:10.1115/1.4052676

Abstract

In this article, we present a dynamic model for a generic drill-string. The model is developed with the intention for component-based simulation with coupling to external subsystems. The performance of the drill-string is vital in terms of efficient wellbore excavation for increased hydrocarbon extraction. Drill-string vibrations limit the performance of rotary drilling; the phenomenon is well known and still a subject of interest in academia and in industry. In this study, we have developed a nonlinear flexible drill-string model based on Lagrangian dynamics to simulate the performance during vibrations. The model incorporates dynamics governed by lateral bending, longitudinal motion, and torsional deformation. The elastic property of the string is modeled by the assumed mode method, representing the elastic deformation, with a finite set of modal coordinates. By developing a bond graph model from the equations of motion, we can ensure correct causality of the model toward interacting subsystems. The model is analyzed through extensive simulations in case studies, comparing the qualitative behavior of the model with state-of-the art models. The flexible drill-string model presented in this article can aid in developing system simulation case studies and parameter identification for offshore drilling operations.

1 Introduction

The drill-string is a vital component of a drilling rig system and crucial in establishing the wellbore for hydrocarbon extraction. A drill-string consists of individual pipes joined together by tool joints, stabilizers to center the pipe in the wellbore, measurement-while-drilling tools, and sensor wiring. This structure can be several kilometers long. A generic model of the drill-string assembly is shown in Fig. 1.

Fig. 1

View large Download slide

A sketch of a drill-string assembly confined in a borehole

The string is suspended in the derrick to the top drive, moving upward when tripping out of hole and downward during drilling and run-in of new pipe stands. The lower section of the string, consisting of drill collars and the bit, is commonly referred to as the bottom hole assembly (BHA). The drill-collar section is a set of heavier pipes to produce enough weight on bit (WOB) when drilling and intervals of stabilizers keeping the string centered in the borehole. The torque on the bit is a consequence of the applied torque at the top, WOB, and the torque from the formation being excavated due to cuttings (cuttings from the rock smashed by the bit due to rotation) and the friction from the rock itself [1].

The drill-string classifies as a flexible slender object [2]. In general, large vibrations are a problem, and in terms of drilling, there exist three distinct vibration modes:

Torsional vibrations, occurring along the string or at the bit where the torque in the drill-string must overcome the static torque from the formation to break free. This results in periodic oscillatory behavior referred to as “stick-slip” [3,4].
Axial, compression, or longitudinal vibrations, also referred to as “bit-bounce.” This occurs when the bit periodically comes loose from the formation [5].
Lateral vibrations are also called whirling motion. Along the length of the string, the whirling motion tends to cause an uneven pattern in the well trajectory and can potentially result in formation fracture in the wellbore [5].

For the drill-string, vibrations can result in extensive wear or failure of the bit and other components in the assembly [5]. The vibration modes have been subject to extensive research for both understanding and mitigation through new control strategies.

Modeling and simulation of drill-string dynamics is a large research field, and subsequently, a large number of models for analyzing vibrations exist. An important contribution in the research field of drill-string vibrations is the work of Kyllingstad and Halsey [3]. The research in Ref. [3] concerns the string torsional vibrations induced by stick-slip. The vibrations are modeled as a transmission line reflecting the torsional waves created by the stick-slip friction.

Commonly found in control applications are low-order ordinary differential equation (ODE) models, confined to vertical well applications. In Ref. [6], a lumped model approach was used, representing each string section as a mass-spring-damper object, for both axial and torsional modes. The focus in Ref. [6] was the analysis of the effect of a downhole bit-rock interaction model. These models are suitable for real-time applications but are confined to uncoupled vibrations.

Christoforou and Yigit [7] have developed a coupled second-order model with axial, torsional, and lateral dependency with the objective of mitigation of vibrations. The model was developed under the assumption of one-mode approximation, and the dependency for the coupled vibrations were created by the bit-rock friction model. The dynamics of the drill-string was a discretized two degree-of-freedom (DOFs) model, where the coupling arose from borehole wall and bit-rock interaction.

Tucker and Wang [4] modeled the drill-string dynamics as an elastic rod with Cosserat theory. This method implies sectioning the rod into finite elements with unit vectors, each representing bending or twist, shear stress, and axial motion along the rod length. The model then consists of a set of nonlinear partial differential equations (PDEs).

Bakhtiari-Nejad and Hosseinzadeh [8] modeled the drill-string dynamics by solving the PDE for axial and torsional motion. This is further lumped by using a Ritz-series expansion with assumed modes to represent the vibrations. The drill-string model interacts with the bit-rock model from Ref. [9].

Jansen [5] analytically derived model equations for the vibrations listed earlier for sections of the drill-string (collar, etc.) and also uses a finite element model (FEM) as foundation for analyzing the vibration stages of a drill-string. Jansen also presented a computer simulation program for a nonlinear drill-string model. A FEM is favorable in terms of close-to-real results in applications but often requires large computational power and significant analysis time [10].

Feng et al. [11] presented the theory of a FEM for a drill-string and a comprehensive derivation of all the interacting dynamics between the drill-string and the wellbore. This study included coupled axial, lateral, and torsional dynamics toward the bit-rock model and fluid in the bottom hole. Through experimental verification, the model was shown to give promising results.

The FEM theory, used in relation with the Euler–Bernoulli beam with absolute-nodal-coordinate-formulation, was presented in the work by Caijin et al. [12] to model the torque and drag dynamics of a drill-string in contact with the wellbore. Using this approach, each local drill-string element had 12 degrees-of-freedom.

In the view of component-based modeling, the framework of bond graph (BG) modeling plays a significant role. The structure creates a unified approach for multidomain systems and allows for graphical modeling. An overview of this tool is found in Ref. [13]. In the study by Pedersen and Polic [14], a bond graph model framework for rotordynamic applications was developed. The resemblance between this and the current work is evident, and as such, the current model development is a formal extension to this bond graph model framework.

In terms of bond graph applications for drill-string modeling, Sarker et al. [15] presented a lumped-section model for axial and torsional motion. The bond graph model was composed of element sections represented by rigid bodies, and the elasticity is modeled by coupled torsional, axial, and bending springs between each body. The model was used to analyze drill-string behavior in horizontal wells and predicted the whirl phenomena close to the BHA together with the coupling effect for axial and torsional motion through the bit-rock model.

In this article, we derive a lumped model for approximating the distributed properties of a drill string in terms of the floating frame of reference formulation, taking advantage of the assumed mode method to describe the elastic deformation. The modeling framework in terms of Lagrangian dynamics is useful for establishing a component-based model of a generic drill-string, undergoing deformation. The mode shape functions are derived by solving the eigenvalue problem for a fixed-free Euler Bernoulli beam, longitudinal bar, and a shaft. Furthermore, we propose to structure the model in a component-based framework, which is useful for system simulation studies. The proposed model and the structure are analyzed in simulations, and the coupling effect is extensively investigated through specific case studies, relevant for a drilling system.^²

2 Modeling of a Drill-String Differential Element

We assume that no rotation of the drive machinery occurs due to the fixed conditions at the dolly guides. The dolly guides are the vertical rails to which the top drive is fixed, in the derrick tower. Furthermore, we assume a vertical well profile, that the assembly consists of drill pipes, and that the drill-string can make contact with the wellbore wall.

With Fig. 1 as a starting point, we consider an inertial frame located in space, a body frame attached to a pipe differential element, and a deformed configuration in frame {x, y, z}₂. These frames are denoted {x, y, z}₀, {x, y, z}₁, and {x, y, z}₂, respectively. The frame and vector definitions are shown in Fig. 2. We assume that the elastic, deformed element is skewed in the horizontal plane, with no rotation around the x and y axes. This assumption simplifies the analysis since both pitch and roll (inclination and azimuth) of the element are neglected. This is reasonable in the current case, considering a vertical well.

Fig. 2

View large Download slide

Orientation of kinematic frames

This formulation is known as the floating reference frame formulation [17]. This gives two sets of coordinates, one describing the rigid body motion, along with the coordinates describing the elastic deformation. In addition, we investigate the effect of coupling the lateral deformation to twisting of the pipe due to torsion. Through the kinematics of the element and the definition of kinetic and potential energies of the system, the equations of motion of the drill string will be expressed as

M (q) \ddot{q} + C (q, \dot{q}) \dot{q} + K q = τ

(1)

where

M (q) \in {I R}^{n_{q} \times n_{q}}

is the system mass matrix,

C (q, \dot{q}) \in {I R}^{n_{q} \times n_{q}}

is the Coriolis-centripetal matrix,

K \in {I R}^{n_{q} \times n_{q}}

is the system stiffness matrix,

τ \in {I R}^{n_{q} \times 1}

is the generalized forces and torques acting on the system,

q \in {I R}^{n_{q} \times 1}

is the system generalized coordinates characterizing the rigid and elastic motion, and M = {m_k,j(q)},

C (q, \dot{q}) = {c_{k, j} (q, \dot{q})}

with k, j = 1, …, n_q, where n_q is the number of generalized coordinates.

We start by defining the position of the pipe segment in space and the transformations occurring in the system, to finally derive the velocity of a point on the pipe element, subjected to elastic deformation. The position of body frame {1} relative to inertial frame {0} is given by

r_{1}^{0} = [x (t), y (t), z (t) + \bar{z}]^{⊤}

⁠, where

\bar{z}

is the arc length along the rigid drill-string in the underformed case, and the bold notation is used for matrices and vectors. The rotation matrices from frame {0} to {1} and {1} to {2} are given as follows:

R_{1}^{0} (θ (t)) = [\begin{matrix} c θ & - s θ & 0 \\ s θ & c θ & 0 \\ 0 & 0 & 1 \end{matrix}], R_{2}^{1} (ψ (\bar{z}, t)) = [\begin{matrix} c ψ & - s ψ & 0 \\ s ψ & c ψ & 0 \\ 0 & 0 & 1 \end{matrix}]

(2)

where cθ = cosθ, sθ = sinθ, and θ(t) is the rotation angle from {0} to {1} around the z₁-axis and

ψ (\bar{z}, t)

is the rotation angle from {1} to {2} around the z₂-axis, given in terms of both space and time.

The position of an arbitrary twisted point

p_{2, t}^{0}

in frame {2} expressed in {0} is given by

p_{2, t}^{0} = r_{1}^{0} + R_{1}^{0} (p_{2, e}^{1} + R_{2}^{1} p_{2, t}^{2})

(3)

where

p_{2, e}^{1} = [u_{e x}, u_{e y}, u_{e z}]^{⊤}

is the position of the origin of frame {2} relative to the {1} frame origin,

p_{2, t}^{2}

is the position of a point undergoing torsion relative to the {2} frame origin, and subscripts e, t refer to elastic and torsion, respectively, see Fig. 2. The twisted pipe material point

p_{2, t}^{2}

has an additional angle

ψ (\bar{z}, t)

due to torsional deformation, as illustrated in x, y planes in Fig. 3.

Fig. 3

View large Download slide

Plane view of the pipe differential element

From Fig. 3, we can define the vector describing the twisted point

p_{2, t}^{2} = [\tilde{r}, 0, 0]^{⊤}

⁠, where

\tilde{r}

is the averaged radius of the pipe. The time derivative of a material point on the pipe element is derived as follows:

{\dot{p}}_{2, t}^{0} = {\dot{r}}_{1}^{0} + (B_{θ} + S_{θ}) \dot{θ} + S_{ψ} \dot{ψ} + R_{1}^{0} {\dot{p}}_{2, e}^{1}

(4)

where

B_{θ} = (\partial R_{1}^{0} / \partial θ) p_{2, e}^{1}

⁠,

S_{θ} = (\partial R_{1}^{0} / \partial θ) R_{2}^{1} p_{2, t}^{2}

⁠,

S_{ψ} = R_{1}^{0} (\partial R_{2}^{1} /

\partial ψ) p_{2, t}^{2}

⁠,

\dot{ψ} = d ψ (\bar{z}, t) / d t = \partial ψ (\bar{z}, t) / \partial t

since

d \bar{z} / d t = 0

⁠,

{\dot{p}}_{2, t}^{2} = 0

since

\tilde{r}

is a constant, and

\dot{R} = (\partial R / \partial θ) \dot{θ}

⁠.

The flexibility of the drill-string is represented by the elastic deformation in

p_{2, e}^{1} (\bar{z}, t)

and

ψ (\bar{z}, t)

⁠. In the assumed mode method, these spatial variables can be described utilizing mode shape functions obtained from solving the eigenvalue problem of lateral bending by means of the Euler–Bernoulli beam model, longitudinal vibration of a bar, and torsional vibration of a shaft. By using the separation of variables technique for an unforced beam, the continuous variables in x, y, and z along with the twist are given by the product of a mode shape dependent only on

\bar{z}

and a function depending on t given as follows:

u_{e, l} (\bar{z}, t) = N_{l} (\bar{z}) h_{l} (t), ψ (\bar{z}, t) = Y (\bar{z}) h_{r} (t)

(5)

where N_l and Y are the mode shape function required for satisfying the boundary conditions, h_l and h_r are time-dependent functions required to satisfy the initial conditions, l ∈ x, y, z, and subscript r denotes rotation. With these definitions, the resulting functional description of the mode shape functions are then derived from a set of boundary conditions [18].

In this work, we apply fixed-free boundary conditions, i.e., u_e,l(0, t) = ψ(0, t) = 0 and assume zero bending moment, shear forces, moment, and stress at

\bar{z} = L

⁠, where L is the length of the drill-string. For lateral bending (l = x, y), the resulting kth mode shape can be given by the normalized length

s = \bar{z} / L

as follows:

N_{k, l} (s) = \cos (β_{k} s) - \cosh (β_{k} s) - b_{k} (\sin (β_{k} s) - \sinh (β_{k} s))

(6)

where l = x, y, b_k = (cos(β_k) + cosh(β_k))/(sin(β_k) + sinh(β_k)), and β_k is obtained from solving the corresponding frequency equation

\cos (β) \cosh (β) + 1 = 0

(7)

yielding β_k, k = 1, …, ∞, and the natural frequency associated with each mode can be found from

ω_{k, x / y} = β_{k}^{2} \sqrt{E I_{x / y} / ρ A L^{4}}

⁠, where E is the Young’s modulus, I_x/y is the second moment of area about x or y axis, and A is the cross-sectional area.

For the longitudinal and torsional coordinates, the mode shape functions are given by the same functional description. The kth mode shape are found by solving the frequency equation cos (γ) = 0, yielding γ_k = (2k − 1)π/2, k = 1, …, ∞, and the natural frequencies are given as

ω_{k, p / s}^{2} = γ_{k}^{2} c_{p / s}^{2}

and the subscript p/s denotes either a speed of sound of a pressure wave in the material (longitudinally) or a shear wave (torsion). The mode shape and associated natural frequency are given as follows:

\begin{aligned} N_{k, z} (s) & = Y_{k} (s) = \sin (\frac{(2 k - 1) π}{2} s) \\ ω_{k, p / s} & = c_{p / s} \frac{(2 k - 1)}{2 L} π, k = 1, \dots, \infty \end{aligned}

(8)

where

c_{p} = \sqrt{E / ρ}

⁠,

c_{s} = \sqrt{G / ρ}

⁠, G is the shear modulus, and ρ is the density of the material. Furthermore, the mode shapes are orthogonal and satisfy:

\int_{0}^{1} N_{k, l} N_{w, l} d s = \int_{0}^{1} Y_{k} Y_{w} d s = 0, k \neq w

(9)

The first three numerical values of β_k, the mode shapes N_k,x, N_k,y for lateral bending, and N_k,z = Y_k for longitudinal and torsional vibrations with the corresponding natural frequencies are presented in Fig. 4.

Fig. 4

View large Download slide

First three natural frequencies and mode shapes fixed-free drill string in lateral, longitudinal, and torsional deformation. The nodes are marked with solid circles and antinodes with a cross.

With the assumed mode method, the elastic motion of the drill string during a time interval can be given as a linear combination of the mode shapes in Eqs. and [19,20]. As such, we can describe the spatial variables

p_{2, e}^{1} (\bar{z}, t)

and

ψ (\bar{z}, t)

with

u_{e, l} (\bar{z}, t) = \sum_{k = 1}^{n_{l}} N_{k, l} (\bar{z}) q_{k l} (t)

(10)

ψ (\bar{z}, t) = \sum_{k = 1}^{n_{r}} Y_{k} (\bar{z}) q_{k} (t)

(11)

where q_kl(t) are independent generalized coordinates in bending and longitudinal motions, and q_k(t) in torsional motion. These coordinates are included in the total set of system states, obtained by solving for q in Eq. .

We can rearrange Eqs. and in the vector form as follows:

p_{2, e}^{1} = [\begin{matrix} N_{x} & 0 & 0 \\ 0 & N_{y} & 0 \\ 0 & 0 & N_{z} \end{matrix}] [\begin{matrix} q_{e x} \\ q_{e y} \\ q_{e z} \end{matrix}] = N q_{e}, ψ (\bar{z}, t) = Y q_{ψ}

(12)

where

N_{l} \in {I R}^{1 \times n_{l}}

⁠,

Y \in {I R}^{1 \times n_{r}}

⁠,

q_{e} = [q_{e x}^{⊤}, q_{e y}^{⊤}, q_{e z}^{⊤}]^{⊤}

with

q_{l} = [q_{1 l}, q_{2 l}, \dots q_{n_{l} l}]^{⊤}

⁠,

q_{ψ} = [q_{1}, q_{2}, \dots q_{n_{r}}]^{⊤}

⁠, and the mode shape matrix N and vector Y are defined as follows:

\begin{aligned} N & = [\begin{matrix} [N_{1, x}, \dots, N_{n_{x}, x}] & 0 & 0 \\ 0 & [N_{1, y}, \dots, N_{n_{y}, y}] & 0 \\ 0 & 0 & [N_{1, z}, \dots, N_{n_{z}, z}] \end{matrix}], \\ Y & = [Y_{1} \dots Y_{n_{r}}] \end{aligned}

(13)

where the total number of generalized coordinates for the mode shapes is then n_modes = n_e + n_r, where n_e = n_x + n_y + n_z.

By using the definition of the spatial coordinates

p_{2, e}^{1}

and ψ, the velocity of the material point on the pipe element can be derived by substituting Eq. into Eq. and using the defined generalized coordinates, yielding

{\dot{p}}_{2, t}^{0} (q, \dot{q}) = L \dot{q}

(14)

where

L = [I H (q) R_{1}^{0} N S_{ψ} Y]

⁠, I is the identity matrix,

H (q) = B_{θ} + S_{θ}

⁠, and

q = [x, y, z, θ, q_{e}^{⊤}, q_{ψ}^{⊤}]^{⊤}

are the drill-string generalized coordinates. The configuration gives a total set of

n_{q} = \dim (r_{1}^{0}) + n_{modes} + 1

generalized coordinates.

Furthermore, inserting Eq. into Eq. , the position and rotation of the pipe element in space and time is then obtained by summation as follows:

\begin{aligned} p_{2, t}^{0} (\bar{z}, t) & = r_{1}^{0} + R_{1}^{0} N (\bar{z}) q_{e} (t) + R_{1}^{0} R_{2}^{1} (Y (\bar{z}) q_{ψ} (t)) p_{2, t}^{2}, \\ ϕ (\bar{z}, t) & = θ (t) + Y (\bar{z}) q_{ψ} (t) \end{aligned}

(15)

at each sample time, and

ϕ (\bar{z}, t)

is the twist of the drill-string at the distance

\bar{z}

⁠.

3 Derivation of System Energy

The kinetic energy due to linear and angular velocity and potential energy in the form of strain energy for the pipe element dz are derived in this section. The kinetic energy T of the material point on the pipe differential element dz is given by

T = \frac{1}{2} \int_{p}^{} {\dot{p}}_{2, t}^{0^{⊤}} {\dot{p}}_{2, t}^{0} d m

(16)

We further assume that the density and the cross-sectional area of the pipe is uniform over the entire length. Integrating Eq. around the pipe circumference for each mass point

p_{2, t}^{0}

and then integrating along the drill-string axis

\bar{z}

give the total kinetic energy. The infinitesimal area is derived according to Fig. 3 as

d A = \tilde{r} t_{w} d α

⁠, where α is an arc angle and t_w is the pipe wall thickness. The element volume is then

d V = d A d \bar{z}

⁠, and the mass element dm is defined as

d m = ρ d V = ρ t_{w} \tilde{r} d α d \bar{z}

⁠. This allows us to rewrite Eq. in terms of the generalized coordinates as follows:

T = \frac{1}{2} {\dot{q}}^{⊤} M \dot{q}

(17)

where

M = M^{⊤}

is the system mass matrix for the pipe element, i.e.,

M = \tilde{r} ρ t_{w} \int_{\bar{z}}^{} \int_{α}^{} [\begin{matrix} I & H & R_{1}^{0} N & S_{ψ} Y \\ \dots & H^{⊤} H & H^{⊤} R_{1}^{0} N & H^{⊤} S_{ψ} Y \\ \dots & \dots & N^{⊤} N & N^{⊤} {R_{1}^{0}}^{⊤} S_{ψ} Y \\ \dots & \dots & \dots & Y^{⊤} S_{ψ}^{⊤} S_{ψ} Y \end{matrix}] d α d \bar{z}

(18)

where

R^{⊤} R = I

has been used. The dimension of the mass matrix is n_q × n_q depending on the chosen number of generalized coordinates. The moment of inertia of the pipe is then the integrated element

I = H^{⊤} H

⁠.

The potential energy due to strain from deformation of the element is due to the property of the pipe to resist deformation. The strain energy of the element, assuming a symmetric cross section, can be defined for each of the distinctive axes in frame {2} [11,12], yielding

U_{x} = \frac{1}{2} E I_{x} (\frac{\partial^{2} u_{e x}}{\partial {\bar{z}}^{2}})^{2}, U_{y} = \frac{1}{2} E I_{y} (\frac{\partial^{2} u_{e y}}{\partial {\bar{z}}^{2}})^{2}

(19)

U_{z} = \frac{1}{2} E A (\frac{\partial u_{e z}}{\partial \bar{z}})^{2}, U_{ψ} = \frac{1}{2} G J (\frac{\partial ψ}{\partial \bar{z}})^{2}

(20)

where J is the torsional constant, which in this case is equal to the second moment of area about z axis, and the products EI_x/y, EA, GJ denote the flexural, axial, and torsional rigidity, respectively. By using the definitions from Eq. , we can express the derivatives with respect to

\bar{z}

as follows:

\frac{\partial p_{2, e}^{1}}{\partial \bar{z}} = \frac{1}{L} N^{^{'}} q_{e}, \frac{\partial^{2} p_{2, e}^{1}}{\partial {\bar{z}}^{2}} = \frac{1}{L^{2}} N^{^{″}} q_{e}, \frac{\partial ψ}{\partial \bar{z}} = \frac{1}{L} Y^{^{'}} q_{ψ}

(21)

where the prime ′ and ″ represent first and second derivatives with respect to

\bar{z}

⁠. Due to the uniform cross-sectional area, we have I_x = I_y, and J = I_x + I_y. The strain energy due to bending, longitudinal stress, and torsion using Eq. is given by

\begin{aligned} U & = \frac{1}{2} q_{e x}^{⊤} \frac{E I_{x}}{L^{3}} \int_{0}^{1} N_{x}^{^{″} ⊤} N_{x}^{^{″}} d s q_{e x} + q_{e y}^{⊤} \frac{E I_{y}}{L^{3}} \int_{0}^{1} N_{y}^{^{″} ⊤} N_{y}^{^{″}} d s q_{e y} \\ + q_{e z}^{⊤} \frac{E A}{L} \int_{0}^{1} N_{z}^{^{'} ⊤} N_{z}^{^{'}} d s q_{e z} + q_{ψ}^{⊤} \frac{G J}{L} \int_{0}^{1} Y^{^{'} ⊤} Y^{^{'}} d s q_{ψ} \end{aligned}

(22)

where N_j are the row vectors from Eq. . Due to the fact that elastic deformation in x and y direction are represented by the same mode shape functions, i.e.,

N_{x}^{^{″}} = N_{y}^{^{″}}

⁠, we can then rewrite the total strain potential energy as follows:

U = \frac{1}{2} (q_{e x}^{⊤} K_{1} q_{e x} + q_{e y}^{⊤} K_{1} q_{e y} + q_{e z}^{⊤} K_{2} q_{e z} + q_{ψ}^{⊤} K_{ψ} q_{ψ})

(23)

The total strain energy in the matrix form can be written as follows:

U_{e} = \frac{1}{2} q^{⊤} [\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & K_{e} & 0 \\ 0 & 0 & 0 & K_{ψ} \end{matrix}] q = \frac{1}{2} q^{⊤} K q

(24)

where

K = K^{⊤}

⁠, and the stiffness matrix K_e = diag(K₁, K₁, K₂). The gravitational potential energy is derived in the form of an external force in Sec. 4 and not included in the total potential energy.

4 External and Applied Forces and Torques

The external and applied forces and torques characterize the interaction between the drill-string with the system on the rig floor, influence of gravity, wellbore dynamics, and the contact between the bit and the rock. Consider the virtual work of a force acting on dm at

p_{2, t}^{0}

given as follows [20]:

δ W = f^{⊤} δ p_{2, t}^{0}

(25)

where δW is the virtual work of the force

f = [F_{x}, F_{y}, F_{z}]^{⊤}

⁠. Since the virtual displacement is given in terms of the generalized coordinates, i.e.,

δ p_{2, t}^{0} = \sum_{k = 1}^{n_{q}} (\partial p_{2, t}^{0} / \partial q_{k}) δ q_{k} = L δ q

⁠, the virtual work in Eq. is rewritten as follows:

δ W = f^{⊤} L δ q

(26)

where

L = \partial p_{2, t}^{0} / \partial q^{⊤}

can be interpreted as a projection matrix relating the forces and torques to the generalized coordinates. Furthermore, the generalized forces acting on the drill string element can be defined as follows:

τ = f^{⊤} L

(27)

and the virtual work can then be rewritten as δW = τδq.

The external forces acting on the body are the gravitational force and impact force from the borehole wall. Since we consider a pipe deformation, which is subject to lateral movement with no rotation about y₂ or x₂, we assume that z₁ and z₂ are oriented vertically downward. The gravitational force acting on dm as

f_{g} = [0, 0, g d m]^{⊤}

⁠, and the generalized gravitational force is found by integrating over the drill-string length as follows:

τ_{g}^{e} = \tilde{r} ρ t_{w} L \int_{0}^{1} \int_{0}^{2 π} f_{g}^{⊤} L d α d s

(28)

where superscript e denotes external.

Inspecting the matrix product in Eq. , we find that the last row of H in L from Eq. is zero and the product

f_{g}^{⊤} R_{1}^{0} N

evaluates to a row vector [0, 0, N_z]. The generalized gravitational forces in Eq. can then be written as follows:

τ_{g}^{e} = 2 π t_{w} \tilde{r} ρ L {[\begin{matrix} 0, & 0, & g, & 0, & 0, 0, g N_{z}, & 0 \end{matrix}]}^{⊤}

(29)

where

τ_{g}^{e}

is the vector of generalized gravitational forces.

The drill-string is subject to lateral motion and hence is constrained by the borehole wall. For the interactions occurring while drilling, we shall assume that the deformation of the borehole wall is negligible, and we only consider the impact response, i.e., the contact forces that occur when the drill string hits the wall.

The impact with the wall is modeled by a set of springs, which set up a reaction force. In this article, the contact dynamics are adopted from the Hertzian contact law in Ref. [10]. Referring to Fig. 1, the impact occurs when the lateral displacement of the drill-string is larger than the clearance variable l_c. The discontinuous function for lateral contact is defined as follows:

F_{c} (z_{i}, t) = {\begin{matrix} k_{c} (χ (z_{i}, t) - l_{c})^{3 / 2}, & χ (z_{i}, t) \geq l_{c} \\ 0, & χ (z_{i}, t) < l_{c} \end{matrix}

(30)

where

χ (z_{i}, t) = ‖ p_{2 e, x / y}^{0} ‖

is the magnitude of the x and y components of

p_{2 e, x / y}^{0}

⁠, i.e., the center of the drill-string differential element, at z_i for a given time instant. Furthermore, the generalized force governing the Hertzian contact law for a point i along the z axis in the borehole can then be given as follows:

τ_{c, i}^{e} = \tilde{r} \int_{0}^{1} \int_{0}^{2 π} {f_{c}}^{⊤} L δ (s - \frac{z_{i}}{L}) d α d s

(31)

where δ(s − z_i/L) is the Dirac delta function and

f_{c} = - F_{c} (z_{i}, t) [p_{2 e, x}^{0}, p_{2 e, y}^{0}, 0]^{⊤} / ‖ p_{2 e, x / y}^{0} ‖

is the vector of contact forces projected in the x and y direction.

4.1 Nonconservative Forces.

Viscous damping is included to simulate the dissipation of energy in the vibrating body. The viscous damping in vertical wells includes the effect of drilling mud between the borehole and drill-string [11]. The damping matrix is defined as

D = diag (R_{1}, R_{2}, \dots, R_{n_{q}})

⁠. Following Ref. [19], that each structure dissipates energy when excited, we define the structural damping R_i for each mode as follows:

R_{q_{e}} = 2 ξ_{e} \sqrt{k_{i} m_{i}}, R_{ψ} = 2 ξ_{ψ} \sqrt{c_{i} I_{i}}

(32)

where ξ is the relative damping ratio, k_i is the equivalent stiffness for axial and lateral modes, c_i is the equivalent torsional stiffness, and I_i is the modal moment of inertia. The vector of damping forces is then

τ_{nc}^{d} = - D \dot{q}

⁠, with subscript nc for nonconservative.

4.2 Actuating Forces.

The actuating forces are the forces at the top of the pipe due to hoisting or lowering, and the applied torque on the pipe. The vector $f_{0} = [F_{x}, F_{y}, F_{z}]^{⊤}$ is the vector of actuating forces expressed in the inertial frame. The two forces perpendicular to the element on the x and y axes are restoring forces, which hold the pipe in place. The definition of the forces acting on the pipe element is illustrated in Fig. 5.

Fig. 5

View large Download slide

Forces distributed on the element

The forcing function directing f₀ to the point z_i on the pipe is the Dirac function δ(s − z_i/L). The virtual work of the actuating forces for F⁰ is derived in a similar manner as Eq. , yielding

{τ_{f}^{a}}^{⊤} = \tilde{r} \int_{0}^{1} \int_{0}^{2 π} {f_{0}}^{⊤} L δ (s - \frac{z_{i}}{L}) d α d s

(33)

Since the Dirac function is used to only evaluate the integral when s = z_i/L, we can rewrite Eq. as follows:

{τ_{f}^{a}}^{⊤} = {f_{0}}^{⊤} Γ (z_{i}, q)

(34)

where

Γ (z_{i}, q) = 2 π \tilde{r} L (z_{i} / L)

⁠.

The applied torque on the pipe is assumed to independently interact with θ and ψ, since the rotation axes are parallel. We then define the virtual work of the torque and the generalized torque as follows:

δ W_{t}^{a} = T_{0} W δ q, {τ_{t}^{a}}^{⊤} = T_{0} Λ (z_{i})

(35)

where T₀ is a scalar, W = [0, 1, 0, Y],

{τ_{t}^{a}}^{⊤}

is the generalized torque vector, and

Λ (z_{i}) = \int_{0}^{1} W δ (s - z_{i} / L) d s

⁠.

5 Equations of Motion

The Lagrangian of the system is defined as L = T − U. Inserting for the kinetic and strain potential energy into the Lagrangian yields

L = \frac{1}{2} {\dot{q}}^{⊤} M (q) \dot{q} - \frac{1}{2} q^{⊤} K q

(36)

The equations of motion for the system are then derived according to the expression

\frac{d}{d t} \frac{\partial L}{\partial {\dot{q}}_{k}} - \frac{\partial L}{\partial q_{k}} = τ_{k}, k = 1, \dots, n_{q}

(37)

where τ_k is a generalized force acting on q_k. This gives a set of n_q ODEs for the system. Performing the partial differentiation of the Lagrangian in Eq. and taking the time derivative of

\partial L / \partial {\dot{q}}_{k}

yield the equation of motion

M (q) \ddot{q} + C (q, \dot{q}) \dot{q} + K q = τ

(38)

where

M \ddot{q}

is due to inertial acceleration of the element,

C (q, \dot{q}) \dot{q}

is the Coriolis-centripetal forces, and τ is the generalized forces and torques. The matrices M and K are given by the expressions in Eqs. and . The representation of the C matrix is not uniquely defined. We derive the Coriolis-centripetal matrix from the elements of Eq. , according to Ref. [20], given as follows:

c_{k, j} = \sum_{i = 1}^{n_{q}} \frac{1}{2} (\frac{\partial m_{k, j}}{\partial q_{i}} + \frac{\partial m_{k, i}}{\partial q_{j}} - \frac{\partial m_{i, j}}{\partial q_{k}}) {\dot{q}}_{i}

(39)

where k, j = 1, …, n_q, and the term between the brackets in Eq. is referred to as the Christoffel symbols of the first kind. Furthermore, the generalized forces and torques in τ are defined as follows:

τ = τ_{g}^{e} + τ_{f}^{a} (f_{0}) + τ_{t}^{a} (T_{0}) + τ_{c}^{e} (f_{c}) + τ_{nc}^{d} (\dot{q})

(40)

6 Bond Graph Model

To formulate the equations of motion in a component-based framework, we utilize the BG methodology [19]. A BG model is derived from Eqs. and according to the generalized momentum p and displacement q in Hamiltonian form [21]. This gives a field representation of the system. The generalized momentum vector is given as

p = \partial L / \partial \dot{q}

⁠, and its time derivative is given as follows:

\dot{p} = \frac{d}{d t} (\frac{\partial L}{\partial \dot{q}}) = \frac{\partial T}{\partial q} - \frac{\partial U}{\partial q} + τ

(41)

from which we can define

e^{c} (q, \dot{q}) = \partial T / \partial q - \partial U / \partial q

as the corresponding effort vector for inertia forces and torques such as Coriolis forces in the system, and restoring forces and torques. Given the generalized momentum vector, and noting that

\partial L / \partial \dot{q} = M \dot{q}

⁠, we can express the generalized velocities as follows:

\dot{q} = M (q)^{- 1} p

(42)

which together with Eq. form the state-space model in terms of BG vector elements.

The gravitational force is represented as an effort source Se and the damping from $τ_{nc}^{d}$ is modeled as an R-field. The external forces acting on the model is distributed with a modulated transformer, denoted MTF, with modulus Γ(q) and Λ. This constitutive relation holds if the transformer is power conservative, yielding $F \dot{x} = τ_{out} {\dot{q}}_{out}$ [19].

The final causal BG model for the drill-string component is shown in Fig. 6. The graph has complete integral causality, constraining the subsystems to supply generalized forces as input. This is the advantage of utilizing the BG methodology, where the complex model is formed and now any relevant subsystem can be attached to the model. Note that the bond $[e_{1}, f_{1}]^{⊤}$ is split into $[e_{6}, f_{6}]^{⊤}$ and $[e_{2}, f_{2}]^{⊤}$ to divide the graph into the rigid and elastic states, and that the restoring forces in $e^{c} (q, \dot{q})$ are separated out by including a compliance element.

Fig. 6

View large Download slide

Bond graph model of the dynamic drill-string model

7 Simulation Results

In this section, we investigate the model performance through four simulation case studies, denoted Case 1–4. Case 1 includes an analysis of the input effect on the individual mode shapes. In Case 2, we perform tracking of defined revolutions-per-minute (RPM) set-points. In Case 3, the effect of reducing the number of modes is investigated, and in Case 4, a friction model is included to analyze stiction on the end of the drill-string. Initially, we define our system with four modes in u_e,l and with two torsional modes in ψ. This gives a set of 18 generalized coordinates for the system. The drill-string parameters used in this section are presented in Table 1.

Table 1

Drill-string structural parameters

L	ρ	OD	ID	E	G	g
(m)	(kg/m²)	(m)	(m)	(GPa)	(GPa)	(m/s²)
1000	7850	0.127	0.107	200	79.3	9.81

The implementation and numerical integration for the BG model in Fig. 6 are done in 20-sim, using the Vode Adams algorithm. The initial conditions for the case studies are $q_{e} (0) = K_{e}^{- 1} τ_{g}^{e}$ ⁠, r(0) = 0, θ(0) = 0, and $q_{ψ} (0) = 0$ ⁠.

The force input in the lateral and vertical plane is set in the inertial reference frame and transformed to the generalized coordinates by $τ = f^{⊤} Γ (q)$ ⁠. Viscous damping is included in the system, given as $T_{f} = \dot{ψ} (\bar{z}, t) c_{v}$ ⁠, where $\dot{ψ} (\bar{z}, t) = \partial ψ (\bar{z}, t) / \partial t$ and c_v is chosen to reflect the damping due to drag forces in the wellbore.

The system is modeled with proportional damping from Eq. (32). Furthermore, the relative damping ratio ξ_i for the modes is tuned in simulation to include damping in the oscillations due to vibration.

The contact forces from $τ_{c, i}^{e}$ define the boundary to the wellbore, starting at z = 300 m being distributed by eight contact points. Furthermore, l_c = 0.1 m, k_c = 2500 · 10³ N/m being drawn from Refs. [7,22].

In Case 1, direct actuation of the top of string torque is done. This is included in the BG model by directly manipulating T₀, which is transformed by Λ in Eq. (35). In addition, we implement a more realistic drive based on the work in Ref. [23], in terms of BG elements depicted in Fig. 7. The controller constitutes then a proportional-integral control law. The rightmost bond in Fig. 7 is then connected to Λ MTF (by T₀) in Fig. 6.

Fig. 7

View large Download slide

Bond graph model of the top drive

Finally, the input signal from the reference source (set-point control) is smoothed through a second-order filter, in the form

H_{f} (s) = \frac{k_{f} ω_{n}^{2}}{s^{2} + 2 ζ ω_{n} + ω_{n}^{2}}

(43)

where ζ is the damping ratio and ω_n is the frequency of the signal.

7.1 Case 1: Modal Contribution.

To evaluate the modal contribution, we set the desired rotational speed in the drive to ω_sp = 5 rad/s (approximately 48 RPM), simulating the model for 20 s, and perform a similar experiment with constant torque T₀ = 100 Nm, i.e., no dynamic drive model. A stiff controller is chosen for the system, with $K_{p} = 10 J \sqrt{G ρ_{s}}$ and $K_{i} = 2 I_{td} = 2 (I_{r} + I_{g} n_{g}^{2})$ ⁠, where I_r = 200 kgm² and I_g = 20 kg/m² is the rotor and gearbox inertias, respectively, n_g = 4 is the gear ratio, and l_s = 1 m is the connection shaft length. We include viscous friction at z = L, i.e., downhole with c_v,L = 0.05 Nms/rad.

The number of modes required depends on how accurately we want to capture the frequencies being excited by the system. Since the longitudinal deformation is uncoupled from the lateral and torsional deformation, we analyze the results of actuating the string in the vertical direction. The root-mean-square (RMS) values of the generalized modal coordinates describing the individual effect on each mode are listed in Table 2.

Table 2

Modal RMS values for case 1, without a drive and with a drive (D, 1–4)

	RMS (q_ex)	RMS (q_ey)	RMS (q_ez)	RMS (q_ψ)
Mode	(m)	(m)	(m)	(rad)
1	0.52	0.53	0.96	0.46
2	0.29	0.27	0.053	0.54
3	0.15	0.14	0.017	N/A
4	0.065	0.066	0.0083	N/A
D,1	0.48	0.5	0.96	1.1
D,2	0.26	0.26	0.057	0.12
D,3	0.15	0.15	0.019	N/A
D,4	0.11	0.1	0.0092	N/A

	RMS (q_ex)	RMS (q_ey)	RMS (q_ez)	RMS (q_ψ)
Mode	(m)	(m)	(m)	(rad)
1	0.52	0.53	0.96	0.46
2	0.29	0.27	0.053	0.54
3	0.15	0.14	0.017	N/A
4	0.065	0.066	0.0083	N/A
D,1	0.48	0.5	0.96	1.1
D,2	0.26	0.26	0.057	0.12
D,3	0.15	0.15	0.019	N/A
D,4	0.11	0.1	0.0092	N/A

Note: ξ_x/y = 0, ξ_z = 0, and ξ_ψ = [0.01, 0.05].

From the values in Table 2, we notice that the first and second modes of q_ex and q_ey are excited the most by actuator signals. Comparing between using the drive or directly manipulating the torque T₀ indicates that the second torsional mode is the largest contributor when directly manipulating T₀, and the first torsional mode is larger when including the top-drive model. The modes are based on unforced response at the end boundary. However, a drill-string is normally in tension, and the drill-collar section is compressed when placed on the rock-bottom. In cases 1–3, we only consider the results in terms of off-bottom rotation, with no interaction with vertical forces at the end boundary.

The twist angle $ψ (\bar{z}, t)$ and the twist angular velocity $\dot{ψ} (\bar{z}, t)$ are shown in Fig. 8. The angle of twist reaches steady state, and the angular velocity along the drill-string oscillates around zero due to the nonlinear coupling with the lateral modes.

Fig. 8

View large Download slide

Twist angle and angular velocity, case 1

7.2 Case 2: RPM Increase.

In this case, we set four individual set-points for ω₀ at 4, 8, 16, and 20 rad/s (38–191 RPM, approximately). The viscous friction is distributed along the drill-string with $c_{v, z_{i}} = 2.0 N m s$ and c_v,L = 4.0 Nms/rad to give the drive more resistance in rotating the drill-string. The results are shown in Fig. 9.

Fig. 9

View large Download slide

Angular velocity at z = 0 and z = L, and top-of-string torque T_i,0. Subscript i denotes the four individual set-points.

The plots show a common first-order response due to external viscous torque opposing the rotation. The twist is not large and (as also seen from the magnitudes in case 1, Fig. 8) with the two curves following each other quite closely, as we have not included stiction between the drill-pipe and the formation. In the small subplot in plot four of Fig. 9, there are small oscillations for increasing angular velocity. This can be linked to the motion of $p_{2, e}^{0}$ ⁠, i.e., the center of the drill-string, which is seen at the boundary z = L in Fig. 10 for ω_sp = 4 rad/s and 20 rad/s.

Fig. 10

View large Download slide

Drill-string center at z = L, for ω_sp = 4 and 20 rad/s

At the end boundary, the center of the drill string is oscillating with larger amplitude and frequency for higher RPMs. At this stage, we have not considered the effect of added mass and the effect of fluid in the annulus. It is reasonable to think that both these effects would provide additional damping.

7.3 Case 3: Evaluation of Model Configurations.

In this case, three different models were constructed. Model 1 is the same as the one from the previous case studies. Model 2 is developed with 14DOFs: three coordinates each in u_l,e and one in ψ. Model 3 consist of 12DOFs, with two coordinates each in u_l,e and ψ. In this case, we set $K_{p} = 5 J \sqrt{G ρ}$ ⁠, $c_{v, z_{i}} = 0.01$ ⁠, c_v,L = 10, and $ξ_{ψ} = 0.01$ ⁠. The result of simulating the three model configurations are shown in Fig. 11.

Fig. 11

Model response for nq = 18, 14, and 12, where ω0 and ωL denote angular velocity at z = 0 and z = L, respectively, and uex is the deformation in the x direction. Superscripts indicate the model number.

View large Download slide

Model response for n_q = 18, 14, and 12, where ω₀ and ω_L denote angular velocity at z = 0 and z = L, respectively, and u_ex is the deformation in the x direction. Superscripts indicate the model number.

The first and second plot indicate that the trajectories correspond to the findings in Table 2, where the first three modes play a significant role in the magnitude of lateral deformation. In model 3, where the third mode is not present, the elastic displacement is seen to be larger when the drill string is accelerating. The amplitude of the lateral coordinate is seen to grow for steady-state angular velocity. The twist of the drill-string for model 2, with one modal coordinate in $ψ (\bar{z}, t)$ ⁠, is seen to be lower than the two other models, where the effect of the second mode increases the magnitude of the twist. Furthermore, the contribution of the fourth mode in lateral vibration is small, seen in the second plot. This is also seen for mode 4 in Table 2, when including a top-drive model.

In terms of simulation performance, we can evaluate the ratio of simulation time to actual time defined as t_rt = t_a/t_sim, where t_a is the actual time. The ratio resulted in t_rt = 6.42 (model 1), t_rt = 3.59 (model 2), and t_rt = 1.61 (model 3). Increasing the number of included modes will increase the size of the model matrices, explaining the difference in simulation speed. Model 3 is then assumed to be better suited for real-time applications, while model 1 might be a better approximation of the real vibration characteristics of the drill-string.

7.4 Case 4: Including Stiction.

To demonstrate the effects of drilling with torque on bit, we implement a friction model to predict the torque as a consequence of applied weight on bit. The applied weight on bit F_n (i.e., the normal force acting on the end of the drill-string) is adjusted manually. Furthermore, $c_{v, z_{i}} = 0.01$ and $ξ_{ψ} = 0.001$ ⁠.

The selection of an appropriate friction model to give the physical interaction with the wellbore is not intuitive. Extensive research has been done on this field (see, e.g., Refs. [6,9,15], and the qualitative response depends on several parameters, such as bit configuration and rock strength. In this study, we limit the case to investigate the model response for model 2 configuration (shown in case 3), in the case of the onset of a friction torque on the drill-string end. For this purpose, we consider the LuGre friction model [24]. This model was also used in terms of stick-slip control in Ref. [25] and showed to capture the stick-slip phenomenon qualitatively well.

The LuGre friction model is given as follows [25]:

\dot{φ} - \dot{ϕ} (L, t) = - σ_{0} | \dot{ϕ} (L, t) | φ / g (\dot{ϕ} (L, t))

(44)

g (\dot{ϕ} (L, t)) = μ_{c} + (μ_{s} - μ_{c}) e^{- (\dot{ϕ} (L, t) / v_{s})^{2}}

(45)

where

\dot{φ}

is the rate of change of internal friction displacement, σ₀ is the elasticity of the bristles between the surfaces, and g is a nonlinear term relating to the Stribeck effect, which sets the deflection, which φ approaches in a steady state [24], μ_c is the Coulomb friction coefficient, μ_s is the static friction coefficient, and v_s is the characteristic velocity. The total torque on bit is then expressed as follows:

T_{f} (L, t) = r_{b} F_{n} (σ_{0} φ + σ_{1} \dot{φ})

(46)

where r_b is the bit radius and F_n is the normal force from the weight on bit.

The LuGre friction model can be implemented as a sub-bond graph, actuating a modulated effort source, which generates the torque on bit in Eq. (46). The sub-bond graph is designed as shown in Fig. 12. The velocity in Eq. (44) is given by adding together the flows of the rightmost 0-junction. The sum $σ_{0} φ + σ_{1} \dot{φ}$ are represented by a C and R-element. The nonlinear term in Eq. (44) including $g (\dot{ϕ} (L, t))$ is modeled as a modulated flow source (MSf-element) to set the additional flow to the 0-junction.

Fig. 12

View large Download slide

LuGre friction model in BG elements

Case 4 parameters are summarized in Table 3, and the LuGre model parameters are drawn from Ref. [25]. The drive parameters are the same as those of case 1. The response of applying a weight on bit of F_n = 35 kN is shown in Fig. 13 for desired drive speeds ω_sp = 5 and 10 rad/s.

Fig. 13

View large Download slide

Drill-string sticking to the formation by friction from torque on bit. The notation sp indicates the angular velocity set-point.

Table 3

Case 4 parameters

L	σ₀	σ₁	μ_c	μ_s	v_s	r_b
(m)	(Nm/rad)	(Nms/rad)	(–)	(–)	(rad/s)	(m)
2000	0.3	$2 \sqrt{σ_{0}}$	0.3	0.35	0.01	0.1

The drive rotates the drill-string, increasing the top-of-string torque until the torque from the drill-string overcomes the breakaway torque defined by the friction model. However, for both desired set-points, the drill-string develops sustained oscillations, seen in plot one and two in the first column of Fig. 13. Increasing the desired angular velocity is seen to result in larger amplitude and lower frequency torsional vibrations. The consequences of applied torque on bit are seen clearly in the first seconds of Fig. 13. The lowest end of the drill-string is not rotating until a given time has passed, and the vibrations develop due to the pipe torque not exceeding the friction torque, seen in the third and fourth plot.

The twist and displacement profile along the length of the drill-string are seen for t ≈ 1 s in second column, plot one and two of Fig. 13. It is noticed that the twist angle is approximately zero for $ω_{{sp}_{1}} = 5 rad / s$ at the end of the drill-string, indicating that this point is still stuck in the formation. The characteristics are seen to change for increasing angular velocity, and the displacement profile along the length is more uneven.

8 Discussion

The qualitative results with the friction model are similar to those of other torsional models with stick-slip friction models. However, the parameters used for the LuGre model had to be adjusted beyond the initial values drawn from Ref. [25] to illustrate the stick-slip cycle. As an initial study of the model, the response in terms of formation sticking were captured.

Comparing the assumed mode method to conventional lumped ODE models, the number of modes rather than the number of point masses is chosen, describing the inertia properties of the system. A better estimate of desired system frequencies is then captured using the assumed mode method, by taking advantage of the knowledge of interacting subsystem frequencies [19].

In comparison to FEMs, the model developed on the basis of assumed modes results in a smaller number of generalized coordinates, and this can have an effect on the simulation time. However, an FEM might be better suited if the geometrical properties of the object are complex, since the model properties are assembled from several local elements, compared to the assumed mode method, where the mode shape functions are defined over the entire length of the object [26].

8.1 Simulation Performance.

The implication of the coupling between the lateral and torsional dynamics is a mass matrix with off-diagonal terms. To facilitate the simulation studies in this work the matrices are assembled offline, and the values for q, $\dot{q}$ are inserted into the matrices before integrating over s during simulation run time. A disadvantage of this approach is increased computational time with increasing number of system states.

In Ref. [16], the sine and cosine terms in Eq. (2) for $R_{2}^{1} (ψ (\bar{z}, t)$ were approximated by a Taylor series expansion. The $\sin (ψ (\bar{z}, t))$ and $\cos (ψ (\bar{z}, t))$ were approximated to be able to construct the matrices offline, with integration over s for all elements. Hence, the accuracy was limited by the order of the expansion. With a fourth-order expansion of the cosine and a third-order expansion of sine, the twist coordinate is approximately limited to $ψ (\bar{z}, t) \in (- π / 2, π / 2)$ ⁠.

In terms of the simulation performance parameter t_rt, the recorded factors in Ref. [16] were 1.47 and 7.82 for models 2 and 3, respectively. Comparing these values to the complete model in this study (see Sec. 7.3), it is noted that the number of $q_{ψ}$ coordinates chosen in the Taylor expansion would dictate the simulation speed of the approximated model.

9 Conclusions

In this article, we have derived a drill-string model with coupled torsion and lateral motion in terms of the floating frame of reference formulation using the assumed mode method and additionally presented the model in a component-based structure in terms of the bond graph methodology.

The model performance was tested in four case studies, analyzing the total modal contribution, showing that modes 1–3 of the generalized coordinates were influenced the most by the top-drive dynamics. Limiting the number of modes through geometry conditions and required accuracy is shown to improve simulation speed toward real-time applications. We included the LuGre friction model and presented a bond graph model implementation of it to analyze the qualitative behavior of drill-string formation sticking, with applied torque at the end boundary of the drill-string. The simulation results show that the friction causes twisting along the drill-string length and that the lateral displacement became larger with increased top-of-string angular velocity.

The proposed model of the distributed characteristics of a drill-string is developed in a compact manner and is flexible toward connecting relevant subsystems in system simulation case studies. A top-drive model and a friction model, for evaluation of critical characteristics such as rotational stiction, were readily integrated into the model structure. The work in this article can then contribute to the design and analysis of control laws in mechatronic systems.

Footnote

2

The work in this article is an extension of the work presented in Ref. [16].

Acknowledgment

The research presented in this article has received funding from the Norwegian Research Council, SFI Offshore Mechatronics, project number 237896.

Conflict of Interest

There are no conflicts of interest.

References

1.

Kapitaniak

,

M.

,

Vaziri Hamaneh

,

V.

,

Chávez

,

J. P.

,

Nandakumar

,

K.

, and

Wiercigroch

,

M.

,

2015

, “

Unveiling Complexity of Drill-String Vibrations: Experiments and Modelling

,”

Int. J. Mech. Sci.

,

101–102

, pp.

324

–

337

.

Google Scholar

Crossref

2.

Ritto

,

T. G.

,

Soize

,

C.

, and

Sampaio

,

R.

,

2009

, “

Non-Linear Dynamics of a Drill-String With Uncertain Model of the Bit-Rock Interaction

,”

Int. J. Non-Linear Mech.

,

44

(

8

), pp.

865

–

876

.

Google Scholar

Crossref

3.

Kyllingstad

,

Å.

, and

Halsey

,

G. W.

,

1988

, “

A Study of Slip/Stick Motion of the Bit

,”

SPE Drilling Eng.

,

3

(

4

), pp.

369

–

373

.

Google Scholar

Crossref

4.

Tucker

,

R.

, and

Wang

,

C.

,

2003

, “

Torsional Vibration Control and Cosserat Dynamics of a Drill-Rig Assembly

,”

Meccanica

,

38

(

1

), pp.

143

–

159

.

Google Scholar

Crossref

5.

Jansen

,

D. J.

,

1993

, “

Nonlinear Dynamics of Oilwell Drillstrings

,”

Ph.D. thesis

,

Delft University

,

Delft, Netherlands

.

6.

Richard

,

T.

,

Germay

,

C.

, and

Detournay

,

E.

,

2007

, “

A Simplified Model to Explore the Root Cause of Stick-Slip Vibrations in Drilling Systems With Drag Bits

,”

J. Sound. Vib.

,

305

(

3

), pp.

432

–

456

.

Google Scholar

Crossref

7.

Christoforou

,

A. P.

, and

Yigit

,

A. S.

,

2003

, “

Fully Coupled Vibrations of Actively Controlled Drillstrings

,”

J. Sound. Vib.

,

267

(

5

), pp.

1029

–

1045

.

Google Scholar

Crossref

8.

Bakhtiari-Nejad

,

F.

, and

Hosseinzadeh

,

A.

,

2017

, “

Nonlinear Dynamic Stability Analysis of the Coupled Axial-Torsional Motion of the Rotary Drilling Considering the Effect of Axial Rigid-Body Dynamics

,”

Int. J. Non-Linear Mech.

,

88

, pp.

85

–

96

.

Google Scholar

Crossref

9.

Detournay

,

E.

,

Richard

,

T.

, and

Shepherd

,

M.

,

2008

, “

Drilling Response of Drag Bits: Theory and Experiment

,”

Int. J. Rock Mech. Mining Sci.

,

45

(

8

), pp.

1347

–

1360

.

Google Scholar

Crossref

10.

Spanos

,

P. D.

,

Chevallier

,

A. M.

, and

Politis

,

N. P.

,

2002

, “

Nonlinear Stochastic Drill-String Vibrations

,”

J. Vib. Acoust.

,

124

(

4

), p.

512

.

Google Scholar

Crossref

11.

Feng

,

T.

,

Vadali

,

M.

,

Ma

,

Z.

,

Chen

,

D.

, and

Dykstra

,

J.

,

2017

, “

A Finite Element Method With Full Bit-Force Modeling to Analyze Drillstring Vibration

,”

ASME J. Dyn. Syst. Meas. Control.

,

139

(

9

), p.

091016

.

Google Scholar

Crossref

12.

Caijin

,

Y.

,

Zaibin

,

C.

,

Wei

,

J.

,

Shiquan

,

J.

, and

Gexue

,

R.

,

2015

, “

A Multibody Dynamic Model of Drillstring for Torque and Drag Analysis

,”

ASME J. Offshore. Mech. Arct. Eng.

,

137

(

3

), p.

031403

.

Google Scholar

Crossref

13.

Borutzky

,

W.

,

2009

, “

Bond Graph Modelling and Simulation of Multidisciplinary Systems—An Introduction

,”

Simul. Model. Pract. Theory

,

17

(

1

), pp.

3

–

21

.

Google Scholar

Crossref

14.

Pedersen

,

E.

, and

Polic

,

D.

,

2014

, “

Bond Graph Modeling of Rotordynamic Systems With a Flexible Shaft Including Shear Correction

,”

Simul. Ser.

,

46

(

8

), pp.

205

–

212

.

15.

Sarker

,

M.

,

Rideout

,

D. G.

, and

Butt

,

S. D.

,

2017

, “

Dynamic Model for 3D Motions of a Horizontal Oilwell BHA With Wellbore Stick-Slip Whirl Interaction

,”

J. Petroleum Sci. Eng.

,

157

, pp.

482

–

506

.

Google Scholar

Crossref

16.

Tengesdal

,

N.

,

Holden

,

C.

, and

Pedersen

,

E.

,

2019

, “

Component-Based Modeling and Simulation of Nonlinear Drill-String Dynamics

,”

Proceedings of the ASME 2019 38th International Conference on Ocean, Offshore and Arctic Engineering. Volume 1: Offshore Technology; Offshore Geotechnics

, vol.

1

,

Glasgow, Scotland, UK

,

June 9–14

.

Google Scholar

Crossref

17.

Shabana

,

A. A.

,

2013

,

Dynamics of Multibody Systems

,

Cambridge University Press

,

Cambridge, UK

.

Google Scholar

Crossref

18.

Rao

,

S. S.

,

2007

,

Vibration of Continuous Systems

,

John Wiley & Sons

,

Hoboken, NJ

.

19.

Karnopp

,

D. C.

,

Margolis

,

D. L.

, and

Rosenberg

,

R. C.

,

2012

,

System Dynamics: Modeling, Simulation, and Control of Mechatronic Systems

,

John Wiley & Sons

,

Hoboken, NJ

.

Google Scholar

Crossref

20.

Egeland

,

O.

, and

Gravdahl

,

J. T.

,

2002

,

Modeling and Simulation for Automatic Control

, Vol.

76

,

Marine Cybernetics Trondheim

,

Norway

.

21.

Karnopp

,

D.

,

1982

, “

Dynamic Forces in Multiport Mechanics: Direct and Lagrangian Formulations

,”

J. Franklin Institute

,

314

(

5

), pp.

265

–

270

.

Google Scholar

Crossref

22.

Spanos

,

P. D.

,

Sengupta

,

A. K.

,

Cunningham

,

R. A.

, and

Paslay

,

P. R.

,

1995

, “

Modeling of Roller Cone Bit Lift-Off Dynamics in Rotary Drilling

,”

ASME J. Energy. Res. Technol.

,

117

(

3

), pp.

197

–

207

.

Google Scholar

Crossref

23.

Kyllingstad

,

A.

, and

Nessjøen

,

P. J.

,

2009

, “

A New Stick-Slip Prevention System

,”

SPE/IADC Drilling Conference and Exhibition

,

Amsterdam, The Netherlands

,

March

, pp.

17

–

19

.

Google Scholar

Crossref

24.

Canudas de Wit

,

C.

,

Olsson

,

H.

,

Astrom

,

K.

, and

Lischinsky

,

P.

,

1995

, “

A New Model for Control of Systems With Friction

,”

IEEE. Trans. Automat. Contr.

,

40

(

3

), pp.

419

–

425

.

Google Scholar

Crossref

25.

Canudas-de Wit

,

C.

,

Rubio

,

F. R.

, and

Corchero

,

M. A.

,

2008

, “

D-OSKIL: A New Mechanism for Controlling Stick-Slip Oscillations in Oil Well Drillstrings

,”

IEEE Trans. Control Syst. Technol.

,

16

(

6

), pp.

1177

–

1191

.

Google Scholar

Crossref

26.

Theodore

,

R. J.

, and

Ghosal

,

A.

,

1995

, “

Finite Element Models for Flexible Multilink Manipulators

,”

Int. J. Rob. Res.

,

14

(

2

), pp.

91

–

111

.

Google Scholar

Crossref

Component-Based Modeling and Simulation of Nonlinear Drill-String Dynamics

Abstract

1 Introduction

2 Modeling of a Drill-String Differential Element

3 Derivation of System Energy

4 External and Applied Forces and Torques

4.1 Nonconservative Forces.

4.2 Actuating Forces.

5 Equations of Motion

6 Bond Graph Model

7 Simulation Results

7.1 Case 1: Modal Contribution.

7.2 Case 2: RPM Increase.

7.3 Case 3: Evaluation of Model Configurations.

7.4 Case 4: Including Stiction.

8 Discussion

8.1 Simulation Performance.

9 Conclusions

Footnote

Acknowledgment

Conflict of Interest

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Component-Based Modeling and Simulation of Nonlinear Drill-String Dynamics

Abstract

1 Introduction

2 Modeling of a Drill-String Differential Element

3 Derivation of System Energy

4 External and Applied Forces and Torques

4.1 Nonconservative Forces.

4.2 Actuating Forces.

5 Equations of Motion

6 Bond Graph Model

7 Simulation Results

7.1 Case 1: Modal Contribution.

7.2 Case 2: RPM Increase.

7.3 Case 3: Evaluation of Model Configurations.

7.4 Case 4: Including Stiction.

8 Discussion

8.1 Simulation Performance.

9 Conclusions

Footnote

Acknowledgment

Conflict of Interest

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only