A multibody dynamics model of contact between the drillstring and the wellbore and the rock penetration process

Abstract

This article models the contact between the drillstring with large slenderness ratio and the extending rigid wellbore based on the multibody dynamics method. The drillstring is modeled as absolute nodal coordinate formulation beams with contact detection points. An algorithm is developed to locate the contact points and calculate the contact forces when the drillstring is sliding and rotating in the wellbore. This provides support force and friction for the moving drillstring and effectively confines it inside the wellbore. A rock penetration model is established based on the rock-breaking velocity equation. The governing equation for the full-hole drilling system including the drillstring and the contact model is established and solved. A rock penetration correction method is proposed to stabilize the computation and to model and simulate the slide drilling process. A field drilling process is modeled and simulated. The simulation result fits the experimental result well.

Introduction

In the past decades, there have been numerous theoretical and simulation investigations of the drilling process on predicting the wellbore and analyzing the drillstring’s motion and deformation.^1–5 Conventional methods include the finite element method (FEM), the differential equation method, and other analytical solutions. Among them, analytical solutions are concise but have over-simplified the stochastic collision between the drillstring and the wellbore; FEM, although playing an important role in drilling process analysis, is not competent to modeling and simulating large deformation, large displacement, and rotation of the drillstring, and the large number of modeling elements greatly impairs the calculation speed.

The multibody dynamics approach has been adopted in drilling system analysis recently. This approach uses the absolute nodal coordinate formulation (ANCF) beam to model the flexible drillstring.^6,7 This method enables modeling and simulation of the large deformation, large displacement, and rotation of the drillstring ⁸ with complex stochastic contact boundaries and slenderness ratio up to 10³. The full-hole multibody dynamics drillstring model is thus developed and the simulation results greatly promote engineering understanding of the drilling process dynamics.

The drillstring–wellbore contact leads to friction and abrasion and may cause severe accidents like stuck-pipe incidents. A lot of studies have been carried out to analyze the contact.^9–12 The bit–rock interaction model is another research focus. Theoretical and experimental work is performed to measure and model the contact force and the rock-breaking process.^13–15 Some statics algorithms are also developed to directly predict the well trajectory in the stable drilling process.^16,17

In this article, section “Modeling of the drillstring–wellbore contact and the bit–rock interaction” introduces a practical drillstring–wellbore contact model and a bit–rock interaction model. Section “Drilling system governing equation and solution” presents a multibody dynamics model of the drilling system, and the system governing equation is presented. When solving the governing equation, the rock penetration process is corrected at the end of each integration step. Simulation results are presented in section “Simulation result.”

Modeling of the drillstring–wellbore contact and the bit–rock interaction

Multibody dynamics model description of the drilling system

A schematic of the drilling system is shown in Figure 1.

Figure 1.

Schematic of drilling system components.

The winch and the top drive control the translational and rotational motion of the drillstring. The drillstring bends and twists inside the wellbore. On the lower end of the drillstring, the bit pushes against the wellbore bottom and penetrates the layer.

In this article, the traveling block, the top drive and the bit are modeled as rigid bodies, and the drill line system is modeled as a spring with variable length and stiffness. The flexible drillstring is modeled with ANCF beam elements.

The drillstring rotates and slides downward in the wellbore and extends the wellbore with a cutting bit at the bottom to the designed direction. The 1000-m-long drillstring, hanging from ground, with a slenderness ratio above 10³, surely would not work without contacting the wellbore. Torque from the top drive or a downhole motor nearby, combined with force from the drillstring, drives the bit to break the rock and thus to extend the wellbore. This makes the wellbore an extending tubular contact surface. Although the drillstring abrades the wellbore and the crustal stress slowly compresses the wellbore, the change in the wellbore section geometry is minor during drilling and is thus dismissed. At the lower end of the wellbore, the well bottom is a moving planar contact surface.

Contact detection points are fixed on the drillstring, so the detection points slide along the wellbore with the drillstring. In coiled tubing drilling, the contact points are evenly distributed on the drill pipe; in other drilling systems, there are tool joints along the drillstring. These tool joints are remarkably thicker than other sections of the drillstring and are the parts to plant detection points. Additional detection points are planted amid the tool joints when necessary. A special detection point against the well bottom is poised on the bit.

Modeling the wellbore

The wellbore is a rigid hole in the ground with a curved centerline and approximately circular cross sections. The wellbore’s centerline is called the well trajectory. A set of wellbore nodes are placed on the well trajectory. In petroleum industry, the well trajectory is described with the measured depth, the inclination, and the azimuth at each node, as shown in Figure 2. At a single point on the well trajectory, the measured depth refers to the arc length of the well trajectory from the starting point of the well, the inclination refers to the angle between the local tangent direction of the trajectory and the vertical downward direction, and the azimuth refers to the angle between the projection of the trajectory on the horizontal plane and the north direction.

Figure 2.

Definition of well trajectory parameters.

Different methods have been developed to transform this set of three parameters to the form of the material coordinate (equals to the measured depth), the position vector, and the tangent direction vector of the well trajectory.^18–20 To describe the wellbore geometry and the contact characteristics, the section radius and the sliding friction coefficient are introduced. Then, the wellbore between two nodes can be interpolated.

In this article, each wellbore node has nine parameters, which are the material coordinate, the position vector, the tangent direction vector of the well trajectory, the radius, and the coefficient of sliding friction of the wellbore, as shown in Figure 3. All these parameters are functions of the material coordinate.

Figure 3.

Wellbore model.

Assume the wellbore has m nodes. The material coordinate at the ith wellbore node is $L_i$ , so that the position vector and the tangent direction vector at the ith wellbore node are

$${\mathbf{x}}_i=[x_i,y_i,z_i{]}^{\mathrm{T}}=\mathbf{x}(L_i),i=1,2,\dots ,m$$

(1)

$${\mathbf{v}}_i=[v_{i,x},v_{i,y},v_{i,z}{]}^{\mathrm{T}}=\mathbf{v}(L_i),i=1,2,\dots ,m$$

(2)

The radius of the wellbore at the ith wellbore node is

$${\rho }_i={\rho }_i(L_i),i=1,2,\dots ,m$$

(3)

The coefficient of sliding friction of the wellbore at the ith wellbore node is

$${\mu }_i={\mu }_i(L_i),i=1,2,\dots ,m$$

(4)

So, the wellbore is modeled using the material coordinate set L, the position vector set X, the tangent direction vector set V, the radius set $\mathbf{\rho }$ , and the coefficient of sliding friction set $\mathbf{\mu }$ at all nodes of the wellbore

$$\mathbf{L}=\{L_1,L_2,\dots ,L_m\}$$

(5)

$$\mathbf{X}(\mathbf{L})=\{{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_m\}$$

(6)

$$\mathbf{V}(\mathbf{L})=\{{\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_m\}$$

(7)

$$\mathbf{\rho }(\mathbf{L})=\{{\rho }_1,{\rho }_2,\dots ,{\rho }_m\}$$

(8)

$$\mathbf{\mu }(\mathbf{L})=\{{\mu }_1,{\mu }_2,\dots ,{\mu }_m\}$$

(9)

In this article, the cubic Hermite interpolation is applied to interpolating the well trajectory between wellbore nodes, as shown in Figure 4.

Figure 4.

Wellbore interpolation.

The dimensionless local material coordinate between the ith and the i + 1th node is s, where

$$s=\frac{2(L-L_i)}{L_{i+1}-L_i}-1,s\in [-1,1]$$

(10)

Thus, the parameters at each node are also the functions of s. The position vector x and the tangent direction vector v at s are

$$\mathbf{x}(s)={\mathbf{x}}_i{\alpha }_i(s)+{\mathbf{x}}_{i+1}{\alpha }_{i+1}(s)+{\mathbf{v}}_i{\beta }_i(s)+{\mathbf{v}}_{i+1}{\beta }_{i+1}(s)$$

(11)

$$\mathbf{v}(s)=\frac{\mathrm{d}\mathbf{x}(s)}{\mathrm{d}s}·\frac{2}{L_{i+1}-L_i}$$

(12)

$${\alpha }_i(s)=\frac{1}{4}{(s-1)}^2(s+2)$$

(13)

$${\alpha }_{i+1}(s)=-\frac{1}{4}{(s+1)}^2(s-2)$$

(14)

$${\beta }_i(s)=\frac{1}{8}{(s-1)}^2(s+1)(L_{i+1}-L_i)$$

(15)

$${\beta }_{i+1}(s)=\frac{1}{8}{(s+1)}^2(s-1)(L_{i+1}-L_i)$$

(16)

Assume the cross section of the wellbore is circular. The radius $\rho$ of the wellbore at s is linearly interpolated

$$\rho (s)=\frac{(1-s)*{\rho }_i}{2}+\frac{(1+s)*{\rho }_{(i+1)}}{2}$$

(17)

where ${\rho }_i$ and ${\rho }_{i+1}$ are the radii of the wellbore at node i and node i + 1, respectively.

Assume the coefficient of sliding friction is continuous along the wellbore and is the linear interpolation of the coefficient of sliding friction at each node. The coefficient of sliding friction $\mu$ at s is

$$\mu (s)=\frac{1-s}{2}{\mu }_i+\frac{1+s}{2}{\mu }_{i+1}$$

(18)

where ${\mu }_i$ and ${\mu }_{i+1}$ are the coefficients of sliding friction at node i and node i + 1, respectively.

Locating the contact point

The contact between the drillstring and the wellbore is to confine the drillstring inside the wellbore. Since the drillstring and the bit are always inside the wellbore, the penetration tolerance of the contact is unlimited. To locate the contact point P of a contact detection point D, we first locate the cross section passing P and then locate P on the cross section.

A bounding box method is used to find at which two nodes the contact point is between, as shown in Figure 5. C is on the well trajectory and marks the wellbore cross section. ${\mathbf{f}}_{\mathrm{n}}$ is the support force.

Figure 5.

Locating contact point.

In petroleum engineering, the curvature of the well trajectory is mostly less than 10°/30 m. The curvature radius of the well trajectory is much larger than the wellbore’s radius. Moreover, the wellbore’s radius is very uniform. Based on these facts, in the neighborhood of P, the wellbore and the drillstring are approximately cylindrical, so that the contact is cylindrical contact. This means D, C, and P are on the line of ${\mathbf{f}}_{\mathrm{n}}$ which is perpendicular to the tangent direction of the wellbore. The cross section of the drillstring at D is in the same plane with the cross section of the wellbore at C.

Set an AABB bounding box for each node. Note $\epsilon$ as the maximum penetration depth of the cylinder into the rigid hole. The edge length of the bounding box for node i is d_i , where

$$d_i=2{\rho }_i+2\epsilon$$

(19)

Then, we can get an extended AABB bounding box enclosing the ith and i + 1th bounding box. Successively check every extended bounding box, as shown in Figure 5. If point D is in the extended bounding box enclosing node i and node i + 1, it is very probable that point C is on the well trajectory between node i and node i + 1.

The precise position s of the contact point C can be located by Newton iteration. Note ${\mathbf{x}}_{\mathrm{P}}$ as the position vector of point P and ${\mathbf{x}}_{\mathrm{C}}$ as the position vector of point C. Note ${\mathbf{e}}_{\tau }$ as the tangent direction of the trajectory at point C. The projection of CP on direction ${\mathbf{e}}_{\tau }$ should be 0

$$f(s)=({\mathbf{x}}_{\mathrm{P}}-{\mathbf{x}}_{\mathrm{C}})·{\mathbf{e}}_{\tau }·\frac{(L_{i+1}-L_i)}{2}=0$$

(20)

Meanwhile, the derivative of this projection for the dimensionless local material coordinate s is non-zero, because

$$\begin{array}{c}\hfill f\prime (s)\hfill \\ \multicolumn{1}{c}{=({\mathbf{x}}_{\mathrm{P}}-{\mathbf{x}}_{\mathrm{C}})·\frac{\mathrm{d}{\mathbf{e}}_{\tau }}{\mathrm{d}s}·\frac{(L_{i+1}-L_i)}{2}-{\mathbf{e}}_{\tau }·{\mathbf{e}}_{\tau }·\frac{{(L_{i+1}-L_i)}^2}{4}}\\ \multicolumn{1}{c}{=({\mathbf{x}}_{\mathrm{P}}-{\mathbf{x}}_{\mathrm{C}})·(x\prime (s),y\prime (s),z\prime (s))-\Vert (x\prime (s),y\prime (s),z\prime (s)){\Vert }^2}\\ \multicolumn{1}{c}{\ne 0(21)}\end{array}$$

(21)

Obviously, $f″(s)$ exists and is continuous. So, the Newton iteration function

$$\phi (s)=s-\frac{f(s)}{f\prime (s)}$$

(22)

is convergent and has second-order convergence rate.

After iteration, if $s<-1$ , the contact detection point is in the prior wellbore sections; if $s>1$ , the contact detection point is in the following wellbore sections.

Drillstring–wellbore contact force calculation

The contact force f is composed of the support force ${\mathbf{f}}_{\mathrm{n}}$ and the friction ${\mathbf{f}}_{\tau }$ , where

$${\mathbf{f}}_{\mathrm{n}}=f_{\mathrm{n}}\mathbf{n},{\mathbf{f}}_{\mathrm{\tau }}=f_{\tau }\mathbf{\tau }$$

(23)

where n is the normal direction of the contact and $\mathbf{\tau }$ is the tangent direction of the contact.

The contact between two cylinders has been modeled with different methods. ²¹ Suppose the support force is a function of the penetration depth and the penetration velocity. Since the penetration tolerance is unlimited, the Hertz contact model is applied, and there is

$$f_n=k{\delta }^e+c\stackrel{·}{\delta }$$

(24)

where k is the contact stiffness coefficient, e is the contact recovery index, c is the contact damping coefficient, $\delta$ is the penetration depth, and $\stackrel{·}{\delta }$ is the penetration velocity. Note ${\gamma }_B$ as the radius of the cylindrical beam. The penetration depth is

$$\delta =\Vert {\mathbf{x}}_{\mathrm{D}}-{\mathbf{x}}_{\mathrm{C}}\Vert +r_{\mathrm{B}}-\rho$$

(25)

Note ${\mathbf{u}}_{\tau }\mathbf{=}\mathbf{u}-(\mathbf{u}·\mathbf{n})\mathbf{n}$ as the tangent velocity of the beam. The direction of the friction is

$$\mathbf{\tau }=\frac{{\mathbf{u}}_{\tau }}{\Vert {\mathbf{u}}_{\tau }\Vert }$$

(26)

The magnitude of the friction is

$$f_{\tau }=\mu (\Vert {\mathbf{u}}_{\tau }\Vert )f_n$$

(27)

where $\mu$ is the friction coefficient and a Coulomb friction model is adopted.

The friction can be decomposed into the axial component and the circumferential component.

The axial component decelerates the axial motion of the drillstring. When pulling up, the friction increases the weight on hook (WOH) and may exceed the maximum WOH. When sliding down, the friction decreases the weight on bit (WOB; axial contact force on the bit from the well bottom) and thus decelerates drilling.

The circumferential component damps the rotation of the drillstring. In rotary drilling process, the circumferential friction increases the torque in the drillstring and may exceed the maximum torque of the drillstring and the top drive. In slide drilling process, when rotating the top drive a small angle to magnitude the bit attitude, the rotation in the drillstring is damped by the friction. The bit attitude control would be difficult. If the top drive rotates too quickly, the torque in the drillstring may accumulate and exceed the maximum torque of the drillstring and the top drive. In some other cases, the sudden release of the friction may lead to thread off, a severe accident in drilling.

Bit–wellbore contact force calculation

The well bottom is a moving contact surface, whose motion is defined in section “Instantaneous rock-breaking velocity calculation.” The track of the moving circular well bottom is the wellbore. Although the bit continuously cuts the well bottom and removes the rock, the bit still hits the static rocks at the well bottom, which is similar to the contact between a rigid ball and a plane. In this article, the bit is modeled as a rigid ball. The contact force ${\mathbf{f}}_{\mathrm{B}-\mathrm{R}}$ on the bit from the well bottom plane is calculated using the Hertz contact model. The WOB is

$$f_{\mathrm{WOB}}=k_{\mathrm{B}-\mathrm{R}}{\delta }_{\mathrm{B}-\mathrm{R}}^e+c_{\mathrm{B}-\mathrm{R}}{\stackrel{·}{\delta }}_{\mathrm{B}-\mathrm{R}}$$

(28)

where $k_{\mathrm{B}-\mathrm{R}}$ is the contact stiffness coefficient, e is the contact recovery index, $c_{\mathrm{B}-\mathrm{R}}$ is the contact damping coefficient, ${\delta }_{\mathrm{B}-\mathrm{R}}$ is the penetration depth, and ${\stackrel{·}{\delta }}_{\mathrm{B}-\mathrm{R}}$ is the penetration velocity. Since the cuttings are carried away by the drilling fluid, the bit continuously hits unbroken rocks, thus the penetration velocity is the velocity of the bit.

The lateral contact force ${\mathbf{f}}_l$ on the bit is the contact force from the side surface of the wellbore, which is expressed in equation (23).

The transformation matrix from the bit body coordinate to the global coordinate is ${\mathbf{A}}_{\mathrm{B}-\mathrm{G}}$ . The contact force is

$${\mathbf{f}}_{\mathrm{B}-\mathrm{R}}=\left[\begin{array}{c}\hfill f_{\mathrm{B}-\mathrm{R},x}\hfill \\ \multicolumn{1}{c}{f_{\mathrm{B}-\mathrm{R},y}}\\ \multicolumn{1}{c}{f_{\mathrm{B}-\mathrm{R},z}}\end{array}\right]={\mathbf{A}}_{\mathrm{B}-\mathrm{G}}\left[\begin{array}{c}\hfill f_{\mathrm{WOB}}\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right]+{\mathbf{f}}_l$$

(29)

The axial contact torque $T_{\mathrm{BIT}}$ on the bit is a function of the WOB, the bit radius $r_{\mathrm{BIT}}$ , and the axial angular velocity ${\omega }_{\mathrm{BIT},\xi }$ of the bit ¹

$$T_{\mathrm{BIT}}={\mu }_{\mathrm{BIT}}{f}_{\mathrm{WOB}}[\tanh ({\omega }_{\mathrm{BIT},\xi })+\frac{b_1{\omega }_{\mathrm{BIT},\xi }}{b_2{\Omega }^2+1}]·\frac{2}{3}{r}_{\mathrm{BIT}}$$

(30)

where ${\mu }_{\mathrm{BIT}}$ is the friction coefficient and $b_1$ and $b_2$ are two parameters.

Instantaneous rock-breaking velocity calculation

The well bottom moving velocity is commonly called the rock penetration velocity. In this article, to be distinguished from the contact penetration velocity, the well bottom moving velocity is called the rock-breaking velocity ${\mathbf{u}}_{\mathrm{B}}$ . ${\mathbf{u}}_{\mathrm{B}}$ is a function of the force and torque on the bit, the bit rotary speed, the bit attitude, the bit parameters, and the layer parameters.^7,13–15 The force and torque on the bit and the bit attitude are functions of the bit load condition and motion state. The bit rotary speed, the bit parameters, and the layer parameters are parameters. For tri-cone roller (TCR) bits, the impact of the bit rotary speed and the torque on bit on ${\mathbf{u}}_{\mathrm{B}}$ is minor. According to Yang, ⁷ in an anisotropic layer, the instantaneous rock-breaking velocity can be simplified

$${\mathbf{u}}_{\mathrm{B}}=\frac{-1}{K}{\mathbf{A}}_{\mathrm{B}-\mathrm{G}}\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \hfill & \hfill \hfill \\ \multicolumn{1}{c}{}& \hfill c_1\hfill & \hfill \hfill \\ \multicolumn{1}{c}{}& \hfill \hfill & \hfill c_2\hfill \end{array}\right]({\mathbf{A}}_{\mathrm{G}-\mathrm{B}}\left[\begin{array}{c}\hfill f_x\hfill \\ \multicolumn{1}{c}{f_y}\\ \multicolumn{1}{c}{f_z}\end{array}\right]-\left[\begin{array}{c}\hfill f_{\mathrm{c},\xi }\hfill \\ \multicolumn{1}{c}{f_{\mathrm{c},\eta }}\\ \multicolumn{1}{c}{f_{\mathrm{c},\zeta }}\end{array}\right])$$

(31)

where K is a coefficient describing the rock drillability, ${\mathbf{A}}_{\mathrm{G}-\mathrm{B}}$ is the transformation matrix from the global coordinate to the bit body coordinate, and $f_{\mathrm{c},\xi }$ is the critical load in the axial direction $\xi$ . $f_{\mathrm{c},\eta }$ and $f_{\mathrm{c},\zeta }$ are the critical loads in the lateral direction $\eta$ and $\zeta$ , respectively. Although the bit is designed to drill ahead, it also slowly breaks the rock on its side, which is modeled by the velocity coefficients $c_1$ and $c_2$ .

Equation (31) is the three-dimensional (3D) drilling velocity equation for normal drilling process. This formulation is also valid for polycrystalline diamond compact (PDC) bits when the bit rotary speed is stable.

The magnitude of ${\mathbf{u}}_{\mathrm{B}}$ is called the rate of penetration (ROP). In the rock penetration process, there are three phases. ¹⁴ In phase I, the WOB is below the threshold and the bit grinds the rock. The ROP is very low and can be dismissed. In phase II, the WOB is large enough to break the rocks, and the ROP increases in proportion to the WOB. In phase III, the WOB is very large, and the situation becomes complicated. The phase II rock penetration process is the normal drilling process which is modeled by equation (31).

The direct application of the instantaneous rock-breaking velocity may lead to calculation difficulties and very bad wellbore prediction quality. A two-stage correction method will be presented in section “Drilling system governing equation and solution.”

Drilling system governing equation and solution

System governing equation of the drilling system

The drilling system shown in Figure 1 is modeled as a rigid-flexible coupling system based on multibody dynamics method, as shown in Figure 6.

Figure 6.

Bodies and constraints.

The crown block, the traveling block, the top drive stator, the top drive rotor, and the bit are modeled as rigid bodies. The crown block is fixed on the derrick. This constraint is modeled as a fixed joint to the ground.

The drill line between the crown block and the traveling block is modeled as a variable-length variable-stiffness extension spring. Suppose there are n pulleys in the traveling block. The length of the fast line (the drill line section connecting the winch) is $l_1$ . The length of the deadline (the drill line section connecting the derrick at the deadline anchor) is $l_2$ . The undeformed length of the drill line between the crown block and the traveling block is $l_3$ . The stiffness coefficient of the spring is

$$k_{\mathrm{s}}=\frac{1}{\frac{1}{4n^2}(\frac{1}{k_1}+\frac{1}{k_2})+\frac{1}{2{\mathrm{nk}}_3}},k_i=\frac{\mathrm{EA}}{l_i},i=1,2,3$$

(32)

where E is Young’s modulus of the drill line, and A is the cross-sectional area of the drill line.

The motion of the winch leads to the length change of the drill line between the crown block and the traveling block, which is modeled as the length change of $l_3$ .

The traveling block and the top drive stator are connected by a fixed joint. The top drive travels along a vertical rail on the derrick and is modeled as constrained by a translational joint with the ground. A revolute joint connects the top drive stator and the top drive rotor. The rotation between the top drive rotor and the stator is modeled as an angular velocity constraint.

To measure the output torque of the top drive, a very stiff torsional spring is added between the top drive rotor and the drillstring along the axis. The motion in the other five dimensions is constrained. The bit is fixed to the bottom of the drillstring.

There are five rigid bodies and six constraints in the drilling system. Based on the ANCF, the coordinate in this article is defined in the ground coordinate system.

The attitude of a rigid body or a beam node is described using the Euler rotation vector $\Phi =[{\phi }_1,{\phi }_2,{\phi }_3{]}^{\mathrm{T}}=\varphi \mathbf{e}$ , where e is the rotation axis and $\varphi$ is the rotation angle. The angular velocity of the rigid body or the beam node in the body coordinate system is

$$\mathbf{\omega }={\mathbf{H}}^{\mathrm{T}}\stackrel{·}{\mathbf{\Phi }}$$

(33)

$$\mathbf{H}=\mathbf{I}+h_2\tilde{\mathbf{\Phi }}+h_3\tilde{\mathbf{\Phi }}\tilde{\mathbf{\Phi }}=h_1\mathbf{I}+h_2\tilde{\mathbf{\Phi }}+h_3{\mathbf{\Phi }\mathbf{\Phi }}^{\mathrm{T}}$$

(34)

$$h_1=\frac{\sin \varphi }{\varphi },h_2=\frac{\varphi -\cos \varphi }{{\varphi }^2},h_3=\frac{\varphi -\sin \varphi }{{\varphi }^3}$$

(35)

The position coordinate of the ith rigid body is ${\mathbf{r}}_i$ . The general coordinate of the ith rigid body is ${\mathbf{q}}_i=[{\mathbf{r}}_i^{\mathrm{T}},{\mathbf{\Phi }}_i^{\mathrm{T}}{]}^{\mathrm{T}}$ . Based on the Lagrange equation, the governing equation of the ith rigid body is

$${\mathbf{M}}_i{\stackrel{··}{\mathbf{q}}}_i-{\mathbf{P}}_i-{\mathbf{Q}}_i+{\displaystyle \sum _{j=1}^6}\frac{\partial {\mathbf{C}}_j}{{\mathbf{q}}_i}{\mathbf{\lambda }}_j=0,i=1,2,3,4,5$$

(36)

$${\mathbf{C}}_j(\mathbf{q},\stackrel{·}{\mathbf{q}},t)=0,j=1,2,\dots ,6$$

(37)

$${\mathbf{M}}_i=\left[\begin{array}{cc}\hfill m_i{\mathbf{I}}_{3\times 3}\hfill & \hfill \hfill \\ \multicolumn{1}{c}{}& \hfill {\mathbf{H}}_i{\mathbf{J}}_i{\mathbf{H}}_i^{\mathrm{T}}\hfill \end{array}\right]$$

(38)

$${\mathbf{P}}_i=\left[\begin{array}{cc}\hfill 0_{3\times 3}\hfill & \hfill \hfill \\ \multicolumn{1}{c}{}& \hfill \frac{\mathrm{d}}{\mathrm{d}t}({\mathbf{H}}_i{\mathbf{J}}_i{\mathbf{H}}_i^{\mathrm{T}})-\frac{1}{2}{\stackrel{·}{\mathbf{\Phi }}}_i^{\mathrm{T}}\frac{\partial ({\mathbf{H}}_i{\mathbf{J}}_i{\mathbf{H}}_i^{\mathrm{T}})}{\partial {\mathbf{\Phi }}_i}\hfill \end{array}\right]{\stackrel{·}{\mathbf{q}}}_i$$

(39)

where ${\mathbf{M}}_i$ is the generalized mass matrix. $m_i$ is the mass. ${\mathbf{J}}_i$ is the rotational inertia matrix. ${\mathbf{Q}}_i$ is the external generalized force. ${\mathbf{C}}_j$ is the jth constraint of the bodies, and ${\mathbf{\lambda }}_j$ is the respective Lagrange multiplier.

An ANCF beam element is applied to modeling the drillstring. The cross section is assumed to be rigid, while shear is considered and the normal direction of the section is not along the centerline of the beam element. Suppose the beam element has n nodes, each with six generalized coordinates. The generalized coordinates of the nodes and the beam element are

$${\overline{\mathbf{q}}}_i=[{\mathbf{r}}_i^{\mathrm{T}},{\mathbf{\Phi }}_i^{\mathrm{T}}{]}^{\mathrm{T}},i=1,2,\dots ,n$$

(40)

$$\overline{\mathbf{q}}=[{\mathbf{r}}_1^{\mathrm{T}},{\mathbf{\Phi }}_1^{\mathrm{T}},{\mathbf{r}}_2^{\mathrm{T}},{\mathbf{\Phi }}_2^{\mathrm{T}},\dots ,{\mathbf{r}}_n^{\mathrm{T}},{\mathbf{\Phi }}_n^{\mathrm{T}}{]}^{\mathrm{T}}$$

(41)

$$\overline{\mathbf{r}}=[{\mathbf{r}}_1^{\mathrm{T}},0_{3\times 3},{\mathbf{r}}_2^{\mathrm{T}},0_{3\times 3},\dots ,{\mathbf{r}}_n^{\mathrm{T}},0_{3\times 3}{]}^{\mathrm{T}}$$

(42)

where ${\mathbf{r}}_i$ and ${\mathbf{\Phi }}_i$ are the position vector and the Euler rotation vector of the ith node of the beam element, respectively.

The generalized coordinate at l in the beam element is interpolated as

$$\overline{\mathbf{q}}={\displaystyle \sum _{i=1}^n}{\mathbf{N}}_i(l){\overline{\mathbf{q}}}_i$$

(43)

where l is the arc length coordinate in the beam element and ${\mathbf{N}}_i(l)$ is the shape function.

Based on the rotation vector $\mathbf{\Phi }$ , a Frenet frame $[{\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3]$ is set on the axis of the beam element. ${\mathbf{e}}_1$ points to the normal direction of the beam cross section. ${\mathbf{e}}_1$ , ${\mathbf{e}}_2$ , and ${\mathbf{e}}_3$ form a right-handed orthogonal Cartesian coordinate system. The normal strain and the two engineering shear strains along the centerline are

$$\begin{array}{cc}\multicolumn{1}{c}{{\epsilon }_{11}}& =\frac{1}{2}({\left(\frac{\partial \mathbf{r}}{\partial s}\right)}^{\mathrm{T}}\left(\frac{\partial \mathbf{r}}{\partial s}\right)-1)\hfill \\ \multicolumn{1}{c}{{\gamma }_{12}}& ={\left(\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}s}\right)}^{\mathrm{T}}{\mathbf{e}}_2,{\gamma }_{13}={\left(\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}s}\right)}^{\mathrm{T}}{\mathbf{e}}_3\hfill \end{array}$$

(44)

The curvature along the beam is

$$\mathbf{\kappa }=[{\kappa }_1,{\kappa }_2,{\kappa }_3{]}^{\mathrm{T}}={\mathbf{H}}^{\mathrm{T}}\frac{\mathrm{d}\mathbf{\Phi }}{\mathrm{d}s}$$

(45)

According to the d’Alembert’s principle, the total virtual work of the inertia force, the elastic force, the constraint force, and the external force (the contact force from the wellbore) are 0

$$\delta W_{\mathrm{inertia}}+\delta W_{\mathrm{elastic}}+\delta W_{\mathrm{external}}+\delta W_{\mathrm{constraint}}=0$$

(46)

The governing equation of the beam element is

$${\mathbf{f}}_{\mathrm{inertia}}+{\mathbf{f}}_{\mathrm{elastic}}+{\mathbf{f}}_{\mathrm{external}}+{\mathbf{f}}_{\mathrm{constraint}}=0$$

(47)

$${\mathbf{f}}_{\mathrm{inertia}}=-\overline{\mathbf{M}}\stackrel{··}{\overline{}}\mathbf{q}-\overline{\mathbf{W}}\stackrel{·}{\overline{}}\mathbf{q}=-{\int }_0^{l_0}[\rho A{\mathbf{N}}^{\mathrm{T}}\stackrel{··}{\overline{}}\mathbf{r}+{\mathbf{N}}^{\mathrm{T}}\mathbf{HJ}\stackrel{·}{\overline{}}\omega ]\mathrm{d}l-{\int }_0^{l_0}{\mathbf{N}}^{\mathrm{T}}\mathbf{H}{\overline{\omega }}^{\mathrm{T}}\mathbf{J}\overline{\omega }\mathrm{d}l$$

(48)

$${\mathbf{f}}_{\mathrm{external}}=\overline{\mathbf{Q}}$$

(49)

$${\mathbf{f}}_{\mathrm{constraint}}=-{\displaystyle \sum _{j=1}^6}{\mathbf{C}}_{j,\overline{\mathbf{q}}}{\mathbf{\lambda }}_j$$

(50)

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathbf{f}}_{\mathrm{elastic}}}& =-{\left(\frac{\partial U_{\mathrm{e}}}{\partial \overline{\mathbf{q}}}\right)}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{}& ={\displaystyle \underset{0}{\overset{l_0}{\int }}}\mathrm{kG}({\gamma }_{12}{\left(\frac{\partial {\gamma }_{12}}{\partial \overline{\mathbf{q}}}\right)}^{\mathrm{T}}+{\gamma }_{13}{\left(\frac{\partial {\gamma }_{13}}{\partial \overline{\mathbf{q}}}\right)}^{\mathrm{T}})\mathrm{d}l\hfill \\ \multicolumn{1}{c}{}& -{\displaystyle \underset{0}{\overset{l_0}{\int }}}[\mathrm{EA}{\epsilon }_{11}{\left(\frac{\partial {\epsilon }_{11}}{\partial \overline{\mathbf{q}}}\right)}^{\mathrm{T}}+{\mathrm{EJ}}_T{\kappa }_1{\left(\frac{\partial {\kappa }_1}{\partial \overline{\mathbf{q}}}\right)}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{}& +{\mathrm{EJ}}_{22}{\kappa }_2{\left(\frac{\partial {\kappa }_2}{\partial \overline{\mathbf{q}}}\right)}^{\mathrm{T}}+{\mathrm{EJ}}_{33}{\kappa }_3{\left(\frac{\partial {\kappa }_3}{\partial \overline{\mathbf{q}}}\right)}^{\mathrm{T}}]\mathrm{d}l\hfill \end{array}$$

(51)

where $l_0$ is the beam length, E is Young’s modulus, G is the shear modulus, and k is a shear coefficient relevant to the shape of the cross section. $\overline{\mathbf{Q}}$ is the external generalized force, including the contact force, the torsional spring, and the gravity force. $\overline{\mathbf{\omega }}$ is the angular velocity in the global coordinate system

$$\overline{\mathbf{\omega }}=[0_{3\times 3},{\overline{\mathbf{\omega }}}_1^{\mathrm{T}},0_{3\times 3},{\overline{\mathbf{\omega }}}_2^{\mathrm{T}},\dots ,0_{3\times 3},{\overline{\mathbf{\omega }}}_n^{\mathrm{T}}{]}^{\mathrm{T}}$$

(52)

Note

$$\overline{\mathbf{P}}=-\overline{\mathbf{W}}\stackrel{·}{\overline{\mathbf{q}}}-{\left(\frac{\partial U_{\mathrm{e}}}{\partial \overline{\mathbf{q}}}\right)}^{\mathrm{T}}$$

(53)

Equation (47) becomes

$$\overline{\mathbf{M}}\stackrel{··}{\overline{\mathbf{q}}}-\overline{\mathbf{P}}-\overline{\mathbf{Q}}+{\displaystyle \sum _{j=1}^6}\frac{\partial {\mathbf{C}}_j}{\partial \overline{\mathbf{q}}}{\mathbf{\lambda }}_j=0$$

(54)

According to equation (31), the rock penetration velocity ${\mathbf{u}}_{\mathrm{B}}$ is a function of the bit–well contact force ${\mathbf{f}}_{\mathrm{B}-\mathrm{R}}$ and the transformation matrixes ${\mathbf{A}}_{\mathrm{G}-\mathrm{B}}$ and ${\mathbf{A}}_{\mathrm{B}-\mathrm{G}}$ . ${\mathbf{A}}_{\mathrm{G}-\mathrm{B}}$ and ${\mathbf{A}}_{\mathrm{B}-\mathrm{G}}$ are functions of the bit generalized coordinate ${\mathbf{q}}_5$ .

The time derivative of the position and the material coordinate of the well bottom node position are

$$\stackrel{·}{\mathbf{x}}={\mathbf{u}}_{\mathrm{B}}({\mathbf{f}}_{\mathrm{B}-\mathrm{R}},{\mathbf{q}}_5)={\mathbf{u}}_{\mathrm{B}}({\mathbf{q}}_5,{\stackrel{·}{\mathbf{q}}}_5,\mathbf{x})$$

(55)

$$\stackrel{·}{L}=\Vert {\mathbf{u}}_{\mathrm{B}}({\mathbf{f}}_{\mathrm{B}-\mathrm{R}},{\mathbf{q}}_5)\Vert =\Vert {\mathbf{u}}_{\mathrm{B}}({\mathbf{q}}_5,{\stackrel{·}{\mathbf{q}}}_5,\mathbf{x})\Vert =\Vert \stackrel{·}{\mathbf{x}}\Vert$$

(56)

The tangent direction vector of the bottom node is pointing to the direction of the rock-breaking velocity

$$\mathbf{v}=\frac{{\mathbf{u}}_{\mathrm{B}}({\mathbf{f}}_{\mathrm{B}-\mathrm{R}},{\mathbf{q}}_5)}{\Vert {\mathbf{u}}_{\mathrm{B}}({\mathbf{f}}_{\mathrm{B}-\mathrm{R}},{\mathbf{q}}_5)\Vert }=\frac{{\mathbf{u}}_{\mathrm{B}}({\mathbf{q}}_5,{\stackrel{·}{\mathbf{q}}}_5,\mathbf{x})}{\Vert {\mathbf{u}}_{\mathrm{B}}({\mathbf{q}}_5,{\stackrel{·}{\mathbf{q}}}_5,\mathbf{x})\Vert }=\frac{\stackrel{·}{\mathbf{x}}}{\Vert \stackrel{·}{\mathbf{x}}\Vert }$$

(57)

Thus, x, L, and v can be acquired from $\stackrel{·}{\mathbf{x}}$ . Substituting equation (31) into equation (55), there is

$$\stackrel{·}{\mathbf{x}}-\frac{1}{K}{\mathbf{A}}_{\mathrm{B}-\mathrm{G}}\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \hfill & \hfill \hfill \\ \multicolumn{1}{c}{}& \hfill c_1\hfill & \hfill \hfill \\ \multicolumn{1}{c}{}& \hfill \hfill & \hfill c_2\hfill \end{array}\right]({\mathbf{A}}_{\mathrm{G}-\mathrm{B}}{\mathbf{f}}_{\mathrm{B}-\mathrm{R}}-\left[\begin{array}{c}\hfill f_{\mathrm{c},\xi }\hfill \\ \multicolumn{1}{c}{f_{\mathrm{c},\eta }}\\ \multicolumn{1}{c}{f_{\mathrm{c},\zeta }}\end{array}\right])=0$$

(58)

The radius $\rho$ and friction coefficient $\mu$ of the new well bottom node are preset according to the bit type and the layer characteristics.

The system governing equation describing this highly coupled process is established using the Lagrange equations. The generalized coordinate is

$$\mathbf{q}=[{\mathbf{q}}_1^{\mathrm{T}},{\mathbf{q}}_2^{\mathrm{T}},\dots ,{\mathbf{q}}_5^{\mathrm{T}},{\overline{\mathbf{q}}}^{\mathrm{T}}{]}^{\mathrm{T}}$$

(59)

The generalized mass matrix is

$$\mathbf{M}=\left[\begin{array}{ccccc}\hfill {\mathbf{M}}_1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \multicolumn{1}{c}{}& \hfill {\mathbf{M}}_2\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \multicolumn{1}{c}{}& \hfill \hfill & \hfill \ddots \hfill & \hfill \hfill & \hfill \hfill \\ \multicolumn{1}{c}{}& \hfill \hfill & \hfill \hfill & \hfill {\mathbf{M}}_5\hfill & \hfill \hfill \\ \multicolumn{1}{c}{}& \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \overline{\mathbf{M}}\hfill \end{array}\right]$$

(60)

The generalized forces are

$$\mathbf{P}=[{\mathbf{P}}_1^{\mathrm{T}},{\mathbf{P}}_2^{\mathrm{T}},\dots ,{\mathbf{P}}_5^{\mathrm{T}},{\overline{\mathbf{P}}}^{\mathrm{T}}{]}^{\mathrm{T}}$$

(61)

$$\mathbf{Q}=[{\mathbf{Q}}_1^{\mathrm{T}},{\mathbf{Q}}_2^{\mathrm{T}},\dots ,{\mathbf{Q}}_5^{\mathrm{T}},{\overline{\mathbf{Q}}}^{\mathrm{T}}{]}^{\mathrm{T}}$$

(62)

The governing equation of the drilling system is

$$\mathbf{M}\stackrel{··}{\mathbf{q}}-\mathbf{P}-\mathbf{Q}+{\displaystyle \sum _{j=1}^6}\frac{\partial {\mathbf{C}}_j}{\partial \mathbf{q}}{\mathbf{\lambda }}_j=0$$

(63a)

$${\mathbf{C}}_j(\mathbf{q},\stackrel{·}{\mathbf{q}},t)=0,j=1,2,\dots ,6$$

(63b)

$$\begin{array}{c}\hfill \stackrel{·}{\mathbf{x}}={\mathbf{u}}_{\mathrm{B}}({\mathbf{f}}_{\mathrm{B}-\mathrm{R}},{\mathbf{q}}_5)\hfill \\ \multicolumn{1}{c}{\stackrel{·}{L}=\Vert \stackrel{·}{\mathbf{x}}\Vert }\\ \multicolumn{1}{c}{\mathbf{v}=\frac{\stackrel{·}{\mathbf{x}}}{\Vert \stackrel{·}{\mathbf{x}}\Vert }}\end{array}\}$$

(63c)

Equation (63) is a differential algebraic equation (DAE). Equation (63a) is the dynamics equation of the bodies and is a combination of equations (36) and (54). Equation (63b) is the constraint equation for the rigid bodies and is from equation (37). Equation (63c) is the rock penetration equation and is a combination of equations (55)–(57). Thus, the rock penetration process is introduced into the system governing equation. The stiffness problem is converted into the DAE solving problem. A stiff DAE solver is essential to solving equation (63).

Solving the governing equation

Since the rock penetration process is very slow, the contact condition of the drillstring is almost unchanged before and after a single integration step. The solution sequence diagram is shown in Figure 7. At the kth integration step, the time is $t^k$ . First, the system generalized coordinate ${\mathbf{q}}^k$ is calculated from ${\mathbf{q}}^{k-1}$ and the history of the well bottom node parameters. Then, ${\stackrel{·}{\mathbf{q}}}^k$ and ${\stackrel{·}{\mathbf{x}}}^k$ are calculated from the history of ${\mathbf{q}}^k$ . Finally, the well bottom node parameters $L^k,{\mathbf{x}}^k,{\mathbf{v}}^k$ are calculated from ${\stackrel{·}{\mathbf{x}}}^k$ .

Figure 7.

Solution sequence diagram.

Suppose the initial wellbore has m nodes. A new wellbore node is added to the well bottom at the end of every integration step, and

$${\mathbf{x}}_{m+k}={\mathbf{x}}^k,{\mathbf{v}}_{m+k}={\mathbf{v}}^k,L_{m+k}=L^k$$

(64)

The governing equation at the kth integration step is

$$\begin{array}{cc}\multicolumn{1}{c}{}& \mathbf{M}({\mathbf{q}}^k){\stackrel{··}{\mathbf{q}}}^k-{\mathbf{P}}_i({\mathbf{q}}^k,{\stackrel{·}{\mathbf{q}}}^k)+{\displaystyle \sum _{j=1}^6}\frac{\partial {\mathbf{C}}_j^k}{\partial \mathbf{q}}{\mathbf{\lambda }}_j^k\hfill \\ \multicolumn{1}{c}{}& -{\mathbf{Q}}_i({\mathbf{L}}^{k-1},{\mathbf{X}}^{k-1},{\mathbf{V}}^{k-1},{\mathbf{\rho }}^{k-1},{\mathbf{\mu }}^{k-1},{\mathbf{q}}^k,{\stackrel{·}{\mathbf{q}}}^k)=0\hfill \end{array}$$

(65a)

$${\mathbf{C}}_j({\mathbf{q}}^k,{\stackrel{·}{\mathbf{q}}}^k,t^k)=0,j=1,2,\dots ,6$$

(65b)

$$\begin{array}{cc}\multicolumn{1}{c}{}& L_{m+k}=L_{m+k-1}+\Vert {\mathbf{u}}_{\mathrm{B}}({\mathbf{f}}_{\mathrm{B}-\mathrm{R}}^k,{\mathbf{q}}_5^k)\Vert ·(t^k-t^{k-1})\hfill \\ \multicolumn{1}{c}{}& {\mathbf{x}}_{m+k}={\mathbf{x}}_{m+k-1}+{\mathbf{u}}_{\mathrm{B}}({\mathbf{f}}_{\mathrm{B}-\mathrm{R}}^k,{\mathbf{q}}_5^k)·(t^k-t^{k-1})\hfill \\ \multicolumn{1}{c}{}& {\mathbf{v}}_{m+k}=\frac{{\mathbf{u}}_{\mathrm{B}}({\mathbf{f}}_{\mathrm{B}-\mathrm{R}}^k,{\mathbf{q}}_5^k)}{\Vert {\mathbf{u}}_{\mathrm{B}}({\mathbf{f}}_{\mathrm{B}-\mathrm{R}}^k,{\mathbf{q}}_5^k)\Vert }\hfill \end{array}\}$$

(65c)

The calculation process in each step is illustrated in Figure 8. In each integration step, the load condition and the motion state of the drillstring is first calculated in equations (65a) and (65b). Then, equation (65c) is calculated and the rock penetration result is acquired.

Figure 8.

One integration step of solving the system governing equation.

In this article, the backward differentiation formula method is applied to solving the DAE. Based on this calculation strategy, the typical time consumption of the rock penetration calculation is about 0.1% of the total calculation time.

Rock penetration correction

During the drilling process, the bit vibrates fiercely, and the direction of the instantaneous rock-breaking velocity calculated in a single integration step might be in any direction. The calculated wellbore will be highly curved, as shown in Figure 9(a), and should be corrected. A two-stage rock penetration correction strategy is applied at the end of each integration step, as shown in Figure 10.

Figure 9.

Rock penetration process calculation with two-stage correction.

Figure 10.

Stage II correction: (a) well trajectory before stage II correction and (b) well trajectory after stage II correction.

Stage I

At each integration step, the average rock-breaking velocity ${\overline{\mathbf{u}}}_{\mathrm{B}}$ of the current integration step and the last several steps is applied instead of ${\mathbf{u}}_{\mathrm{B}}$ . Thus, the new node position vector is

$${\mathbf{x}}_{m+k+1}={\mathbf{x}}_{m+k}+{\overline{\mathbf{u}}}_{\mathrm{B}}·(t^{k+1}-t^k)$$

(66)

with the bit direction as the new node tangent direction

$${\mathbf{v}}_{m+k+1}=\mathbf{v}({\mathbf{q}}_{\mathrm{bit}})$$

(67)

Stage II

After every N integration steps, set the node tangent direction to be along the vector from the first node of the N nodes to the last node

$${\mathbf{v}}_{m+k+N}=\frac{{\mathbf{x}}_{m+k+N}-{\mathbf{x}}_{m+k}}{\Vert {\mathbf{x}}_{m+k+N}-{\mathbf{x}}_{m+k}\Vert }$$

(68)

Then, delete the first N − 1 nodes generated in the last N steps and reserve the Nth node, as shown in Figure 10.

The two-stage correction of the wellbore is essential for slide drilling system. The first stage stabilizes the calculation, and the second stage corrects the average rock-breaking velocity.

In stage I, the tangent direction of the new bottom node calculated in each step could only be the direction of the bit, as illustrated in Figure 11. If the direction of ${\mathbf{u}}_{\mathrm{B}}$ is adopted directly, in the first integration step, the angle between the bottom node tangent direction and the bit axis direction can be very large; in the second integration step, the bit would vibrate to the other direction and punch the wellbore heavily. However, in the second integration step, the bit would hit the side of the wellbore instead of the bottom, thus the rock is not broken in this step according to the rock penetration algorithm. In some cases, the colliding force could be very large that the calculation will diverge.

Figure 11.

Initial contact between the bit and the wellbore bottom without stage I correction.

In stage II, the wellbore passes a low-pass filter, so that the zigzag in the well trajectory is diminished, as shown in Figure 9. A simulation result of the stage II correction is shown in section “Wellbore correction effectiveness.” In practical drilling process, there is a procedure after drilling called reaming, which makes the wellbore smoother and straighter. After reaming, the high-frequency oscillation in the well trajectory tangent direction is eliminated. So, the low-pass filter helps better approximate the real well trajectory. Moreover, in the compound drilling process (the process uses a slide drilling system but rotates the top drive to get a straight wellbore), none of the bit direction, ${\mathbf{u}}_{\mathrm{B}}$ or ${\overline{\mathbf{u}}}_{\mathrm{B}}$ , is pointing to the average rock penetration direction. The well trajectory will become a twisted helix without stage II correction, as shown in section “Wellbore correction effectiveness.”

In addition, stage II correction reduces the wellbore nodes and thus reduces the storage cost and accelerates the contact point locating process.

Simulation result

Test oil well description

A full-hole multibody dynamics model of a slide drilling system in Bohai Sea, like the system shown in Figure 1, is established. The bottom hole assembly includes a PDC bit, a mud motor with a 1.25° bent housing, a float valve, a stabilizer, a measurement while drilling (MWD) system, and 14 heavyweight drill collars.

A rotary drilling process begins when the measured depth reaches 730 m and ends at the measured depth of 760 m. The initial inclination is 51° and the initial azimuth is 297°. The bit is initially off the bottom. At 250 s, the bit hits the well bottom and begins to break the rocks. The ROP is 32.7 m/h (9.08 mm/s). The WOH remains around 34 T. The angular velocity of the top drive is 53 r/min (5.55 rad/s). The torque of the top drive remains around $8\mathrm{kN}\mathrm{m}$ . In such a drilling process, the newly generated well trajectory is a helix with very small radius and will become a straight line after reaming.

This compound drilling process is modeled and simulated, and the result is compared to the experimental result. The friction coefficient and the layer parameters of the wellbore are determined so as to make sure the simulated WOH and top drive torque are similar to the experimental results.

Wellbore correction effectiveness

In Figure 12, the helical red line is the original prediction of the well trajectory without stage II correction. Stage II correction makes the wellbore smoother and produces a new helix of much longer thread pitch, which is the blue line with pentagons.

Figure 12.

The simulated well trajectory before and after stage II correction.

Well trajectory prediction quality examination

The MWD system records the inclination and azimuth while drilling and offers data to compare with the simulation result. Figures 13 and 14 show the horizontal and vertical projections of the well trajectory, respectively. The simulated well trajectory fits well with the experiment.

Figure 13.

Horizontal projection of the well trajectory.

Figure 14.

Vertical projection of the well trajectory.

Figures 15 and 16 show the inclination and the azimuth of the well trajectory, respectively. The wellbore in compound drilling process has a helical well trajectory. However, reaming enlarges the wellbore and turns the experiment wellbore from a thin helical hole to a thick straight hole. Additionally, the sampling rate of the MWD system is lower than the simulation result and reveals less information of the fluctuation of the well trajectory tangent direction. So, the experiment inclination and azimuth are nearly straight lines (the red lines with pentagons), while the simulation results are oscillating (the blue lines). The average inclination and the azimuth in simulation well fit the experiment data.

Figure 15.

Inclination of the wellbore.

Figure 16.

Azimuth of the wellbore.

Conclusion

In this article, the contact between the drillstring and the wellbore and the rock penetration process are modeled under the multibody dynamics method framework to model the drilling process. A rigid-flexible coupled multibody dynamics model of the drilling system is applied. The oil derrick components are modeled with rigid bodies, the drillstring is modeled with ANCF beams, and the bit is also considered a rigid body. Geometry and friction coefficient of the wellbore are interpolated from nodal values. Hertz contact forces between the contact points on the drillstring and the wellbore boundary provide support force and friction to the drillstring and effectively confine the drillstring inside the wellbore. A practical algorithm is developed to locate the contact points with respect to the Hermite nodes while the drillstring is carrying the contact points and sliding in the wellbore. On each step of the multibody system equation integration, the wellbore is extended at the bottom according to the 3D drilling velocity equation. A wellbore trajectory correction method is developed out of the detailed examination of the bit–bottom contacts within several integration steps.

A compound drilling process of 30 m in a 730-m-depth oil well in Bohai Sea is simulated. The average WOB, top torque, as well as the newly drilled well trajectory in the simulation results show very good correspondence with the collected field drilling data. The proposed method proves effective in drilling process simulation.