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

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.


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 1 ). 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. 2 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. 3 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, 6 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., 7 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. 8 ). 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. 9 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. 10 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. 11 For deviated wellbore formulations, a two-dimensional (2D) lumped model (LM) was proposed in Zhao et al., 12 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 13 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 twoelement 3D LM for a straight horizontal well was given in Xie et al. 14 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., 15 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: 1. A real-time oriented derivation of a drill string model for directional wells, where arbitrary wellbores can be simulated. 2. 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 16 for a comparison of this method to well-established formulations). 3. 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. 4. Introduction of Generalized-a 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-a.
This work also includes a field-validated model from Gjerstad et al. 17 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. 18 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.

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 fẽ i g, in coordinate free form, yieldingp and its coordinate form is given as p 0 0i = c f h ½ > , where fc(s), f (s), h(s)g 2 C 1 are parameterized by s. In general,p c ki and p i=c ki 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.
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 is the Euclidean norm. A frame local to the origin ofp 0i is derived from the unit tangent, normal and bi-normal vectors for C at the origin ofp 0i (s). The unit tangent vector is given as where 0 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 wherec 1 Á dc 1 =ds = 0, where we have used the definition of curvature at i given as k = dc 1 =ds k k, 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 k = 0. The unit normal and bi-normal vector given as Finally, the representation of a skew-symmetric matrix is denoted by a hat notation yieldinĝ with x 3 y =xy, andx =Àx > , and Equation (7) will be used in the following sections.

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 wheref i ,T i are the external acting forces and torques,f H i , T H i are the inertia forces and torques, andf c i ,T c i 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 byp k , on body i can be expressed as the sum of the partial contribution from each generalized coordinate q r , r = 1, . . . , n q , in the system. Hence, the virtual displacement is defined as where dq r is the virtual displacement of the generalized coordinate q r . Equation (9) is then a vector tangent to the positionp k (q, t) with q being a vector of generalized coordinates. 18 Moreover, assuming that each time derivative of q r , being _ q r , is independent, the system will have n q degrees of freedom. In addition, to satisfy the constraint p k =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.
where n p is the number of particles in i.
Multiplying the virtual displacements with Equation (8) yields where whereṽ c 0i is the linear velocity of the mass center of i relative to the inertial frame 0,ṽ 0i is the angular velocity of i relative to 0, and ∂ṽ c 0i =∂ _ q r = ∂p c 0i =∂q r is obtained by first taking the time derivative ofp c 0i (q, t), being and furthermore, the partial derivative of Equation (13) with respect to _ q r . Summing over the number of bodies i in Equation (11) yields the d'Alembert principle of equation of motion, expressed as 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 by considering that the variations in dq r are arbitrary. In Kane and Levinson, 19 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, 19 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ṽ whereṽ c ir =∂ṽ c 0i =∂ _ q r is the partial linear velocity, v ir =∂ṽ 0i =∂ _ 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 19 denoted Kane's equations, can be obtained from Equation (15) describing the dynamic equilibrium for body i where F ir is all of the generalized active forces and torques acting on body i, and F H ir is the generalized inertia forces and torques for body i, in the direction of each generalized velocity _ q r . The projection of the inertia forces and torques on _ q r can be expressed by the terms in Equation (15) where m i is the mass of i,ã c 0i is the linear acceleration of the mass center,ã 0i is the angular acceleration of i, andĨ c i is the central inertia dyadic of i.
Performing a summation over n for Equations (17)-(20) we get 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 whereṽ c ir ,ṽ ir 2 R 3 3 1 , and the dimensions of J i is 6 3 n q . Hence, if the equation of motion are for a link, J i is the link Jacobian as derived in Egeland and Sagli. 20 The active and inertia forces and torques are then expressed as a column vector yielding where a i=c i , a i 0i 2 R 3 3 1 are the coordinate form linear and angular accelerations of the mass center of i, respectively, v i 0i 2 R 3 3 1 is the angular velocity of i andv i 0i 2 R 3 3 1 is the skew symmetric form as given in Equation (7), I i=c i R 3 3 3 is the inertia dyadic, f i i , t i i 2 R 3 3 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 where t = P n i 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 _ q r , equal to Equation (21).
For systems exhibiting compliance and damping effects, such as a drill string, Equation (25) can be extended to where K, D 2 R n q 3 n q are the stiffness and damping matrix included by superposition, respectively.

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. The representation is based upon the work in Tengesdal et al. 21 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, . . . , 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 ic 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 binormal 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ã i and is denoted as the fixed base frame, with positiveã 3 axis point downward (ã 3 =Àẽ 3 ), frame i located on C is defined by the unit tangent ic 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ã 3 axis, and horizontal reach is defined along theã 1 axis. In addition, the horizontal reach in the direction ofã 2 can be specified for the model. 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 _ φ 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 drill string model can ultimately be described by the equations of motion similar to what was introduced in section 2 as 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 t 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.

Kinematics
The static wellbore curve C is defined by a vector in the coordinates of 0 relative to an inertial frame and the arclength parameter s as where the bold notation denotes a column vector, i = 1, . . . , 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 2 ½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 = il 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 0i (s i ) in the directions of ic 1 , ic 2 , ic 3 . In addition, each body frame is subject to rotation of about ib 1 by f 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 e 1 ib 2 , e 2 ib 3 , whereb 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 where superscript c denotes the geometric center, ic 1 is the unit tangent vector of frame i computed from Equation (3), ic 3 is the unit normal vector computed according to Equation (5), ic 2 = ic 1 3 ic 3 is the bi-normal, and j i ; z i , and x i are the magnitudes in each of the directions. The curvature k at s i is then computed according to Equation (4). Noting that ib we express the position in the coordinates of i as  where s f i = sin (f i ), and c f i = cos (f i ). The position of b i relative to frame 0, in the coordinates of frame 0, given in coordinate form is given as where R 0 i =fr jk g, R 0 i 2 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 denoted special orthogonal group.
3.1.1. Remark. For completely straight wells, we havec 3 andc 2 undefined since k = 0. Hence, for a strictly vertical sectionp = L vrã3 , ic 3 are aligned with theã 1 axis, and for strictly horizontal wellsp=sL hrã1 orp=sL hrã2 , ic 3 =Àã 3 . Since the curve C is static, we assume that the angular velocityṽ 0i of i relative to 0 is zero, and the linear velocity of the origin of b i relative to the base frame 0, is derived as sinceṽ 0i = 0 due to (∂p 0i =∂s)∂s=∂t = 0, andṽ 0i = 0. The generalized coordinates define the configuration of the system, comprising of where q 2 R n q 3 1 , n q is the number of generalized coordinates. Time differentiation of q r in Equation (35) yields the generalized velocities _ q r . Furthermore, partial differentiation of the linear and angular velocities in Equation (33) with respect to the generalized velocities _ q r yields the partial velocities. 19 These velocities define the projection of the external and inertial forces onto each _ q r since the latter has no directional quantity.
A velocity vector can be expressed as and the partial linear and angular velocity associated with b i is derived as from which we can define the projection Jacobian as 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 20 as where _ J = ∂J i =∂q ð Þ_ q is describing the Coriolis effects if drill string eccentricity is taken into account.

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.
The drill string element inertia forces of the center of mass yieldf where a ! c ib i = _ v c ib i , andm 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 22 where d o 2 i, dp=dc is the outer diameter for drill pipe or collar (notation dp/dc) for element i, d i i, dp=dc is the inner diameter, r s is the density of the pipe material, r m is the mud density, and C a is the added mass coefficient.
The inertia torques are given as Furthermore, the central inertia dyadic of the element is given as where I c yy , I c xx , I c zz 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 The restoring torques about the centroid axis are given asT where G is the shear modulus, and J i is the ith 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 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 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 ib 2 and ib 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 where z i and x i are the local perturbations in the direction ofc 2 andc 3 , respectively. The lumped bending moment can be derived considering the potential energy of a rotational spring, located in each frame i, given as where k b i (E, I x=y i , l i ) is the bending stiffness coefficient. Using the small angle approximation, i.e., sin u'u, we assume that the perturbed lateral motion from the initial wellbore configuration is small. Note that u 1, 0 = z 1 =l 1 , and u 2, 0 = x 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 where 02 n q 3 n q , and K2 n q 3 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 23 as where the gravitational forces acting on i in inertial vertical direction can be expressed as 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 0 i 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 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 which is applied at i = 1, and f 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 The applied forces from the block and drive are then 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 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 and c r i is an estimated coefficient. The viscous damping forces and torques can be rewritten using Equations (58), (59) and (36), yielding where D i 2 6 3 6 .

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 binormal plane, i.e., the local ic 2 À ic 3 plane. 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 ofp c ib i in Equation (30). The angle c i = arctan (p c ib i , ½2 =p c ib i , ½1 ) 6 ¼ f 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, jjp c ib i , ½1, 2 jj À r 0 = 0 where r 0 = (d w À d o i, dp=dc )=2 is the nominal clearance for each element on the curve C. The contact normal force arising from impact can be expressed as Figure 6. Drill pipe or collar in contact with wellbore wall. Note that all the forces and torques are lumped to the ith element. The green dot-dashed curve represents the wellbore curve C.
where d c, i = jj c ib i , ½1, 2 jj, with d c, i being the displaced length of the ith center. The second term on the left-hand side of Equation (61) is the formation damping, assumed proportional to the magnitude of _ d 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 26 whereṽ s ib i is the slip velocity. The friction force arising from wellbore contact is dependent on the unit vector of the slip velocity, given as where m k is the coefficient of kinetic friction (also referred to as the Coulomb friction coefficient), and jjf f , i jj 4 m s jjf n, i jj, where m s . m 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 where the rightmost term denotes the viscous torque proportional to element angular velocity. The contact forces and torques are summarized as

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 which can be rewritten using Equations (39) and (60) as Adding in the forces resulting from linear bending in Equation (51), and using the generalized velocities in Equation (39), we get 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 where M = P n i J > i M i J i 2 R n q 3 n q is the system inertia matrix, C(q, _ q) = P n i J > i M i _ J i 2 R n q 3 n q is the Coriolis matrix, D( _ q) = P n i J > i D ( _ q) i J i 2 R n q 3 n q is the system damping matrix, and t = t a + t g + t e + t c are the generalized forces and torques defined as where t a is the applied forces and torques, t g is the external gravitational forces, t c is the generalized contact forces and torques, and t e represents the generalized tension forces.

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 c i = 0 in this section. When the drill string motion reaches steady state, we have € q(t) = _ q(t) = 0 and Equation (67)

reduces to
Kq eq + t a (q eq ) + t e (q eq ) + t 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.

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 b = 0. The Case 1 parameters are presented in Table 1.
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.
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.

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.
A static force, f e , applied at the bit (s n ) in the direction of ic 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 nc 2 À nc 3 plane. The supporting forces acting on element b n in Equation (67) are given as 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.
The Generalized-a (G a ) scheme (see, e.g., Chung and Hulbert 27 ) is implemented in the LM for nonlinear analysis. The latter method is an extension of the classical method by Newmark 28 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. 29 Furthermore, the parameter p ' 2 ½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).
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.
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 freehanging 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, 30 however, without axial contraction between the model elements.

Simulation and results
In this section, a simulation study is performed, and the qualitative properties of the model are presented. The 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 n1 Support 10 9 , 10 9 , 10 9 N·m À1 k n2 Support 10 9 , 10 9 , 10 9 N·m À1 g Acc. of gravity 9:81·f10 À6 ,10 À4 ,1g m·s À2 f e Static force 0, 4*,   numerical solution forward in time t is obtained using the G a 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 withp 0i (s i ) = L hr s iã3 is used in this section, and the wellbore and drill string parameters are given in Table 3.
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 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 where F s i 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 6 ¼ 0. Hence, for the ith 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 31 with a slip velocity dead zone given as v tr .

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 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 f s = Ð t n g _ f m dt 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 where _ f d is the desired angular velocity at the top-ofstring, 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. 32 is implemented, and the drive parameters are shown in Table 4.

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 4 ).
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. 33 The mud velocity is given by the pipe velocity of the ith element as where 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. 34 It is defined as where DP 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 where g 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 = 8t w, i =(r m v 2 m, i ), yielding where d w À d o, i = D h, i , and t w, i is the wall shear stresses, which in this paper is computed from the work in Gjerstad et al. 17 (section 4.1.1; Appendix I).
An estimate of the bottom hole pressure (BHP) is then obtained as where h n (s n ) = arccos (c 1, 1 =jjc 1 jj) is the inclination at i along C, c 1, 1 is the first component ofc 1 , the first term on the left-hand side is the hydrostatic pressure measured at the true vertical depth (TVD) measured alongã 3 to the position of the element, and the second term is the pressure loss due to inclination.

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. 35 The rig floor position due to heave is given as 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.
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 DP 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. 12 for axial velocities in the vicinity of 2 m/s. From plots two and three (counting top-down) in Figure 12, the clearance variable d 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.

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 m 2 f0:2, 0:3g. Note that BHA eccentricity of e 2 = 0:1 mm is included. The corresponding drive outputs are presented in Figure  14 The tripping velocity is 0:3m 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. 36 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).

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 23 could facilitate further model validation. The rest of this section includes a discussion on model effects and a study of two numerical integration methods.

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 lengthto-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.

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 28 (see, e.g., Wilson and Noynaert, 4 Butlin and Langley 37 ), 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 stickslip behavior may require a relatively short time step (h'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, 39  where v 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-a method, the tangent stiffness and damping matrix are formed from the partial derivative of t(q, _ q, t), with respect to q and _ 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 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 v co, max '226 rad Á s À1 , and h'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'2:8= ffiffiffiffiffiffi k co p . Hence, the Generalized-a algorithm is better suited for shorter drill string element configurations for Equation (67), while simulating with larger time steps.

Brief comparison of solvers.
This section gives a comparison of applying an RK4 and a G a for a periodic block position trajectory. As RK4 is among one of the most wellknown 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 = (jjp RK4 1 n, b n jj À jjp S n, b n jj) ð85Þ where S 2 fG a1 , . . . , G a5 , RK4 1 , . . . , RK4 5 g, and the superscripts denote the solver configurations for the time steps h 2 f10 À4 , 10 À3 , 10 À2 , 2 Á 10 À2 , 10 À1 g. The block position is given by a periodic function u h = sin (2p(t À 5)=20 + (t À 5) 2 =400) to evaluate the error propagation and convergence. The bit position error between the RK4 and G a with Equation (85) using RK4 with h = 10 À4 as the base case for comparison, is shown in Figure 15.
The square of the input to the block, along with the norm of the bit velocity for RK4 1 and G a1 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 a1 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.
The RMSE is computed from 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 a solver uses p ' = 0:1. The G a 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 a appears as this method is unconditionally stable for larger time steps than RK4.

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. Wellknown 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-a 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-a 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.

Hydraulic model parameters
For the hydraulic model calculations, the initial parameters in this paper are drawn from Gjerstad et al. 17 and are shown in Table 6. Note that the well used to calibrate the model is at bit depths from 4.5 to 6 km, compared with the 3.2 km vertical depth of the well in this paper. As such, the parameters used in the model from our paper serve as an initial guess and require further calibration. Note that RK4 for h = 0:02 and 0:1s is numerically unstable, and did not converge. RTF: real-time factor; RMSE: root-mean-square of the error. Figure 15. Comparing e 2 n for h ∈ f10 À4 ,10 À3 ,10 À2 ,10 À1 g with G α and for RK4 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 dashdotted line h = 10 À1 , respectively.