A multibody dynamic model of the drilling system with drilling fluid

Abstract

This article is intended to present a multibody dynamic model of the drilling system, consisting of drillstring and drilling fluid. The drillstring is a complex rigid–flexible coupling system, including rigid bodies, Euler–Bernoulli beam elements, constraints and dynamic loads, and its dynamic model is established using the absolute nodal coordinate formulation. The drilling fluid, composed of internal, annulus, and under-bit fluids, is modeled as one-dimensional compressible fluid; the relative flow of the drilling fluid is modeled using the Arbitrary Lagrangian–Eulerian description; the force of the drillstring acting on the drilling fluid is introduced through the drilling fluid transport motion; meanwhile, the reaction force acting on the drillstring is taken as an external load. The contact between the drillstring and drilling fluid is simulated based on Hertz contact theory, and the rock penetration model is built based on the rock-breaking velocity equation. Based on this model, the coupled vibration of the drillstring and the effects of the drilling fluid flow rate and density on the drilling process are investigated through several examples.

Introduction

In drilling engineering, due to the interaction of the drillstring with the drilling fluid, as well as the dynamic contact between the drillstring and the wellbore, the drilling system is a complex nonlinear fluid–solid coupling system. Therefore, it is necessary to establish a model with high computational accuracy and efficiency that can simulate the dynamic behavior of the drillstring sliding and rotating in the wellbore filled up with drilling fluid.

Several methods have been widely implemented to investigate the dynamic characteristics of the drillstring system, such as the finite element method (FEM),^1–4 the differential equation method, ⁵ and other analytical solution.

Although these methods are relatively mature, there are still some deficiencies in practical applications due to the low computational efficiency or the over-simplification of the system. For example, the FEM has some difficulty in simulating the large deformation of the drillstring with large slenderness ratio and the dynamic contact between the drillstring and the wellbore. ⁶ Recently, the multibody dynamic method is developed to investigate the drilling system.^7,8 In the typical model, the slender flexible drillstring is simulated with absolute nodal coordinate formulation (ANCF) ⁹ beam element. Meanwhile, the contact between the drillstring and the wellbore is modeled based on Hertz contact theory. ⁶ Compared with FEM, the multibody dynamic method can well simulate not only the three-dimensional rotation, large deformation, and large displacement of the drillstring but also the dynamic contact of the system. Furthermore, this method also has significant advantage in the calculation speed. ¹⁰

Since the slender drillstring is rotating and sliding in the wellbore filled up with drilling fluid, the interactions of the drillstring with the drilling fluid may have some influence on the dynamic behavior of the drilling system. Although there has been considerable studies in modeling and analyzing the dynamic characteristics of the system, such as the coupled vibration of the drillstring system,^11,12 the study that takes into account the coupling of the drillstring and drilling fluid is still lacking,^1,13 especially for investigating the effects of the drilling fluid on the drilling process.

In this article, a new simulation model of the drilling system considering the coupling of the drillstring and the drilling fluid is established using the multibody dynamic approach. In this model, the motion of the drilling fluid is described using the Arbitrary Lagrangian–Eulerian (ALE) method, ¹⁴ the drillstring is modeled as ANCF beam with several contact detection points, and then the contact between the drillstring and the wellbore is simulated based on Hertz contact theory. Furthermore, the rock penetration process has also been considered. Therefore, this model can be used to simulate the coupled vibration and the drilling process of the system with emphasis on the effects of drilling fluid.

The following context is organized in this way. Section “The dynamic model of the drilling system” establishes the dynamic equations of the drilling system using the multibody system dynamic approach. Then, the effects of the drilling fluid on the coupled vibration and the drilling process of the system are investigated in section “The effects of drilling fluid on the dynamic behavior of the drilling system.” Conclusions are given in section “Conclusion,” subsequently.

The dynamic model of the drilling system

The drilling system is a complex fluid–solid coupling dynamic system, as shown in Figure 1. The system has two subsystems: the drillstring and the drilling fluid. The drillstring subsystem consists of several rigid and flexible bodies with constraints between each other. And, the drilling fluid subsystem is composed of internal, annulus, and under-bit fluids. To establish the dynamic model of the drilling system, the dynamic equations of the drilling fluid and drillstring are derived, respectively.

Figure 1.

The schematic representation of the drilling system.

The dynamic equations of the drilling fluid

The relative flow of the drilling fluid is modeled based on framework of the ALE description. ¹⁴ The force of the drillstring acting on the drilling fluid is introduced through the drilling fluid transport motion, and the reaction force acting on the drillstring is taken as an external load.

To establish the dynamic equations of the drilling fluid, following assumptions are adopted:

1. The drilling fluid is a one-dimensional compressible flow.

2. Volumetric strain is proportional to the pressure and the increment area of the channel cross section is proportional to pressure.

3. The formula of pressure loss for steady flow is also available for the unsteady flow.

4. The axial stretch effect of the channel is considered.

5. There is no heat exchange in the subsystem.

The generalized coordinates and mass equation of the drilling fluid

Fluid flux $\varphi$ is used to describe the amount of the fluid mass passing through one flow section from a certain time. Generally, the time when the pipe is filled with fluid whose pressure is zero is chosen as the reference time (or the initial time). The fluid density ${\rho }_0$ and flow section area $A_0$ at the reference time are adopted as constant parameters. Then, the equivalent displacement u is used to describe the motion of the fluid, defined as

$$u=\frac{\varphi }{{\rho }_0{A}_0}$$

(1)

Thus

$$\frac{\partial u}{\partial t}=\frac{1}{{\rho }_0{A}_0}\frac{\partial \varphi }{\partial t}=\frac{\rho A}{{\rho }_0{A}_0}\upsilon$$

(2)

where $\rho ,A,\mathrm{and}\upsilon$ are the real-time density, section area, and velocity of the drilling fluid, respectively.

As shown in Figure 2, assuming that the channel length changed from initial length $\mathrm{d}x$ to current length ds, and at this moment, the density and flow section area are $\rho$ and A, and the inflow mass at section 1 is $\varphi$ , the outflow mass at section 2 is $\varphi +(\partial \varphi /\partial x)\mathrm{d}x$ .

Figure 2.

The schematic representation of the fluid mass conservation.

According to the mass conservation

$${\rho }_0{A}_0\mathrm{d}x-\rho A\mathrm{d}s=\frac{\partial \varphi }{\partial x}\mathrm{d}x$$

(3)

And then

$$\frac{\partial u}{\partial x}=\frac{1}{{\rho }_0{A}_0}\frac{\partial \varphi }{\partial x}=1-\frac{\rho A}{{\rho }_0{A}_0}\frac{\mathrm{d}s}{\mathrm{d}x}$$

(4)

Equations (2) and (4) are the derivative expressions of u. Adopting the equivalent displacement to describe the motion of the fluid can automatically satisfy the mass equation, and therefore is convenient.

The momentum equations of the fluid and the reaction force acting on the string

Since the drilling fluid is composed of internal, annulus, and under-bit fluids, therefore, the momentum equations of these three kinds of fluids will be derived, respectively.

As shown in Figure 3, the mass center of the internal fluid always coincides with the section center of the drillstring

$${\mathbf{r}}_{\mathrm{in}}={\mathbf{r}}_p$$

(5)

where ${\mathbf{r}}_{\mathrm{in}}$ and ${\mathbf{r}}_p$ are the section centroid position vectors of the drilling fluid and drillstring, respectively.

Figure 3.

The schematic representation of the momentum balance of the internal fluid.

According to the derivation rule of material derivative, the velocity and acceleration of the internal fluid are expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{\frac{\mathrm{d}{\mathbf{r}}_{\mathrm{in}}}{\mathrm{d}t}}& =\frac{\partial }{\partial t}{\mathbf{r}}_p+v\frac{\partial }{\partial s}{\mathbf{r}}_p={\stackrel{·}{\mathbf{r}}}_p+v{\mathbf{t}}_p\hfill \\ \multicolumn{1}{c}{\frac{{\mathrm{d}}^2{\mathbf{r}}_{\mathrm{in}}}{\mathrm{d}{t}^2}}& =\frac{\partial }{\partial t}\frac{\mathrm{d}{\mathbf{r}}_{\mathrm{in}}}{\mathrm{d}t}+v\frac{\partial }{\partial s}\frac{\mathrm{d}{\mathbf{r}}_{\mathrm{in}}}{\mathrm{d}t}\hfill \\ \multicolumn{1}{c}{}& ={\stackrel{··}{\mathbf{r}}}_p+2v{\stackrel{·}{\mathbf{t}}}_p+(\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial s}+\frac{\stackrel{·}{\lambda }v}{1+\lambda }){\mathbf{t}}_p+v^2\frac{\partial {\mathbf{t}}_p}{\partial s}\hfill \end{array}$$

(6)

where ${\stackrel{·}{\mathbf{r}}}_p$ and ${\stackrel{··}{\mathbf{r}}}_p$ are the velocity and acceleration vectors of the drillstring section center, respectively; ${\mathbf{t}}_p$ is the axial direction vector of the drillstring, ${\stackrel{·}{\mathbf{t}}}_p$ ; $\partial {\mathbf{t}}_p/\partial s$ denote the time and location gradient of ${\mathbf{t}}_p$ , respectively; v is the velocity of the internal fluid and $\partial v/\partial t,\partial v/\partial s$ denote the time and location gradient of v; $\lambda =\mathrm{d}s/\mathrm{d}x-1$ denotes the axial strain of drillstring; and $\stackrel{·}{\lambda }$ is the axial strain rate.

There are several forces acting on the internal fluid: the gravity, the pressure force of the end face, the pressure, and friction force caused by drillstring inter wall. They are expressed as follows:

• Pressure force on section 1

  $$[\mathrm{pA}-\frac{1}{2}\frac{\partial (\mathrm{pA})}{\partial s}\mathrm{d}s]({\mathbf{t}}_p-\frac{1}{2}\frac{\partial {\mathbf{t}}_p}{\partial s}\mathrm{d}s)$$

• Pressure force on section 2

  $$[\mathrm{pA}+\frac{1}{2}\frac{\partial (\mathrm{pA})}{\partial s}\mathrm{d}s]({\mathbf{t}}_p+\frac{1}{2}\frac{\partial {\mathbf{t}}_p}{\partial s}\mathrm{d}s)$$

• Gravity

  $$\rho A\mathbf{g}\mathrm{d}s$$

• Pressure force on the drillstring inter wall

  $$-\mathrm{pA}\frac{\partial {\mathbf{t}}_p}{\partial s}\mathrm{d}s$$

• Friction force on the drillstring inter wall

  $$\pi D_i\tau {\mathbf{t}}_p\mathrm{d}s$$

where $D_i$ is the inner diameter of the channel, and $\tau$ is the shear stress of the fluid boundary. The other normal force acting on the fluid by drillstring is ${\mathbf{F}}_n\mathrm{d}s$ .

According to Newton’s second law, the acceleration of internal fluid is balanced with external force

$$\begin{array}{c}\rho A[{\stackrel{··}{\mathbf{r}}}_p+2v{\stackrel{·}{\mathbf{t}}}_p+(\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial s}+\frac{\stackrel{·}{\lambda }v}{1+\lambda }){\mathbf{t}}_p+v^2\frac{\partial {\mathbf{t}}_p}{\partial s}]\mathrm{d}s\hfill \\ =[\mathrm{pA}-\frac{1}{2}\frac{\partial (\mathrm{pA})}{\partial s}\mathrm{d}s]({\mathbf{t}}_p-\frac{1}{2}\frac{\partial {\mathbf{t}}_p}{\partial s}\mathrm{d}s)\hfill \\ -[\mathrm{pA}+\frac{1}{2}\frac{\partial (\mathrm{pA})}{\partial s}\mathrm{d}s]({\mathbf{t}}_p+\frac{1}{2}\frac{\partial {\mathbf{t}}_p}{\partial s}\mathrm{d}s)\hfill \\ +\rho A\mathbf{g}\mathrm{d}s+\mathrm{pA}\frac{\partial {\mathbf{t}}_p}{\partial s}\mathrm{d}s-\pi D_i\tau {\mathbf{t}}_p\mathrm{d}s+{\mathbf{F}}_n\mathrm{d}s\hfill \end{array}$$

(7)

Three-dimensional dynamic equation of the internal fluid is obtained by simplifying equation (7), expressed as

$$\begin{array}{c}\frac{\partial }{\partial t}\left(\rho \mathrm{Av}\frac{\mathrm{d}s}{\mathrm{d}x}{\mathbf{t}}_p\right)+\rho \mathrm{Av}\frac{\partial }{\partial t}\left(\frac{\mathrm{d}s}{\mathrm{d}x}{\mathbf{t}}_p\right)\hfill \\ +\frac{\partial }{\partial x}[(\rho A{v}^2+\mathrm{pA}){\mathbf{t}}_p]+\rho A\frac{\mathrm{d}s}{\mathrm{d}x}({\stackrel{··}{\mathbf{r}}}_p-\mathbf{g})\hfill \\ +p_{c}A\frac{\mathrm{d}s}{\mathrm{d}x}{\mathbf{t}}_p-\mathrm{pA}\frac{\partial {\mathbf{t}}_p}{\partial x}-\frac{\mathrm{d}s}{\mathrm{d}x}{\mathbf{F}}_n=\mathbf{0}\hfill \end{array}$$

(8)

where $p_c$ is the flow pressure loss, $p_{c}A=\pi D\tau$ .

Since

$$p\frac{\partial A}{\partial x}<<A\frac{\partial p}{\partial x}$$

The additional pressure caused by the area change in channel can be ignored.

The one-dimensional dynamic equations of the internal fluid can be obtained by projecting equation (8) to the direct ${\mathbf{t}}_p$

$$\begin{array}{l}\frac{\partial }{\partial t}\left(\rho Av\frac{\mathrm{d}s}{\mathrm{d}x}\right)+\rho Av\frac{\partial }{\partial t}\left(\frac{\mathrm{d}s}{\mathrm{d}x}\right)+\frac{\partial (\rho A{v}^2+pA)}{\partial x}\\ +A\frac{\mathrm{d}s}{\mathrm{d}x}\left[\rho t_p^{\mathrm{T}}({\stackrel{\mathrm{··}}{r}}_p-g)+p_c\right]=0\end{array}$$

(9)

The first term represents the time rate of change of fluid momentum, the second term represents the additional inertia force caused by channel expansion, the third term represents the gradient of fluid hydrodynamic pressure, and the fourth term represents the external forces including relative motion, gravity, and friction. When the effect of channel expansion is not considered, equation (9) can be degenerated to classical one-dimensional flow equation

$$\frac{\partial }{\partial t}(\rho \mathrm{Av})+\frac{\partial (\rho A{v}^2+\mathrm{pA})}{\partial x}+\rho A{\mathbf{t}}_p^{\mathrm{T}}({\stackrel{··}{\mathbf{r}}}_p-\mathbf{g})+p_{c}A=0$$

Considering that the pressure and friction forces by the drillstring inner wall is the reciprocal force between the drillstring and the drilling fluid, rejecting the two forces in equation (8), the remaining forces are the reaction forces of the internal fluid acting on the drillstring, expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathbf{F}}_p}& =\frac{\mathrm{d}s}{\mathrm{d}x}(\mathrm{pA}\frac{\partial {\mathbf{t}}_p}{\partial s}+{\mathbf{F}}_n-p_{c}A{\mathbf{t}}_p)\hfill \\ \multicolumn{1}{c}{}& =\frac{\partial }{\partial t}\left(\rho \mathrm{Av}\frac{\mathrm{d}s}{\mathrm{d}x}{\mathbf{t}}_p\right)+\rho \mathrm{Av}\frac{\partial }{\partial t}\left(\frac{\mathrm{d}s}{\mathrm{d}x}{\mathbf{t}}_p\right)\hfill \\ \multicolumn{1}{c}{}& +\frac{\partial }{\partial x}[(\rho A{v}^2+\mathrm{pA}){\mathbf{t}}_p]+\rho A\frac{\mathrm{d}s}{\mathrm{d}x}({\stackrel{··}{\mathbf{r}}}_p-\mathbf{g})\hfill \\ \multicolumn{1}{c}{}& =\frac{\partial }{\partial t}(\rho \mathrm{Av})\mathbf{r}{\prime }_p+2\rho \mathrm{Av}\stackrel{·}{\mathbf{r}}{\prime }_p+\rho A\frac{\mathrm{d}s}{\mathrm{d}x}({\stackrel{··}{\mathbf{r}}}_p-\mathbf{g})\hfill \\ \multicolumn{1}{c}{}& +\frac{\partial }{\partial x}[(\rho A{v}^2+\mathrm{pA}){\mathbf{t}}_p]\hfill \end{array}$$

(10)

As shown in Figure 4, although the mass center of the annulus fluid not always coincides with the section centroid of drillstring, the projections of the two vectors along the borehole direction are equal. In addition, the static moment of wellbore cross section relative to coordinate system o-xyz equals to the sum of that of annulus fluid section and drillstring section. Then, the following relationship can be obtained

$$\begin{array}{c}{\mathbf{t}}_w^{\mathrm{T}}{\mathbf{r}}_{\mathrm{ex}}={\mathbf{t}}_w^{\mathrm{T}}{\mathbf{r}}_p={\mathbf{t}}_w^{\mathrm{T}}{\mathbf{r}}_w\hfill \\ {\mathbf{r}}_{\mathrm{ex}}{A}_o+{\mathbf{r}}_p(A_w-A_o)={\mathbf{r}}_w{A}_w\hfill \end{array}$$

(11)

where ${\mathbf{r}}_{\mathrm{ex}}$ , ${\mathbf{r}}_p$ , and ${\mathbf{r}}_w$ are the section centroid location vectors of the annulus, drillstring, and borehole, respectively; ${\mathbf{t}}_w$ is the direction vector of the borehole axis; $A_w$ is the cross-sectional area of the borehole; and $A_{\mathrm{o}}$ is the channel area of the annulus fluid. Therefore, the section centroid location of the annulus fluid is

$${\mathbf{r}}_{\mathrm{ex}}=-(\frac{A_w}{A_o}-1){\mathbf{r}}_p+\frac{A_w}{A_o}{\mathbf{r}}_w$$

(12)

where ${\mathbf{r}}_w$ is decided by projecting the drillstring centroid to the borehole axis. So, there is relationship as

$$\Delta {\mathbf{r}}_w={\mathbf{t}}_w{\mathbf{t}}_w^{\mathrm{T}}\Delta {\mathbf{r}}_p$$

(13)

Figure 4.

The schematic representation of the momentum balance of the annulus fluid.

As well as

$$\begin{array}{c}\frac{\partial }{\partial t}{\mathbf{r}}_w={\mathbf{t}}_w{\mathbf{t}}_w^{\mathrm{T}}\frac{\partial }{\partial t}{\mathbf{r}}_p\hfill \\ \frac{\partial }{\partial s}{\mathbf{r}}_w={\mathbf{t}}_w{\mathbf{t}}_w^{\mathrm{T}}\frac{\partial }{\partial s}{\mathbf{r}}_p\hfill \end{array}$$

(14)

According to the derivation rule of the material derivative, the velocity and acceleration of the annulus fluid are obtained as

$$\begin{array}{c}\frac{\mathrm{d}{\mathbf{r}}_{\mathrm{ex}}}{\mathrm{d}t}=(\frac{\partial }{\partial t}+v_o\frac{\partial }{\partial s}){\mathbf{r}}_{\mathrm{ex}}={\mathbf{E}}_w(\frac{\partial }{\partial t}+v_o\frac{\partial }{\partial s}){\mathbf{r}}_p\hfill \\ \frac{{\mathrm{d}}^2{\mathbf{r}}_{\mathrm{ex}}}{\mathrm{d}{t}^2}={(\frac{\partial }{\partial t}+v_o\frac{\partial }{\partial s})}^2{\mathbf{r}}_{\mathrm{ex}}=(\frac{\partial }{\partial t}+v_o\frac{\partial }{\partial s})\hfill \\ \left[{\mathbf{E}}_w(\frac{\partial }{\partial t}+v_o\frac{\partial }{\partial s}){\mathbf{r}}_p\right]\hfill \end{array}$$

(15)

where

$${\mathbf{E}}_w=\frac{-(A_w-A_o)\mathbf{I}+A_w{\mathbf{t}}_w{\mathbf{t}}_w^{\mathrm{T}}}{A_o}$$

is the flow transformation matrix of the annulus fluid, ${\mathbf{E}}_w$ is an invertible matrix, and

$${\mathbf{E}}_w^{-1}=\frac{-A_o\mathbf{I}+A_w{\mathbf{t}}_w{\mathbf{t}}_w^{\mathrm{T}}}{A_w-A_o}$$

Comparing equations (6) and (15), we find the transformation matrix ${\mathbf{E}}_w$ is the only difference in the velocity and acceleration between internal and annulus fluid; therefore, when ${\mathbf{t}}_w$ coincides with ${\mathbf{t}}_p$ , the expressions of the velocity and acceleration of these two fluids are exactly the same.

The forces acting on the annulus fluid are gravity, pressure force on the end face, pressure, and friction force caused by the drillstring outer wall and borehole inter wall. They are expressed as follows:

• Pressure force on section 1

  $$(p-\frac{1}{2}\frac{\partial p}{\partial s}\mathrm{d}s)A_w{\mathbf{t}}_w-[(p{A}_p)-\frac{1}{2}\frac{\partial (p{A}_p)}{\partial s}\mathrm{d}s]({\mathbf{t}}_p-\frac{1}{2}\frac{\partial {\mathbf{t}}_p}{\partial s}\mathrm{d}s)$$

• Pressure force on section 2

  $$(p+\frac{1}{2}\frac{\partial p}{\partial s}\mathrm{d}s)A_w{\mathbf{t}}_w-[(p{A}_p)+\frac{1}{2}\frac{\partial (p{A}_p)}{\partial s}\mathrm{d}s]({\mathbf{t}}_p+\frac{1}{2}\frac{\partial {\mathbf{t}}_p}{\partial s}\mathrm{d}s)$$

• Gravity

  $$\rho A_o\mathbf{g}\mathrm{d}s$$

• Pressure force on the outer wall of drillstring

  $$p(A_w-A_o)\frac{\partial {\mathbf{t}}_p}{\partial s}\mathrm{d}s$$

• Pressure force on the inter wall of borehole

  $$-p{A}_w\frac{\partial {\mathbf{t}}_w}{\partial s}\mathrm{d}s$$

• Friction force on the outer wall of drillstring

  $$\pi D_o{\tau }_p{\mathbf{t}}_p\mathrm{d}s$$

  where $D_o$ is the external diameter of the drillstring, and ${\tau }_p$ is the fluid shear stress on the string outer wall.

• Friction force on the inter wall of borehole

  $$\pi D_w{\tau }_w{\mathbf{t}}_w\mathrm{d}s$$

where $D_w$ is the diameter of the wellbore, and ${\tau }_w$ is the fluid shear stress on the inter wall of borehole.

The normal force acting on the fluid by drillstring is ${\mathbf{F}}_p$ . The normal force acting on the fluid by the borehole is ${\mathbf{F}}_w$ .

According to Newton’s second law, the acceleration of annulus fluid is balanced with external force

$$\begin{array}{c}\rho A_o{\mathbf{E}}_w[{\stackrel{··}{\mathbf{r}}}_p+2v{\stackrel{·}{\mathbf{t}}}_p+(\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial s}+\frac{\stackrel{·}{\lambda }v}{1+\lambda }){\mathbf{t}}_p+v^2\frac{\partial {\mathbf{t}}_p}{\partial s}]\mathrm{d}s\hfill \\ =(p-\frac{1}{2}\frac{\partial p}{\partial s}\mathrm{d}s)A_w{\mathbf{t}}_w-[(p{A}_p)-\frac{1}{2}\frac{\partial (p{A}_p)}{\partial s}\mathrm{d}s]({\mathbf{t}}_p-\frac{1}{2}\frac{\partial {\mathbf{t}}_p}{\partial s}\mathrm{d}s)\hfill \\ -(p+\frac{1}{2}\frac{\partial p}{\partial s}\mathrm{d}s)A_w{\mathbf{t}}_w+[(p{A}_p)+\frac{1}{2}\frac{\partial (p{A}_p)}{\partial s}\mathrm{d}s]({\mathbf{t}}_p+\frac{1}{2}\frac{\partial {\mathbf{t}}_p}{\partial s}\mathrm{d}s)\hfill \\ +\rho A\mathbf{g}\mathrm{d}s-p(A_w-A_o)\frac{\partial {\mathbf{t}}_p}{\partial s}\mathrm{d}s+p{A}_w\frac{\partial {\mathbf{t}}_w}{\partial s}\mathrm{d}s-\pi D_o{\tau }_p{\mathbf{t}}_p\mathrm{d}s\hfill \\ -\pi D_w{\tau }_w{\mathbf{t}}_w\mathrm{d}s+{\mathbf{F}}_p\mathrm{d}s+{\mathbf{F}}_w\mathrm{d}s\hfill \end{array}$$

(16)

Simplified

$$\begin{array}{c}\rho A_o{\mathbf{E}}_w[{\stackrel{··}{\mathbf{r}}}_p+2v{\stackrel{·}{\mathbf{t}}}_p+(\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial s}+\frac{\stackrel{·}{\lambda }v}{1+\lambda }){\mathbf{t}}_p+v^2\frac{\partial {\mathbf{t}}_p}{\partial s}]\frac{\mathrm{d}s}{\mathrm{d}x}\hfill \\ =\frac{\partial (p{A}_p)}{\partial s}\frac{\mathrm{d}s}{\mathrm{d}x}{\mathbf{t}}_p-\frac{\partial (p{A}_w)}{\partial s}\frac{\mathrm{d}s}{\mathrm{d}x}{\mathbf{t}}_w+\rho A_o\mathbf{g}\frac{\mathrm{d}s}{\mathrm{d}x}-p_{c}A\frac{\mathrm{d}s}{\mathrm{d}x}{\mathbf{t}}_p\hfill \\ +{\mathbf{F}}_p\frac{\mathrm{d}s}{\mathrm{d}x}+{\mathbf{F}}_w\frac{\mathrm{d}s}{\mathrm{d}x}\hfill \end{array}$$

(17)

Combining the frictions caused by drillstring and the borehole, the three-dimensional dynamic equations of the annulus fluid are expressed as

$$\begin{array}{c}\frac{\partial }{\partial t}\left(\rho A_{o}v\frac{\mathrm{d}s}{\mathrm{d}x}{\mathbf{t}}_p\right)+\rho \mathrm{Av}\frac{\partial }{\partial t}\left(\frac{\mathrm{d}s}{\mathrm{d}x}{\mathbf{t}}_p\right)+\frac{\partial }{\partial x}(\rho A_o{v}^2{\mathbf{t}}_p)-\frac{\partial (p{A}_o)}{\partial x}{\mathbf{t}}_p\hfill \\ +p_c{A}_o\frac{\mathrm{d}s}{\mathrm{d}x}+\rho A_o({\stackrel{··}{\mathbf{r}}}_p-{\mathbf{E}}_w^{-1}\mathbf{g})\frac{\mathrm{d}s}{\mathrm{d}x}-{\mathbf{E}}_w^{-1}({\mathbf{F}}_p+{\mathbf{F}}_w)\frac{\mathrm{d}s}{\mathrm{d}x}=0\hfill \end{array}$$

(18)

where $p_c$ is the flow pressure loss.

Because

$$p\frac{\partial A}{\partial x}<<A_o\frac{\partial p}{\partial x}$$

the additional pressure caused by the change of channel area can be ignored.

The one-dimensional dynamic equations of the annulus fluid are obtained by projecting equation (18) to the direct ${\mathbf{t}}_p$ and assuming ${\mathbf{t}}_p^{\mathrm{T}}{\mathbf{E}}_w^{-1}({\mathbf{F}}_p+{\mathbf{F}}_w)=0$ and then expressed as

$$\begin{array}{l}\frac{\partial }{\partial t}\left(\rho A_{o}v\frac{\mathrm{d}s}{\mathrm{d}x}\right)+\rho A_{o}v\frac{\partial }{\partial t}\left(\frac{\mathrm{d}s}{\mathrm{d}x}\right)+\frac{\partial (\rho A_o{v}^2+p{A}_o)}{\partial x}\\ +A_o\frac{\mathrm{d}s}{\mathrm{d}x}[\rho t_p^{\mathrm{T}}({\stackrel{\mathrm{··}}{r}}_p-E_w^{-1}g)+p_c]=0\end{array}$$

(19)

In this equation, the first term represents the time rate of change of fluid momentum, the second term represents the additional force caused by channel expansion, the third term represents the gradient of fluid hydrodynamic pressure, and the fourth term represents the external forces including relative motion, gravity, and friction. When the effect of channel expansion is not considered, equation (19) can also be degenerated to classical one-dimensional flow equation

$$\frac{\partial }{\partial t}(\rho \mathrm{Av})+\frac{\partial (\rho A{v}^2+\mathrm{pA})}{\partial x}+A[\rho {\mathbf{t}}_p^{\mathrm{T}}({\stackrel{··}{\mathbf{r}}}_p-{\mathbf{E}}_w^{-1}\mathbf{g})+p_c]=0$$

(20)

Investigating the one-dimensional flow equations of the internal fluid and the annulus fluid, we find the main difference is that the internal fluid is projecting to ${\mathbf{t}}_p$ , while the annulus is ${\mathbf{E}}_w^{-1}{\mathbf{t}}_p$ . However, the other terms are the same. Therefore, when the axes of the borehole and string are completely coincident, that is, ${\mathbf{t}}_p={\mathbf{t}}_w$ , and then, ${\mathbf{E}}_w^{-1}{\mathbf{t}}_p={\mathbf{t}}_p$ , the one-dimensional flow equations of these two type fluids are the same.

Altogether, the pressure and friction force on the drillstring outer is the reciprocal forces between the drillstring and the fluid, the same as the reciprocal forces between the borehole and fluid. By rejecting the interaction forces mentioned above in equation (19), the remaining forces are the reaction forces by the annulus fluid acting on the drillstring, expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{\mathbf{F}}& =(p{A}_p\frac{\partial {\mathbf{t}}_p}{\partial s}-p_c{A}_p{\mathbf{t}}_p)\frac{\mathrm{d}s}{\mathrm{d}x}+{\mathbf{F}}_p\frac{\mathrm{d}s}{\mathrm{d}x}\hfill \\ \multicolumn{1}{c}{}& ={\mathbf{E}}_w[\frac{\partial }{\partial t}(\rho \mathrm{Av})\mathbf{r}{\prime }_p+2\rho \mathrm{Av}\stackrel{·}{\mathbf{r}}{\prime }_p+\rho A({\stackrel{··}{\mathbf{r}}}_p-{\mathbf{E}}_w^{-1}\mathbf{g})\frac{\mathrm{d}s}{\mathrm{d}x}]\hfill \\ \multicolumn{1}{c}{}& +{\mathbf{E}}_w\frac{\partial }{\partial x}[(\rho A{v}^2+\mathrm{pA}){\mathbf{t}}_p]+(p{A}_w\frac{\partial {\mathbf{t}}_p}{\partial s}-p_c{A}_p{\mathbf{t}}_p)\frac{\mathrm{d}s}{\mathrm{d}x}\hfill \end{array}$$

(21)

where

$${\mathbf{E}}_w=\frac{-(A_w-A_o)\mathbf{I}+A_w{\mathbf{t}}_w{\mathbf{t}}_w^{\mathrm{T}}}{A_o}\mathrm{and}{\mathbf{E}}_w^{-1}=\frac{-A_o\mathbf{I}+A_w{\mathbf{t}}_w{\mathbf{t}}_w^{\mathrm{T}}}{A_w-A_o}$$

The momentum equations of the under-bit fluid can be obtained by replacing ${\mathbf{t}}_p$ with ${\mathbf{t}}_w$ in equation (8) and ignoring the effects of channel expansion; therefore, the three-dimensional dynamic equations of the under-bit fluid can be expressed as

$$\frac{\partial }{\partial t}(\rho \mathrm{Av}){\mathbf{t}}_w+\frac{\partial }{\partial x}[(\rho A{v}^2+\mathrm{Ap})]{\mathbf{t}}_w+A(-\rho \mathbf{g}+p_c{\mathbf{t}}_w)=0$$

(22)

The one-dimensional dynamic equations are obtained by projecting equation (22) to ${\mathbf{t}}_w$

$$\frac{\partial }{\partial t}(\rho \mathrm{Av})+\frac{\partial }{\partial x}(\rho A{v}^2+\mathrm{Ap})+A(-\rho {\mathbf{t}}_w^{\mathrm{T}}\mathbf{g}+p_c)=0$$

(23)

Obviously, there is no reaction force acting on the drillstring caused by under-bit fluid.

The constitutive equation of the drilling fluid

The constitutive equation of one-dimensional compressible fluid is

$$\frac{1}{E}=\frac{1}{\rho }\frac{\mathrm{d}\rho }{\mathrm{d}p}+\frac{1}{A}\frac{\mathrm{d}A}{\mathrm{d}p}$$

(24)

where E is the equivalent elastic modulus of drilling fluid, $(1/\rho )(\mathrm{d}\rho /\mathrm{d}p)$ is the fluid compression coefficient represented by $\alpha$ , while $(1/A)(\mathrm{d}A/\mathrm{d}p)$ is the channel elasticity coefficient represented by $\beta$ .

Integrating equation (24) once and according to equation (4), the following equation can be obtained

$$\begin{array}{cc}\multicolumn{1}{c}{p}& =E\ln \left(\frac{\rho A}{{\rho }_0{A}_0}\right)\hfill \\ \multicolumn{1}{c}{}& =E[\ln (1-\frac{\partial u}{\partial x})-\ln \left(\frac{\mathrm{d}s}{\mathrm{d}x}\right)]\hfill \\ \multicolumn{1}{c}{}& \approx -E(u\prime +{\epsilon }_p)\hfill \end{array}$$

(25)

Furthermore, the expansion of the channel area can be expressed as

$$A=A_0{e}^{\beta p}$$

(26)

Considering that

$$p\frac{\partial A}{\partial x}=(\beta p)A\frac{\partial p}{\partial x}<\frac{p}{E}A\frac{\partial p}{\partial x}<<A\frac{\partial p}{\partial x}$$

the additional pressure caused by the elastic transmutation of the channel area is so small that it is neglected. Therefore, a simplified constitutive equation is used

$$\{\begin{array}{c}\epsilon =u\prime +{\epsilon }_p\hfill \\ \multicolumn{1}{c}{p=-E\epsilon }\end{array}$$

(27)

where the equivalent elastic modulus E is coalescence of the fluid compression coefficient $\alpha$ and the channel elasticity coefficient $\beta$ .

The fluid compression coefficient is related to the fluid type, temperature, and pressure. However, the change of $\alpha$ is small when pressure is in the range of 0–49 MPa; therefore, the compression coefficient can be replaced by that of the water in the state of 50°C and 50 MPa.

The channel elasticity coefficient can be easily derived based on the elastic theory, the elasticity coefficient of internal flow channel is

$$\beta =\frac{2}{E_s}(\frac{D_o^2+D_i^2}{D_o^2-D_i^2}+{\mu }_s)$$

The elasticity coefficient of the annulus flow channel is

$$\beta =\frac{2}{(D_w^2-D_o^2)}[\frac{D_w^2}{E_f}(1+{\mu }_f)+\frac{D_o^2}{E_s}(\frac{D_o^2+D_i^2}{D_o^2-D_i^2}-{\mu }_s)]$$

where $E_s,{\mu }_s$ , and $E_f,{\mu }_f$ are the elasticity modulus and Poisson’s ratio of the steel and the formation, respectively.

Flow pressure loss of the fluid

There is flow resistance term in the dynamic equation of the drilling fluid, and the resistance, presented as flow pressure loss, is mainly caused by on-way resistance and related with the viscosity property, the state, and the rheological model of the drilling fluid. In drilling engineering, the drilling fluid viscosity is measured by direct-indicating viscometer. The rheological properties of the drilling fluid are complex. The most widely used models are power law, Bingham, Casson, and Herschel–Bulkley (HB). In this article, the power law model assumed that the fluid shear stress $\tau$ and velocity gradient follow an exponential law, which is adopted as follows

$$\tau =K{\left(\frac{\mathrm{d}v}{\mathrm{d}x}\right)}^n$$

(28)

where K is the consistency index, and n is the flow behavior index. They can be obtained as ¹⁵

$$\begin{array}{c}n=\frac{1}{\lg 2}\lg \left(\frac{{\varphi }_{600}}{{\varphi }_{300}}\right)\hfill \\ K=0.511\frac{{\varphi }_{300}}{{511}^n}\mathrm{Pa}·s^n\hfill \end{array}$$

(29)

where ${\varphi }_{300}$ and ${\varphi }_{600}$ are torsion frames of the rotational viscometer which rotate at the speeds 300 and 600 r/min, respectively.

Then, the pressure loss of the internal fluid can be expressed as

$$p_c=\{\begin{array}{cc}\frac{32}{\mathrm{Re}}\frac{\rho v^2}{D_i}\hfill & \mathrm{Re}<3470-1370n\hfill \\ \multicolumn{1}{c}{\frac{(\lg n+3.93)}{25{\mathrm{Re}}^{(0.25-\frac{\lg n}{7})}}\frac{\rho v^2}{D_i}}& \mathrm{Re}<3470-1370n\hfill \end{array}$$

(30)

where

$$\begin{array}{c}v=v_i\hfill \\ \mathrm{Re}=\frac{8\rho v^2}{K}{\left[\frac{D_i}{(6+\frac{2}{n})v}\right]}^n\hfill \end{array}$$

And, the pressure loss of the annulus fluid can be obtained as

$$p_c=\{\begin{array}{cc}\frac{48}{\mathrm{Re}}\frac{\rho v^2}{(D_w-D_o)}\hfill & \mathrm{Re}<3470-1370n\hfill \\ \multicolumn{1}{c}{\frac{(\lg n+3.93)}{25{\mathrm{Re}}^{(0.25-\frac{\lg n}{7})}}\frac{\rho v^2}{(D_w-D_o)}}& \mathrm{Re}<3470-1370n\hfill \end{array}$$

(31)

where

$$\begin{array}{c}v=v_e-(1-K_c)v_p\hfill \\ \mathrm{Re}=\frac{12\rho v^2}{K}{\left[\frac{D_w-D_o}{(8+\frac{4}{n})v}\right]}^n\hfill \end{array}$$

In the equations above, $v_p$ is the velocity along the axial direction of the drillstring, $v_i$ and $v_o$ are the velocity of the internal and the annulus fluids relative to the drillstring, and $K_C$ is the adhesion factor that represents the equivalent flow coefficient caused by the annulus fluid adhesion to the drillstring and can be obtained by the following equations

$$K_c=\{\begin{array}{c}1-\frac{1}{\ln {\sigma }^2}-\frac{1}{1-{\sigma }^2}\mathrm{Laminar\; flow}\hfill \\ \multicolumn{1}{c}{1-\frac{1-\sqrt{{\sigma }^3-{\sigma }^2+\sigma }}{1-{\sigma }^2}\mathrm{Turbulent\; flow}}\end{array}$$

(32)

where $\sigma =D_o/D_w$ .

The dynamic equations of the drilling fluid

By coupling the equation of momentum and constitutive equation, the dynamic equations of the drilling fluid are obtained as

$$\{\begin{array}{c}{\rho }_i^{(0)}{A}_i^{(0)}\{{\stackrel{··}{u}}_i(1+\lambda )+2{\stackrel{·}{u}}_i\stackrel{·}{\lambda }+(1-{u\prime }_i)[{\mathbf{t}}_p^{\mathrm{T}}({\stackrel{··}{\mathbf{r}}}_p-\mathbf{g})+f_i]\}+\frac{\partial }{\partial x}({\rho }_i{A}_i{v}_i^2+p_i{A}_i)=0\hfill \\ \multicolumn{1}{c}{{\rho }_e^{(0)}{A}_e^{(0)}\{{\stackrel{··}{u}}_e(1+\lambda )+2{\stackrel{·}{u}}_e\stackrel{·}{\lambda }+(1-{u\prime }_e)[{\mathbf{t}}_p^{\mathrm{T}}({\stackrel{··}{\mathbf{r}}}_p-{\mathbf{E}}_w^{-1}\mathbf{g})+f_e]\}+\frac{\partial }{\partial x}({\rho }_e{A}_e{v}_e^2+p_e{A}_e)=0}\\ \multicolumn{1}{c}{{\rho }_b^{(0)}{A}_b^{(0)}\{{\stackrel{··}{u}}_b+(1-{u\prime }_b)[{\mathbf{t}}_p^{\mathrm{T}}({\stackrel{··}{\mathbf{r}}}_p-\mathbf{g})+f_b]\}+\frac{\partial }{\partial x}({\rho }_b{A}_b{v}_b^2+p_b{A}_b)=0}\end{array}$$

(33)

and the reaction force acting on the drillstring is

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathbf{F}}_p}& ={\rho }_i^{(0)}{A}_i^{(0)}[{\stackrel{··}{u}}_i(1+\lambda ){\mathbf{t}}_p+2(1+\lambda ){\stackrel{·}{u}}_i{\stackrel{·}{\mathbf{t}}}_p+2{\stackrel{·}{u}}_i\stackrel{·}{\lambda }{\mathbf{t}}_p]\hfill \\ \multicolumn{1}{c}{}& +{\rho }_e^{(0)}{A}_e^{(0)}{\mathbf{E}}_w[{\stackrel{··}{u}}_e(1+\lambda ){\mathbf{t}}_p+2(1+\lambda ){\stackrel{·}{u}}_e{\stackrel{·}{\mathbf{t}}}_p+2{\stackrel{·}{u}}_e\stackrel{·}{\lambda }{\mathbf{t}}_p]\hfill \\ \multicolumn{1}{c}{}& +\frac{\partial }{\partial x}\left[({\rho }_i{A}_i{\stackrel{·}{u}}_i^2\frac{1+\lambda }{1-{u\prime }_i}-E_i{A}_i\frac{{u\prime }_i+\lambda }{1+\lambda }){\mathbf{t}}_p\right]\hfill \\ \multicolumn{1}{c}{}& +{\mathbf{E}}_w\frac{\partial }{\partial x}\left[({\rho }_e{A}_e{\stackrel{·}{u}}_e^2\frac{1+\lambda }{1-{u\prime }_e}-E_e{A}_e\frac{{u\prime }_e+\lambda }{1+\lambda }){\mathbf{t}}_p\right]\hfill \\ \multicolumn{1}{c}{}& +{\rho }_i^{(0)}{A}_i^{(0)}(1-u{\prime }_i)({\stackrel{··}{\mathbf{r}}}_p-\mathbf{g})+{\rho }_e^{(0)}{A}_e^{(0)}(1-{u\prime }_o)({\mathbf{E}}_w{\stackrel{··}{\mathbf{r}}}_p-\mathbf{g})\hfill \\ \multicolumn{1}{c}{}& -p_e{A}_w\frac{\partial {\mathbf{t}}_w}{\partial x}+2\pi R_w{\tau }_w\frac{\partial s}{\partial x}{\mathbf{t}}_w\hfill \end{array}$$

(34)

In the above equations, the subscripts i, e, and b represent internal, annulus, and under-bit fluids, respectively.

Dividing the channel into several segments and creating fluid node at the section, the neighboring two nodes construct one element. The equivalent displacement u and fluid strain $\epsilon$ at the two nodes are selected as the generalized coordinate $\mathbf{q}$ , which is

$$\mathbf{q}=[\begin{array}{cccc}{u}_1\hfill & {\epsilon }_1\hfill & u_2\hfill & {\epsilon }_2\hfill \end{array}{]}^{\mathrm{T}}$$

(35)

The equivalent displacement gradient of the element can be obtained at the nodes

$$\begin{array}{c}u{\prime }_1={\epsilon }_1-{\lambda }_1\hfill \\ u{\prime }_2={\epsilon }_2-{\lambda }_2\hfill \end{array}$$

(36)

Hermite interpolation is applied to describe the equivalent displacement of the element

$$u=N_1{u}_1+N_{2}u{\prime }_1+N_3{u}_2+N_{4}u{\prime }_2$$

(37)

where

$$\begin{array}{c}{N}_1=\frac{1}{4}(\xi -1{)}^2(2+\xi )N_2=\frac{L}{8}(\xi -1{)}^2(\xi +1)\hfill \\ N_3=\frac{1}{4}(\xi +1{)}^2(2-\xi )N_4=\frac{L}{8}(\xi +1{)}^2(\xi -1)\hfill \\ \xi =\frac{2x}{L}-1\hfill \end{array}$$

By substituting equation (36) into equation (37), the interpolation formula using generalized coordinates can be obtained as

$$u=\mathbf{Nq}-N_2{\lambda }_1-N_4{\lambda }_2$$

(38)

where

$$\mathbf{N}=[\begin{array}{cccc}{N}_1\hfill & N_2\hfill & N_3\hfill & N_4\hfill \end{array}]$$

is the interpolation function matrix of the element. Furthermore, the partial derivatives of the generalized gradient are obtained as

$$\begin{array}{c}u\prime =\mathbf{N}\prime \mathbf{q}-N{\prime }_2{\lambda }_1-N{\prime }_4{\lambda }_2\hfill \\ \stackrel{·}{u}\prime =\mathbf{N}\prime \stackrel{·}{\mathbf{q}}-N{\prime }_2{\stackrel{·}{\lambda }}_1-N{\prime }_4{\stackrel{·}{\lambda }}_2\hfill \\ \stackrel{·}{u}=\mathbf{N}\stackrel{·}{\mathbf{q}}-N_2{\stackrel{·}{\lambda }}_1-N_4{\stackrel{·}{\lambda }}_2\hfill \\ \stackrel{··}{u}=\mathbf{N}\stackrel{··}{\mathbf{q}}-N_2{\stackrel{··}{\lambda }}_1-N_4{\stackrel{··}{\lambda }}_2\hfill \end{array}$$

(39)

where

$$\begin{array}{c}\mathbf{N}\prime =[\begin{array}{cccc}N{\prime }_1\hfill & N{\prime }_2\hfill & N{\prime }_3\hfill & N{\prime }_4\hfill \end{array}]\hfill \\ N{\prime }_1=\frac{3}{2L}({\xi }^2-1)N{\prime }_2=\frac{1}{4}(\xi -1)(3\xi +1)\hfill \\ N{\prime }_3=\frac{3}{2L}(1-{\xi }^2)N{\prime }_4=\frac{1}{4}(\xi +1)(3\xi -1)\hfill \end{array}$$

is the gradient matrix of element interpolation function.

Using Galerkin integration method and adopting the shape function as the weight function, the element integral equations of the internal fluid can be obtained as

$$\begin{array}{l}{\displaystyle \underset{0}{\overset{L}{\int }}{N}^{\mathrm{T}}{\rho }_i^{(0)}{A}_i^{(0)}[{\stackrel{··}{u}}_i(1+\lambda )+2{\stackrel{.}{u}}_i\stackrel{.}{\lambda }+(1-{u^{\prime }}_i)t_p^{\mathrm{T}}({\stackrel{\mathrm{··}}{r}}_p-g)+(1-{u^{\prime }}_i)f_i]\mathrm{d}xx}\\ +{\displaystyle \underset{0}{\overset{L}{\int }}{N}^{\mathrm{T}}\frac{\partial }{\partial x}({\rho }_i{A}_i{v}_i^2+p_i{A}_i)\mathrm{d}x}=0\end{array}$$

(40)

The weak solution of equation (40) can be obtained by integrating the second term using partial integration once

$$\begin{array}{l}{\displaystyle \underset{0}{\overset{L}{\int }}{N}^{\mathrm{T}}{\rho }_i^{(0)}{A}_i^{(0)}[(1+\lambda ){\stackrel{··}{u}}_i+2{\stackrel{.}{u}}_i\stackrel{.}{\lambda }+(1-{u^{\prime }}_i)t_p^{\mathrm{T}}({\stackrel{\mathrm{··}}{r}}_p-g)+(1-{u^{\prime }}_i)f]\mathrm{d}x}\\ -{\displaystyle \underset{0}{\overset{L}{\int }}{N^{\prime }}^{\mathrm{T}}({\rho }_i{A}_i{v}_i^2+p_i{A}_i)\mathrm{d}x}+{N^{\mathrm{T}}[{\rho }_i{A}_i{v}_i^2+p_i{A}_i]|}_0^L=0\end{array}$$

(41)

In the equations above, u is the only requested one-order continuous that satisfies the property of Hermite interpolation

$${{\mathbf{N}}^{\mathrm{T}}[p_i{A}_i]|}_0^L$$

(42)

Equation (42) represents the integral boundary of the pressure and has three cases:

1. Middle node. Integral boundary term is counteracted by the neighborhood elements.

2. Displacement boundary node. Integral boundary term is replaced by the constrained force.

3. Force boundary node. Boundary term is applied using force.

So, when the dynamic equations of all the fluid elements are assembled, the boundary term (42) would not be retained.

Finally, the finite element equation of the drilling internal fluid can be expressed as

$$\begin{array}{c}{\mathbf{M}}_f^{(i)}{\stackrel{··}{\mathbf{q}}}_i+{\mathbf{K}}_f^{(i)}{\mathbf{q}}_i-{\mathbf{W}}_f^{(i)}+{\mathbf{Q}}_{\lambda }^{(i)}+{\mathbf{Q}}_p^{(i)}+{\mathbf{Q}}_c^{(i)}-{\mathbf{Q}}_v^{(i)}+{\mathbf{Q}}_{\mathrm{end}}^{(i)}=0\hfill \\ {\mathbf{M}}_f^{(i)}=\underset{0}{\overset{L}{\int }}{\rho }_i^{(0)}{A}_i^{(0)}{\mathbf{N}}^{\mathrm{T}}\mathbf{N}(1+\lambda )\mathrm{d}x\hfill \\ {\mathbf{K}}_f^{(i)}=\underset{0}{\overset{L}{\int }}\frac{\mathrm{EA}{\mathbf{N}\prime }^{\mathrm{T}}{\mathbf{N}\prime }^{\mathrm{T}}}{(1+\lambda )\mathrm{d}x}\hfill \\ {\mathbf{W}}_f^{(i)}=\underset{0}{\overset{L}{\int }}\rho A{\mathbf{N}}^{\mathrm{T}}({\mathbf{g}}^{\mathrm{T}}{\mathbf{t}}_b)(1-u{\prime }_i)\mathrm{d}x\hfill \\ {\mathbf{Q}}_{\lambda }^{(i)}=\underset{0}{\overset{L}{\int }}2{\mathbf{N}}^{\mathrm{T}}\stackrel{·}{u}\stackrel{·}{\lambda }\mathrm{d}x\hfill \\ {\mathbf{Q}}_p^{(i)}=\underset{0}{\overset{L}{\int }}\rho A{\mathbf{N}}^{\mathrm{T}}({{\stackrel{··}{\mathbf{r}}}_b}^{\mathrm{T}}{\mathbf{t}}_b)(1-u{\prime }_i)\mathrm{d}x\hfill \\ {\mathbf{Q}}_c^{(i)}=\underset{0}{\overset{L}{\int }}{\mathbf{N}}^{\mathrm{T}}(1+\lambda )p_c\mathrm{d}x\hfill \\ {\mathbf{Q}}_v^{(i)}=\underset{0}{\overset{L}{\int }}{\mathbf{N}\prime }^{\mathrm{T}}{\rho }_i^{(0)}{A}_i^{(0)}{\stackrel{·}{u}}_i^2\frac{1+\lambda }{1-u{\prime }_i}\mathrm{d}x\hfill \\ {\mathbf{Q}}_{\mathrm{end}}^{(i)}={\left[\begin{array}{cccc}{({\rho }_i{A}_i{v}_i^2)|}_{x=0}\hfill & 0\hfill & -{({\rho }_i{A}_i{v}_i^2)|}_{x=L}\hfill & 0\hfill \end{array}\right]}^{\mathrm{T}}\hfill \end{array}$$

(43)

where ${\mathbf{M}}_f^{(i)}$ is the element mass matrix, ${\mathbf{K}}_f^{(i)}$ is the element stiffness matrix, ${\mathbf{Q}}_{\lambda }^{(i)}$ is the generalized additional force caused by the expansion of the channel, ${\mathbf{W}}_f^{(i)}$ is the generalized gravity, and ${\mathbf{Q}}_p^{(i)}$ is the generalized force caused by the transport motion of the string. ${\mathbf{Q}}_c^{(i)}$ is the generalized force of the viscous force. ${\mathbf{Q}}_v^{(i)}$ is the generalized force caused by the fluid convection, and ${\mathbf{Q}}_{\mathrm{end}}^{(i)}$ represents the integral boundary.

Similarly, the finite element equations of the annulus and under-bit fluids can be obtained using weighted residual method

$$\{\begin{array}{c}{\mathbf{M}}_f^{(i)}{\stackrel{··}{\mathbf{q}}}^{(i)}+{\mathbf{K}}_f^{(i)}{\mathbf{q}}^{(i)}-{\mathbf{W}}_f^{(i)}+{\mathbf{Q}}_{\chi }^{(i)}+{\mathbf{Q}}_b^{(i)}+{\mathbf{Q}}_c^{(i)}-{\mathbf{Q}}_v^{(i)}+{\mathbf{Q}}_{\mathrm{end}}^{(i)}=0\hfill \\ \multicolumn{1}{c}{{\mathbf{M}}_f^{(e)}{\stackrel{··}{\mathbf{q}}}^{(e)}+{\mathbf{K}}_f^{(e)}{\mathbf{q}}^{(e)}-{\mathbf{W}}_f^{(e)}+{\mathbf{Q}}_{\chi }^{(e)}+{\mathbf{Q}}_b^{(e)}+{\mathbf{Q}}_c^{(e)}-{\mathbf{Q}}_v^{(e)}+{\mathbf{Q}}_{\mathrm{end}}^{(e)}=0}\\ \multicolumn{1}{c}{{\mathbf{M}}_f^{(b)}{\stackrel{··}{\mathbf{q}}}^{(b)}+{\mathbf{K}}_f^{(b)}{\mathbf{q}}^{(b)}-{\mathbf{W}}_f^{(b)}+{\mathbf{Q}}_c^{(b)}-{\mathbf{Q}}_v^{(b)}+{\mathbf{Q}}_{\mathrm{end}}^{(b)}=0}\end{array}$$

(44)

And the reaction forces of the drilling fluid acting on drillstring can be expressed as

$$\begin{array}{c}{\mathbf{F}}_p=\underset{0}{\overset{L}{\int }}{\rho }_i^{(0)}{A}_i^{(0)}{\left(\frac{\partial {\mathbf{r}}_p}{\partial {\mathbf{q}}_p}\right)}^{\mathrm{T}}\hfill \\ [{\stackrel{··}{u}}_i(1+\lambda ){\mathbf{t}}_p+2{\stackrel{·}{u}}_i(1+\lambda ){\stackrel{·}{\mathbf{t}}}_p+2{\stackrel{·}{u}}_i\stackrel{·}{\lambda }{\mathbf{t}}_p+(1-u{\prime }_i)({\stackrel{··}{\mathbf{r}}}_p-\mathbf{g})]\mathrm{d}x\hfill \\ +\underset{0}{\overset{L}{\int }}{\rho }_e^{(0)}{A}_e^{(0)}{(\frac{1}{\alpha }-1)}^2{\left(\frac{\partial {\mathbf{r}}_p}{\partial {\mathbf{q}}_p}\right)}^{\mathrm{T}}\hfill \\ [{\stackrel{··}{u}}_e(1+\lambda ){\mathbf{t}}_p+2{\stackrel{·}{u}}_e(1+\lambda ){\stackrel{·}{\mathbf{t}}}_p+2{\stackrel{·}{u}}_e\stackrel{·}{\lambda }{\mathbf{t}}_p+(1-u{\prime }_e)({\stackrel{··}{\mathbf{r}}}_p-{\mathbf{E}}_w^{-1}\mathbf{g})]\mathrm{d}x\hfill \\ +\underset{0}{\overset{L}{\int }}{\rho }_e^{(0)}{A}_e^{(0)}\frac{(2\alpha -1)}{{\alpha }^2}{\left({\mathbf{t}}_w^{\mathrm{T}}\frac{\partial {\mathbf{r}}_p}{\partial {\mathbf{q}}_p}\right)}^{\mathrm{T}}{\mathbf{t}}_w^{\mathrm{T}}\hfill \\ [{\stackrel{··}{u}}_e(1+\lambda ){\mathbf{t}}_p+2{\stackrel{·}{u}}_e(1+\lambda ){\stackrel{·}{\mathbf{t}}}_p+2{\stackrel{·}{u}}_e\stackrel{·}{\lambda }{\mathbf{t}}_p+(1-u{\prime }_e)({\stackrel{··}{\mathbf{r}}}_p-{\mathbf{E}}_w^{-1}\mathbf{g})]\mathrm{d}x\hfill \\ -\underset{0}{\overset{L}{\int }}{\left(\frac{\partial \lambda }{\partial {\mathbf{q}}_p}\right)}^{\mathrm{T}}({\rho }_i{A}_i{v}_i^2+p_i{A}_i)\mathrm{d}x\hfill \\ -\underset{0}{\overset{L}{\int }}{S}_1{\left(\frac{\partial \lambda }{\partial {\mathbf{q}}_p}\right)}^{\mathrm{T}}({\rho }_e{A}_e{v}_e^2+p_e{A}_e)\mathrm{d}x\hfill \\ -\underset{0}{\overset{L}{\int }}{S}_2{\left(\frac{\partial {\mathbf{t}}_p}{\partial {\mathbf{q}}_p}\right)}^{\mathrm{T}}{\mathbf{t}}_w({\rho }_e{A}_e{v}_e^2+p_e{A}_e)\mathrm{d}x\hfill \\ +\underset{0}{\overset{L}{\int }}{\left(\frac{\partial {\mathbf{r}}_p}{\partial {\mathbf{q}}_p}\right)}^{\mathrm{T}}[\frac{A_w(A_w-A_e)}{A_e}{p}_e\frac{\partial {\mathbf{t}}_w}{\partial x}+{\rho }_e^{(0)}{A}_e^{(0)}(1-u{\prime }_e)f_w{\mathbf{t}}_w]\mathrm{d}x\hfill \end{array}$$

(45)

where

$$\begin{array}{c}\alpha =\frac{A_o}{A_w}\hfill \\ S_1=\frac{1}{{\alpha }^2}[(1-\alpha {)}^2+(2\alpha -1)({\mathbf{t}}_p^{\mathrm{T}}{\mathbf{t}}_w{)}^2]\hfill \\ S_2=\frac{1}{{\alpha }^2}(1+\lambda )(2\alpha -1)({\mathbf{t}}_w^{\mathrm{T}}{\mathbf{t}}_p)\hfill \end{array}$$

(46)

The dynamic equations of the drillstring system

The drillstring is modeled as a complex rigid–flexible coupling system which consists of rigid bodies, flexible beam elements, constraints, and dynamic loads. The dynamic equations of rigid body and flexible beam are divided, respectively.

Dynamic equations of rigid bodies

The top drive, the stabilizer, and the bit are modeled as rigid bodies in this article. The generalized coordinate $\mathbf{q}$ of a rigid body can be defined as the position coordinate ${\mathbf{r}}_b$ of its mass center in the global coordinate system, and unit quaternion ${\mathbf{\theta }}_b$ that is defined based on Euler’s rotation theorem

$${\mathbf{q}}_b={\left[\begin{array}{cc}{\mathbf{r}}_b^{\mathrm{T}}\hfill & {\mathbf{\theta }}_b^{\mathrm{T}}\hfill \end{array}\right]}^{\mathrm{T}}$$

(47)

Based on the Lagrange equation, the dynamic equation of the rigid body can be expressed as ¹⁶

$$\{\begin{array}{c}{\mathbf{M}}_b{\stackrel{··}{\mathbf{r}}}_b-{\mathbf{F}}_b=\mathbf{0}\hfill \\ \multicolumn{1}{c}{4{\mathbf{G}}_b^{\mathrm{T}}{\mathbf{J}}_b{\mathbf{G}}_b{\stackrel{··}{\mathbf{\theta }}}_b-8{\stackrel{·}{\mathbf{G}}}_b^{\mathrm{T}}{\mathbf{J}}_b{\mathbf{G}}_b{\stackrel{·}{\mathbf{\theta }}}_b+2{\mathbf{\theta }}_b{\sigma }_b-{\mathbf{Q}}_b=\mathbf{0}}\\ \multicolumn{1}{c}{{\mathbf{\theta }}_b^{\mathrm{T}}{\mathbf{\theta }}_b-1=0}\end{array}$$

(48)

where ${\mathbf{M}}_b$ and ${\mathbf{J}}_b$ are the mass matrix and inertia matrix of the body, ${\mathbf{F}}_b$ and ${\mathbf{Q}}_b$ are the generalized external force and moment, ${\sigma }_b$ is the Lagrange multiplier, and ${\mathbf{G}}_b$ is

$${\mathbf{G}}_b=\left[\begin{array}{cccc}-{\theta }_{b,1}\hfill & {\theta }_{b,0}\hfill & {\theta }_{b,3}\hfill & -{\theta }_{b,2}\hfill \\ \multicolumn{1}{c}{-{\theta }_{b,2}}& -{\theta }_{b,3}\hfill & {\theta }_{b,0}\hfill & {\theta }_{b,1}\hfill \\ \multicolumn{1}{c}{-{\theta }_{b,3}}& {\theta }_{b,2}\hfill & -{\theta }_{b,1}\hfill & {\theta }_{b,0}\hfill \end{array}\right]$$

Dynamic equations of ANCF beam

The ANCF beam elements which can describe the motion, rotation, and deformation of the drillstring are used to model the flexible components of the drillstring. Although the theory of the ANCF beam has been proposed by Zhao and Ren, ¹⁷ it is preferred to introduce this element briefly for the integrity of this article.

A two-node beam element with cross-sectional area A and initial length L is shown in Figure 5, and the generalized coordinates of the nodes I and II are

$${\mathbf{q}}^{\mathrm{I}}=\left[\begin{array}{c}{\mathbf{r}}^{\mathrm{I}}\hfill \\ \multicolumn{1}{c}{{\mathbf{\theta }}^{\mathrm{I}}}\\ \multicolumn{1}{c}{{\epsilon }^{\mathrm{I}}}\end{array}\right],{\mathbf{q}}^{\mathrm{II}}=\left[\begin{array}{c}{\mathbf{r}}^{\mathrm{II}}\hfill \\ \multicolumn{1}{c}{{\mathbf{\theta }}^{\mathrm{II}}}\\ \multicolumn{1}{c}{{\epsilon }^{\mathrm{II}}}\end{array}\right]$$

(49)

where $\mathbf{r}$ is the global displacement vector of the node, $\epsilon$ is the normal strain at axial line of the beam, and $\mathbf{\theta }$ is the unit quaternion, which satisfied the following constraint equation

$$\mathbf{C}\left({\mathbf{q}}_e\right)=\left[\begin{array}{c}{\mathbf{\theta }}^{\mathrm{I}}·{\mathbf{\theta }}^{\mathrm{I}}-1\hfill \\ \multicolumn{1}{c}{{\mathbf{\theta }}^{\mathrm{II}}·{\mathbf{\theta }}^{\mathrm{II}}-1}\end{array}\right]=0$$

(50)

Figure 5.

The schematic representation of ANCF beam element.

Then, the generalized elemental coordinate is

$${\mathbf{q}}_e=\left[\begin{array}{c}{\mathbf{q}}^{\mathrm{I}}\hfill \\ \multicolumn{1}{c}{{\mathbf{q}}^{\mathrm{II}}}\end{array}\right]$$

(51)

Assuming that the arc-length coordinate of an arbitrary point on the axis of the beam element is x, then the isoparametric coordinate $\xi$ is

$$\xi =\frac{x}{L}$$

(52)

According to the definition of normal strain, the derivatives of the centerline at nodes are

$${\left(\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}\xi }\right)}^{\mathrm{I}}=L(1+{\epsilon }^{\mathrm{II}}){\mathbf{t}}^{\mathrm{I}},{\left(\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}\xi }\right)}^{\mathrm{II}}=L(1+{\epsilon }^{\mathrm{II}}){\mathbf{t}}^{\mathrm{II}}$$

(53)

where $\epsilon$ is the axial strain. Then, the displacement vector $\mathbf{r}$ of the point is expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{\mathbf{r}\left(\xi \right)}& =N_1\left(\xi \right){\mathbf{r}}^{\mathrm{I}}+N_2\left(\xi \right){\left(\frac{d\mathbf{r}}{d\xi }\right)}^{\mathrm{I}}+N_3\left(\xi \right){\mathbf{r}}^{\mathrm{II}}+N_4\left(\xi \right){\left(\frac{d\mathbf{r}}{d\xi }\right)}^{\mathrm{II}}\hfill \\ \multicolumn{1}{c}{}& =N_1\left(\xi \right){\mathbf{r}}^{\mathrm{I}}+N_2\left(\xi \right)L(1+{\epsilon }^{\mathrm{I}}){\mathbf{t}}^{\mathrm{I}}\hfill \\ \multicolumn{1}{c}{}& +N_3\left(\xi \right){\mathbf{r}}^{\mathrm{II}}+N_4\left(\xi \right)L(1+{\epsilon }^{\mathrm{II}}){\mathbf{t}}^{\mathrm{II}}\hfill \end{array}$$

(54)

where $N_i$ is the shape function and can be given by Hermite interpolation

$$\begin{array}{c}{N}_1(\xi )=1-3{\xi }^2+2{\xi }^3\hfill \\ N_2(\xi )=\xi (1-\xi {)}^2\hfill \\ N_3(\xi )={\xi }^2(3-2\xi )\hfill \\ N_4(\xi )={\xi }^2(\xi -1)\hfill \end{array}$$

(55)

To describe the rotation of the beam element, the material reference frame $[\mathbf{m},\mathbf{n},\mathbf{t}]$ is set on the node of the element, and $\mathbf{t}$ is the normal direction of the section. Due to the assumption of the Euler–Bernoulli beam theory, the tangent vector $\mathbf{t}(\xi )$ can be obtained as

$$\mathbf{t}\left(\xi \right)=\frac{\mathbf{r}\prime }{|\mathbf{r}\prime |}$$

(56)

The symbol “′” denotes the derivative with respect to $\xi$ . As shown in Figure 5, an intermediate reference frame $[\overline{\mathbf{m}},\overline{\mathbf{n}},\overline{\mathbf{t}}]$ with $\overline{\mathbf{t}}=\mathbf{t}$ at $\xi$ can be obtained by rotating $[{\mathbf{m}}^{\mathrm{I}},{\mathbf{n}}^{\mathrm{I}},{\mathbf{t}}^{\mathrm{I}}]$ an angle $\mathrm{arc}\cos ({\mathbf{t}}^{\mathrm{I}}·\mathbf{t})$ with respect to the vector ${\mathbf{t}}^{\mathrm{I}}\times \mathbf{t}$ , the corresponding quaternion is expressed as ¹⁸

$$\overline{\mathbf{\theta }}\left(\xi \right)={\mathbf{\theta }}^{\mathrm{I}}\mathbf{\pi }$$

(57)

where

$$\pi =\left({\pi }_0,{\pi }_1,{\pi }_2,{\pi }_3\right)=\left(\sqrt{\frac{1+t·t^{\mathrm{I}}}{2}},0,\frac{-t·n^{\mathrm{I}}}{\sqrt{2\left(1+t·t^{\mathrm{I}}\right)}},\frac{t·m^{\mathrm{I}}}{\sqrt{2\left(1+t·t_1^{\mathrm{I}}\right)}}\right)$$

Obviously, the final $\mathbf{\theta }\mathbf{(}\mathbf{\xi }\mathbf{)}$ can be obtained by rotating $\overline{\mathbf{\theta }}(\xi )$ an angle $\beta (\xi )$ with respect to t

$$\mathbf{\theta }\mathbf{\left(}\mathbf{\xi }\mathbf{\right)}\mathbf{=}\overline{\mathbf{\theta }}\mathbf{\left(}\mathbf{\xi }\mathbf{\right)}\mathbf{(}\mathbf{cos}\frac{\mathbf{\beta }\mathbf{\left(}\mathbf{\xi }\mathbf{\right)}}{2}\mathbf{,}\mathbf{sin}\frac{\mathbf{\beta }\mathbf{\left(}\mathbf{\xi }\mathbf{\right)}}{2}\mathbf{,}0\mathbf{,}0\mathbf{)}$$

(58)

Then, the kinetic energy of the element is

$$T=\frac{1}{2}\rho L\underset{0}{\overset{1}{\int }}[A(\stackrel{·}{\mathbf{r}}·\stackrel{·}{\mathbf{r}})+4{\stackrel{·}{\theta }}^{\mathrm{T}}{\overline{\mathbf{E}}}^{\mathrm{T}}\mathbf{J}\overline{\mathbf{E}}\stackrel{·}{\theta }]d\xi =\frac{1}{2}{\stackrel{·}{\mathbf{q}}}_e^{\mathrm{T}}\mathbf{M}{\stackrel{·}{\mathbf{q}}}_e$$

(59)

where $\overline{\mathbf{E}}$ is

$$\overline{\mathbf{E}}=\left[\begin{array}{cccc}-{\theta }_1\hfill & {\theta }_0\hfill & {\theta }_3\hfill & -{\theta }_2\hfill \\ \multicolumn{1}{c}{-{\theta }_2}& -{\theta }_3\hfill & {\theta }_0\hfill & {\theta }_1\hfill \\ \multicolumn{1}{c}{-{\theta }_3}& {\theta }_2\hfill & -{\theta }_1\hfill & {\theta }_0\hfill \end{array}\right]$$

$\mathbf{M}$ is the mass matrix and can be expressed as

$$\mathbf{M}=\rho L\underset{0}{\overset{1}{\int }}[A{\left(\frac{\partial \mathbf{r}}{\partial {\mathbf{q}}_e}\right)}^{\mathrm{T}}\left(\frac{\partial \mathbf{r}}{\partial {\mathbf{q}}_e}\right)+4{\left(\frac{\partial \mathbf{\theta }}{\partial {\mathbf{q}}_e}\right)}^{\mathrm{T}}{\overline{\mathbf{E}}}^{\mathrm{T}}\mathbf{J}\overline{\mathbf{E}}\left(\frac{\partial \mathbf{\theta }}{\partial {\mathbf{q}}_e}\right)]\mathrm{d}\xi$$

The elastic potential work is

$$U=\frac{1}{2}L\underset{0}{\overset{1}{\int }}(\mathrm{EA}{\epsilon }^2+G{J}_T{\kappa }_1^2+G{J}_{22}{\kappa }_2^2+G{J}_T{\kappa }_3^2)d\xi$$

(60)

where E and G are the elastic and shear moduli, respectively, and ${\kappa }_1$ , ${\kappa }_2$ , and ${\kappa }_3$ are the torsion and two bending curvatures.

Therefore, the governing equations of the beam can be obtained by the Lagrange equation

$$\{\begin{array}{cc}\multicolumn{1}{c}{\mathbf{M}{\stackrel{··}{\mathbf{q}}}_e-\mathbf{P}-\mathbf{Q}+{\left(\frac{\partial \mathbf{C}}{\partial {\mathbf{q}}_e}\right)}^{\mathrm{T}}\mathbf{\lambda }\mathbf{=}}& 0\hfill \\ \multicolumn{1}{c}{\mathbf{C}\mathbf{\left(}{\mathbf{q}}_{\mathbf{e}}\mathbf{\right)}\mathbf{=}}& 0\hfill \end{array}$$

(61)

where $\mathbf{\lambda }$ is the corresponding Lagrange multiplier, and $\mathbf{Q}$ is the external generalized force. $\mathbf{P}=-\mathbf{V}{\stackrel{·}{\mathbf{q}}}_e-\partial U/\partial {\mathbf{q}}_e$ and $\mathbf{V}=\partial (\mathbf{M}{\stackrel{·}{\mathbf{q}}}_e)/{\mathbf{q}}_e-(1/2)(\partial (\mathbf{M}{\stackrel{·}{\mathbf{q}}}_e)/{\mathbf{q}}_e{)}^{\mathrm{T}}$ .

Boundary conditions

Rock penetration equation

A three-dimensional rate of penetration (ROP) equation that considering the cutting bit anisotropy and the formation anisotropy is adopted to describe the rock penetration process.

By defining the coordinate system $o-\xi \eta \zeta$ of the bit and the coordinate system $o-x\prime y\prime z\prime$ of the layer, $\xi$ is the axial direction of the bit, and $z\prime$ is the normal direction of the layer. In the rock penetration process, the relationship of the force $\mathbf{F}=[F_{\xi }{F}_{\eta }{F}_{\zeta }{]}^{\mathrm{T}}$ of bit on rock and the penetration velocity $\mathbf{v}=[v_{x\prime }{v}_{y\prime }{v}_{z\prime }]$ is ⁷

$$\left[\begin{array}{c}{v}_{x\prime }\hfill \\ \multicolumn{1}{c}{v_{y\prime }}\\ \multicolumn{1}{c}{v_{z\prime }}\end{array}\right]=\left[\begin{array}{ccc}1-h_1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1-h_2\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]·\mathbf{A}·\left[\begin{array}{c}{F}_{\xi }\hfill \\ \multicolumn{1}{c}{c_1{F}_{\eta }}\\ \multicolumn{1}{c}{c_2{F}_{\zeta }}\end{array}\right]·n_{\xi }$$

(62)

where $h_1$ , $h_2$ represent the formation anisotropic parameters, and $c_1$ , $c_2$ represent the bit anisotropic parameters. $\mathbf{A}$ is the transformation matrix from bit body coordinate to layer coordinate, and $n_{\xi }$ is the axial cutting parameter, which is related to rotary speed and the weight on the bit.

In the drilling simulation, the moving-down velocity and the rotary speed of the top drive are the main input parameters, while the bit weight and the actual drilling speed are the output parameters induced by the dynamic model. In order to determine the bit weight and the drilling speed, the forces acting on the bit by the wellbore should be calculated first. Therefore, the Hertz contact theory is adopted to model the interaction between the bit and the wellbore like an ellipsoid sphere in a bucket. Then, the bit weight could be easily counted. Furthermore, the drilling speed could be accomplished with equation (62), in which the other parameters are determined by the characterization of the stratum. It should be emphasized that the growth of the wellbore is implemented at the end of each integral step using the drilling speed multiplied with the step size which is along the direction of the bit axis.

Contact

The wellbore is modeled as a rigid hole with circular cross section, and a set of wellbore node is placed on the curved well centerline to describe the well trajectory, each wellbore node has nine parameters: the material coordinate, the position vector, the tangent direction, the radius, and the sliding friction coefficient. In petroleum engineering, the well trajectory is described with the measure depth (the arc length of the well trajectory from the start point of the well to the wellbore node), the inclination (the angle between the local tangent direction of the wellbore node with the vertical direction), and the azimuth (the angle between the projection of the trajectory on the horizontal plane with the north direction). And the parameters of the wellbore node can be transformed from the measure depth, the inclination, and the azimuth using the cubic Hermite interpolation. ⁶

Several detecting points are set on the axial line of the beam element. Therefore, the detection between the drillstring and the wellbore is transformed to the detection between detecting points and rigid cylinder, as shown in Figure 6.

Figure 6.

The contact model between drillstring and wellbore.

According to the Hertz contact theory, the contact force ${\mathbf{f}}_c$ is composed of support force $f_n$ and tangential friction force $f_t$ , which can be expressed as ¹⁹

$${\mathbf{f}}_c=f_n\mathbf{n}+f_t\mathbf{t}$$

(63)

where $\mathbf{n}$ is the unit normal direction of the contact, and $\mathbf{\tau }$ is the corresponding unit tangential vector.

The support force $f_n$ is a function of the penetration depth $\delta$ and the penetration velocity $\stackrel{·}{\delta }$ , expressed as

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

(64)

where k is the contact stiffness, e is the contact recovery index, and c is the contact damping.

The magnitude of friction force can be obtained by Coulomb model

$$f_{\tau }=\mu f_n$$

(65)

where $\mu$ is the friction coefficient. And, the friction direction is parallel to the velocity direction of the contact point.

Finally, the contact force is exerted on the nodes of the beam element. According to the virtual work principle, the generalized contact force acting on the element can be expressed as

$${\mathbf{Q}}_p={\mathbf{N}}^{\mathrm{T}}({\xi }_p){\mathbf{f}}_c$$

(66)

where $\mathbf{N}({\xi }_p)$ is the shape function of beam element, ${\xi }_p$ is the element coordinate of the detection point p, and ${\mathbf{f}}_c$ is the force exerted on the detecting point.

The dynamic equations of the drilling system

The string is divided into several segments, each segment is one element, and the two ends of the segment are selected as the node, then the dynamic simulation model of the string can be established. The nodes of internal and annulus fluids are appended to the string nodes that have the same location. The dynamic equations of the whole system can be obtained by assembling equations (44), (48), and (61). Using the multibody dynamic system method, these three equations can be expressed as

$$\{\begin{array}{cc}\mathbf{M}({\stackrel{·}{\mathbf{q}}}_s,{\mathbf{q}}_s,t){\stackrel{··}{\mathbf{q}}}_s+\mathbf{Q}({\stackrel{·}{\mathbf{q}}}_s,{\mathbf{q}}_s,t)+{\mathbf{A}}_s^{\mathrm{T}}{\mathbf{\lambda }}_A+{\left(\frac{\partial {\mathbf{C}}_s}{\partial {\mathbf{q}}_s}\right)}^{\mathrm{T}}{\mathbf{\lambda }}_C\hfill & =\mathbf{0}(a)\hfill \\ \multicolumn{1}{c}{{\mathbf{C}}_s({\mathbf{q}}_s,t)}& =\mathbf{0}(b)\hfill \\ \multicolumn{1}{c}{{\mathbf{A}}_s({\mathbf{q}}_s,t){\stackrel{·}{\mathbf{q}}}_s+{\mathbf{B}}_s({\mathbf{q}}_s,t)}& =\mathbf{0}(c)\hfill \end{array}$$

(67)

where equation (a) is the dynamic equation of the system, equation (b) is the holonomic constraint equation of the system, and equation (c) is the non-holonomic constraint equation; ${\mathbf{q}}_s$ is the generalized coordinate vector; $\mathbf{M}$ is the system mass matrix; ${\mathbf{\lambda }}_A$ is the corresponding Lagrange multiplier of the non-holonomic constraint; and ${\mathbf{\lambda }}_C$ is the corresponding Lagrange multiplier of the holonomic constraint.

The equations of the system form a set of differential algebraic equations (DAEs), which can be solved by the implicit first-order backward differentiation formula (BDF) method. ⁶

Effects of drilling fluid on the dynamic behavior of the drilling system

Coupled vibration of the drillstring

Due to the contact of the drillstring with the wellbore, the contact of the bit and formation, as well as the unbalance mass and the initial curvature of the drillstring, the coupled vibration will occur when the drillstring rotates in the wellbore. Typical forms are axial, torsional, and lateral vibrations. A coupling model of drilling system with a constant top drive rotary speed is established to investigate the influence of the drilling fluid on the coupled vibration.

Without loss of generality, the wellbore is assumed to be vertical and its measure depth is 2000 m, the dimension of the drillstring is shown in Figure 7, where “OD” and “ID” denote the outside and inside diameter, and the physical parameters of the drillstring and drilling fluid are shown in Table 1. The rotary speed of the top drive is 80 r/min. Through the mesh convergence test, the drillstring is discretized into 225 elements. And, the element number of internal and annulus fluids is equal to drillstring.

Figure 7.

The schematic representation of the analytical system.

Table 1.

The physical parameters of the analytical system.

Subsystem	Parameters	Value	Unit
Drillstring	Elasticity modulus	210	GPa
	Density	7850	kg/m3
	Poisson’s ratio	0.28
Drilling fluid	Flow rate	30	L/s
	Density	1350	kg/m3
	${\varphi }_{600}$	100
	${\varphi }_{300}$	60
	Stiffness	1 × 108	N/m
Contact	Damping	1 × 103	N s/m
	Friction coefficient	0.0005

For the convenience of discussion, another simulation model without drilling fluid, named as “drillstring model,” is established for comparison.

The axial vibration

The axial vibration can be characterized by the fluctuation of the bit pressure. The bit pressures of the two models are shown in Figure 8. It can be found that the fluctuation of the bit pressure approximates periodic fluctuation, and the difference of the two models is not obvious.

Figure 8.

The bit pressure of the two models.

The torsional vibration

The torsional vibration is usually caused by the friction between the drillstring and the wellbore. When the driver rotates at constant speed, the bit speed is 0–2 times that of top driver. As shown in Figure 9, the bit speed curves of the two models are almost consistent.

Figure 9.

The speed of the bit.

The whirling motion

The typical transverse vibration is the whirling motion. The unbalance mass or other interference force may cause the string revolve around the borehole axis and rotate on its axis.

Figures 10(a) and 11(a) show the deformation of the drillstring at the end of the simulation, and Figures 10(b)–(d) and 11(b)–(d) are the center whirl orbits of drillstring at measure depths 600, 1200, and 1900 m.

Figure 10.

Whirling motion of the coupling model (a) Deformation of the drillsting (t = 50 s). (b–d) Whirl orbits of drillstring.

Figure 11.

Whirling motion of the drillstring model (a) Deformation of the drillsting (t = 50 s). (b–d) Whirl orbits of drillstring.

Due to the lower drillstring subjected to compression, the whirling motion occurs more easily compared with the upper drillstring, as shown in Figure 10. According to the dynamic response of the system, we found that the whirling motion of the coupling model is more likely to occur at the constant rotary speed and the drillstring deformation is larger. In contrast, the drillstring remains in the center of the borehole when the effect of drilling fluid is ignored. Of particular note is that the helical deformation appeared in the collar of the coupling model and the collar outward against the borehole well.

In short, the drilling fluid has little influence on the axial and torsional vibration of the drillstring but leads to the whirling motion is easier to occur, and this influence needs more research in the future.

Effects of drilling fluid on the drilling process

Several simulation models with different drilling fluid parameters are established to investigate the effects of the flow rate $\varphi$ and density $\rho$ of drilling fluid on the drilling process. To simulate the deflection process, a concentrated force, named as thrust force, is applied on the ribs of the rotary steering tool (RST), the magnitude is 1.96 × 10⁴ N, and the direction is along the direction of deflection, the rotary speed of the top drive is 50 r/min, and the introducing speed is 0.01 m/s. The drilling process begins when the measure depth is 2000 m and ends when the measure depth is 2030 m.

The well trajectory and the inclination angle of the system with different parameters are shown in Figures 12 and 13. In the following section, the effects of the flow rate and density of the drilling fluid will be discussed, respectively.

Figure 12.

The well trajectory of the system with different drilling fluid parameters.

Figure 13.

The inclination angle of the system with different drilling fluid parameters.

Comparing the well trajectory and inclination angle of the simulation model with different flow rates with that of the simulation model without drilling fluid, it is found that the well trajectory of the system with low flow rate is close to the system without drilling fluid. As the flow rate increases, the difference becomes greater.

On the other hand, the inclination angle increases with the increase in the density of the drilling fluid, especially when the flow velocity increases at the same time.

For a further discussion, the effect of the drilling fluid on the drilling process and the concept of average buildup rate are used. Assuming that the new well trajectory approximates to a circle arc, then the average buildup rate is defined as

$$\lambda =30\frac{{\alpha }_2-{\alpha }_1}{L_2-L_1}$$

(68)

where ${\alpha }_1$ and ${\alpha }_2$ represent the inclination of the start point and end point of the new well trajectory, $L_1$ and $L_2$ represent the measure depth of the corresponding points. The average buildup rate of all the simulation models is shown in Table 2.

Table 2.

The average buildup rate of the simulation models.

Flow rate (L/s)	Density (kg/m3)	$\lambda (\circ /30\mathrm{m})$
0	0	5.594
30	1350	5.907
30	2350	6.034
60	1350	6.101
60	2350	6.251

It is found that the average buildup rate of the system increases with the increase in flow rate and density of the drilling fluid. But viewed as a whole, the influence of flow rate is more significantly than the density.

The difference of average buildup between the model with $\varphi =30\mathrm{L}/\mathrm{s}$ , $\rho =1350\mathrm{kg}/{\mathrm{m}}^3$ and the system without drilling fluid is just 5.6%. However, when the flow rate increases to 60 L/s, the difference of the two models is about 11.7%, which means that the influence of the drilling fluid could not be neglected when the flow rate of the drilling fluid is larger.

Conclusion

A fluid–solid coupling model of the drilling system based on the multibody dynamic system method is proposed for simulating the dynamic behavior of the drilling system with emphasis on the effects of drilling fluid. In this model, the relative flow of the drilling fluid is modeled using ALE description, and the dynamic equation of the slender flexible drillstring is established by the ANCF. Meanwhile, the model takes into account the contact between the drillstring and the wellbore as well as the rock penetration process, so this model can be used to simulate the coupled vibration and drilling process.

The effects of the drilling fluid on the dynamic behavior of the system are investigated. The simulation results show that the drilling fluid has significant influence on the whirling motion of the drillstring: the whirling motion will appear easily when the drillstring rotates in the wellbore filled with drilling fluid. Furthermore, the contact between the drillstring and wellbore is frequent and continuous. The deflecting capacity of the system increases with the increase in flow rate and the density of the drilling fluid. For the drilling system with high flow rate, the influence of the drilling fluid on the drilling process should be considered.

The simulation examples show that the proposed model is capable of describing the dynamic behavior of the drilling system and provides a potential accurate and efficient simulation model in drilling engineering.