## Services on Demand

## Article

## Indicators

## Related links

## Share

## Journal of the Brazilian Society of Mechanical Sciences and Engineering

*Print version* ISSN 1678-5878

### J. Braz. Soc. Mech. Sci. & Eng. vol.32 no.3 Rio de Janeiro July/Sept. 2010

#### http://dx.doi.org/10.1590/S1678-58782010000300008

**TECHNICAL PAPERS**

**Stochastic dynamics of a drill-string with uncertain weight-on-hook **

**T. G. Ritto ^{I}; C. Soize^{II}; Rubens Sampaio^{III} **

^{I}PUC-Rio Department of Mechanical Engineering 22453-900 Rio de Janeiro, RJ, Brazil. Université Paris-Est Lab. de Modélisation et Simulation Multi-Echelle 77454 Marne-la-Vallée, France. thiagoritto@gmail.com

^{II}Université Paris-Est Lab. de Modélisation et Simulation Multi-Echelle 77454 Marne-la-Vallée, France. christian.soize@univ-paris-est.fr

^{III}PUC-Rio Department of Mechanical Engineering 22453-900 Rio de Janeiro, RJ, Brazil. rsampaio@puc-rio.br

**ABSTRACT**

A drill-string is a slender structure that turns and drills into the rock in search of oil. There are many sources of uncertainties in this complex dynamical system. However, this article is concerned only with uncertainties in the weight-on-hook, which is the supporting force exerted by the hook at the top. A probabilistic model is constructed for the random variable related to the weight-on-hook using the Maximum Entropy Principle, and the random response of the system is computed through Monte Carlo simulations. The idea is to understand how the performance of the system (which is measured by the rate of penetration) is affected by the uncertainties of the weight-on-hook. The continuous system analyzed is discretized by means of the Finite Element Method and a computer code is developed to do the simulations.

**Keywords:** random weight-on-hook, uncertainty modeling, stochastic dynamics, drill-string dynamics, reduced models

**Introduction **

In a drilling operation there are many sources of uncertainties, such as the material properties of the column and of the drilling fluid; the dimensions of the system, especially of the borehole; the fluid-structure interaction; and the bit-rock interaction. To have an improved computational model with better predictability capacity, uncertainties should be quantified and included in the computational model, which then becomes a stochastic computational model (if the probability theory is used to quantify the uncertainties).

This paper is concerned with uncertainties in the weight-on-hook because it is one of the three parameters that are continuously controlled in a drilling operation. (The other two parameters are the rotational speed at the top and the inlet fluid velocity). Figure 1 shows the general scheme of the system analyzed. The motor torque is taken into account as a constant rotational speed at the top (Ω_{x}) and the weight-on-hook (which is the supporting force, *w*) is constant. The bit reacts to the dynamics of the column with the torque *t _{bit}* and the force ƒ

_{bit}.

The probability theory, which is a powerful tool to model uncertainties, is used in the present analysis. In a general way, there are two types of uncertainties: (1) one type related to the randomness of some parameter (aleatory uncertainties) and (2) other type related to the lack of knowledge of a given phenomenon (epistemic uncertainties), which is related to modeling errors. One way to take into account epistemic uncertainties is to use the nonparametric probabilistic approach (Soize, 2000) because it is able to model both parameter uncertainties as well as model uncertainties.

The analysis of the dynamics of flexible structures considering the deterministic and the stochastic problem (Sampaio and Ritto, 2008) is having an increasing attention of the scientific community over the years. Ritto et al. (2008), for example, analyze parameter and model uncertainties in the boundary condition of a beam.

In Tucker and Wang (1999), Khulief and AL-Naser (2005), Sampaio et al. (2007), Piovan and Sampaio (2009), for instance, the deterministic model of a drill-string system is analyzed. There are few articles treating the stochastic problem of the drill-string dynamics; see, for instance, Spanos et al. (1997, 2009), Ritto et al. (2009, 2010). In Kotsonis and Spanos (1997) a random weight-on-bit is analyzed in a simple two-degrees-of-freedom drillstring model, and in Spanos et al. (2009) lateral forces at the bit are modeled as stochastic. In Ritto et al. (2009) a probabilistic model (nonparametric probabilistic approach) is proposed for the bit-rock-interaction model, and in Ritto et al. (2010) a robust optimization problem is performed, where model uncertainties are modeled using the nonparametric probabilistic approach.

In the present paper, the focus of the stochastic analysis is only on the uncertain weight-on-hook, although there are other sources of uncertainties, such as the ones mentioned above. It is important to analyze a complex problem incrementally, in a way that the final result is not obscured by the influence of many factors at the same time. The parametric probabilistic approach is used to model the uncertainties in the weight-on-hook, and the probability density function of the random variable related to the weight-on-hook is constructed by the means of the Maximum Entropy Principle (Shannon, 1948; Jaynes, 1957a,b).

The strategy goes the following way. After the construction of the probability density function of the random weight-on-hook, a random generator is used to generate a value for the weight-on-hook for each Monte Carlo simulation (Rubinstein, 2007). The number of Monte Carlo simulations is chosen analyzing a convergence curve. The random response can finally be analyzed by the computations of some statistics, such as mean and confidence intervals.

In the first Section of the paper, the discretized system is presented and the element matrices are depicted, then, in the second Section, the reduced-order model is developed. The probabilistic model for the weight-on-hook is constructed in the third Section and the stochastic dynamical system is presented in the fourth Section. The numerical results are presented in the fifth Section and, finally, the concluding remarks are presented in the last Section.

**Nomenclature **

A | = | cross sectional area, m^{2} |

conv | = | convergence function, m^{2} |

D | = | diameter, m |

E | = | Young modulus, Pa |

g | = | acceleration of gravity, m/s^{2} |

G | = | shear modulus, Pa |

I | = | moment of inertia of the transversal section, m^{4} |

S | = | Shannon entropy measure |

t | = | time, s |

u | = | displacement in x-direction, m |

Z | = | regularizing function |

ƒ | = | force vector, N, N.m |

N | = | shape functions, m |

u | = | displacement vector, m, rad |

[C] | = | damping matrix, N.s/m, N.s.m |

[K] | = | stiffnes matriz, N/M, N.m |

[M] | = | mass matrix, kg, kg.m^{2} |

**Greek Symbols**

= | assumes value 1 if x belongs to B and 0 otherwise | |

δ | = | coeffictient of variation |

ρ | = | density, kg/m^{3} |

Ω_{x} | = | rotational speed at x = 0, rad/s |

σ | = | standard deviation, N |

θ_{x} | = | rotational about x-axis, rad/s |

[Φ] | = | modal basis, m, rad |

**Subscripts **

br | bit-rock |

f | fluid |

g | geometric (for [K]) and gravity (for f) |

inside | |

o | outside |

p | polar |

S | static response |

u | displacement in x-direction |

θ_{x} | rotational about x-axis |

**Final Discretized System **

The present work is concerned with the stochastic response of the system when the weight-on-hook (axial force) is random. To focus the attention of the analysis on the stochastic problem, and to speed up the numerical simulations, a simplification of the model found in Ritto et al. (2009) is used for the present analysis. The lateral displacements are assumed to be small and, therefore, they are neglected; and the system is linearized about the prestressed state **u**_{s}. This state is computed through the equation, **u**_{s} = [*K*]^{-1} (**f**_{g} + **f**_{c} + **f**_{ƒ} ) where [*K*] is the stiffness matrix of the system, **f**_{g} is the gravity force, **f**_{c} is the concentrated reaction force at the bit, and **f**_{f} is the fluid axial force. The final discretized system considering the prestressed state is written as:

in which the response **ū(= ū - ū*** _{s}*) is represented in a subspace

*V*⊂

_{m}*, where*

^{m}*m*equals the number of degrees of freedom of the system. [

*M*], [

*C*], and [

*K*] are the classical mass, damping and stiffness matrices, [

*K*(

_{g}**u**

_{S})] is the geometric stiffness matrix (due to the finite strain formulation),

**g**is the force due to the Dirichlet boundary condition (imposed rotational speed at the top) and

**f**

_{br}are the forces due to the bit-rock interactions.

The proportional damping matrix is constructed a posteriori [C] = [*C*] = *a*[*M*]+*b*([*K*]+[*K _{g}* (

**u**

_{S})]), where

*a*and

*b*are positive constants. The finite element approximation of the displacement fields is written as

where *u* is the axial displacement, *θ _{x }*is the rotation about the

*x*-axis, ξ =

*x/l*is the element coordinate,

_{e}*l*is the element length,

_{e}**N**are the shape functions

and the element node displacements are

where the exponent *T* means transposition.

The element forces and matrices are depicted in the sequence. The mass [*M*]^{(e)} and stiffness [*K*]^{(e)} element matrices are written as

and

where *ρ* is the density, *A* is the cross sectional area, *I _{p}* is the polar moment of inertia,

*E*is the elasticity modulus,

*G*is the shear modulus, and the space derivative (

*d/d*ξ) is denoted by ('). The element geometric stiffness matrix is written as

where and L_{p4} = ∫_{A}(*y*^{4} + z^{4}). The gravity element force and the fluid element force are written as

and

where *g* is the gravity acceleration, *M*_{ƒ} is the fluid mass per unit length, *A _{i} *and

*A*are the cross-sectional areas corresponding to the inner and outer diameters of the column,

_{0}*D*is the diameter,

*U*is the flow speed outside the column,

_{0}*C*

_{ƒ}is a fluid viscous damping coefficient and

*p*is the pressure inside the column (see Paidoussis et al. (2008) and Ritto et al. (2009) for details). The concentrated force at the bit is written as

_{i}where *ƒ*_{c} is the initial reaction force at the bit. The bit-rock interaction force is written as

in which *ƒ*_{bit} is the axial force and *t _{bit}* is the torque about the

*x*-axis. These two functions can be written as (Tucker and Wang, 2003)

where Z(_{bit}) is the regularizing function such that

In the above equation, *a*_{1}, ... , *a*_{5} are positive constants that depend on the bit and rock characteristics as well as on the weight-on-bit (*ƒ*_{bit}).

**Reduced-Order Model **

To accelerate the computations of the dynamical system response, a reduced-order model is constructed. One way to construct it is to project the nonlinear dynamical equation (Eq. (1)) on a subspace *V _{n}* ⊂

*, with*

^{m}*n*<<

*m*. In Trindade et al. (2005), the Karhunen-Loève decomposition is used to reduce a coupled axial/bending beam dynamics subjected to impacts. In the present paper, the basis used for the reduction corresponds to a basis of normal modes, which are obtained from the following generalized eigenvalue problem

where *Φ*_{i} is the *i*-th normal mode and *ω*_{i} is the *i*-th natural frequency. Using the representation

and substituting it in the equation of motion yields

where [*Φ*] is a (*m* × *n*) real matrix composed by *n* normal modes obtained using the prestressed configuration, Eq. (14). Projecting the equation on the subspace spanned by these normal modes yields

which can be rewritten as

in which

are the reduced matrices. The advantage of using a reduced-order model is to end up with a diagonal reduced matrix of size 22, instead of a banded finite element matrix of size 114 (values used in the numerical simulations).

**Probabilistic Model of the Weight-on-Hook **

The uncertainties in the weight-on-hook are modeled using the parametric probabilistic approach; therefore, the weight-on-hook *w* is modeled as a random variable *W*. To be coherent with the physics of the problem, the probability density function of *W* is constructed by means of the Maximum Entropy Principle (Shannon, 1948; Jaynes, 1957a,b). This principle consists in finding the probability density function that maximizes the entropy, given some available information. To compute the distribution, an optimization problem with constraints is solved, and it is guaranteed that all available information are respected and that no unphysical conditions appear (Kapur and Kesavan, 1992).

The available information is derived from the mechanical properties of the weight-on-hook. These properties are: (1) the column must penetrate the soil, i.e., *W* must be lower than the weight *w _{2}* of the column; (2) buckling must not occur, i.e.,

*W*must be greater than the buckling limit

*w*; (3) the probability must go to zero when

_{1}*W*approaches

*w*; (4) the probability must go to zero when

_{1}*W*approaches

*w*.

_{2}Conditions (1) and (2) are expressed by setting the support of the probability density function as ]*w _{1},w_{2}*[. Conditions (3) and (4) are expressed by

*E*{ln(

*W-w*)} = and,

_{1}*E*{ln(

*w*)} = with ||< +∞ and, ||< +∞, where

_{2}-W*E*{ } is the mathematical expectation. The reason why the logarithm, (ln), is used is because it imposes a weak decreasing of the probability density function in

*w*

_{1}^{+}and

*w*

_{2}^{-}. To facilitate the calculus, we introduce a normalized random variable X with values in ]0,1[, such that:

The expected value of *W* is written as

We introduce the notation = *E{W}* and = *E{X}*. The second moment of *W* may be written as

The available information is re-expressed in terms of the new random variable:

with |c_{1}|< +∞ and |c_{2}|< +∞. The optimization problem for the maximum Entropy Principle is finally written as:

where *C* is the space of admissible probability density functions *p _{X}* satisfying the constraints given by Eq. (23) and the entropy measure

*S*is given by (Shannon, 1948):

The probability density function, solution of the optimization problem defined by Eq. (24), is the Beta probability density function which may be written as:

where the Gamma function for *y* > 0. Also, *α* > 2 and *β* > 2 so that Eq. (23) holds. The random generator of independent realizations of the random variable *X* is already implemented in many computer codes. The mean value of *X* is given by

and coefficient of variation is given by

where *δ*_{X } = σ_{X}/, in which σ_{X} is the standard deviation.

The random weight-on-hook defined by Eq. (20) depends on the four parameters (*w*_{1}, *w*_{2}, α, β). For applications, *w*_{1} and *w*_{2} will be fixed. Parameters α and β have no physical meaning, consequently, we express them as function of the physical meaningful parameters *w* and δ. After some manipulations we obtain:

and

It should be noticed that with this scheme the random weight-on-hook can be computed for any fixed values of *w*_{1} and *w*_{2}. A specialist should be the one who gives these limits depending on the drill-string system analyzed. If the support ]*w*_{1},*w*_{2}[ is defined, and the mean and coefficient of variation are given, it is very simple to compute *W*.

**Stochastic Dynamical System **

Using the probabilistic model of the weight-on-hook, the deterministic reduced model defined by Eq. (18) is replaced by the following stochastic equations:

where **Q** is the random response and **F**_{W} is a vector for which the only nonzero component is related to the axial d.o.f. of the first node **F**_{W}(1) = (*W - *). Note that was subtracted because the response is calculated in the prestressed configuration.

**Numerical Results **

The data used in the simulations is found in the appendix. The drill-string is discretized with 56 finite elements, and for the construction of the reduced dynamical model, 10 torsional modes, 10 axial modes and also the two rigid body modes of the structure (axial and torsional) are used, hence 22 modes in the total for the reduced model. The time integration is done using an explicit Runge-Kutta algorithm with a time step controller to keep the error within a given accuracy.

**Convergence of the Stochastic Solution **

Let [**U**(t, s)] be the response of the stochastic dynamical system calculated for each realization *s*. The mean-square convergence analysis with respect to the number *n _{s}* of independent realizations is carried out studying the function conv(

*n*) defined by

_{s}Figure 2 shows that 500 simulations are sufficient to reach the mean-square convergence.

**Response of the Stochastic System **

The stochastic system's response is analyzed in this Section. Note that the dispersion of the response is all due to the random weight-on-hook, which is modeled as shown in Section four (Probabilistic model of the weight-on-hook). We may identify experimentally the parameters of the probabilistic model of the random weight-on-hook (mean and coefficient of variation) directly by measuring the weight-on-hook, or indirectly by measuring the dynamical response of the bit, for instance. In both cases, the identification procedure can be done using, for example, the Maximum Likelihood method (Aldrich, 1997; Spall, 2005). As we know that there are other sources of uncertainties for the problem, the experimentally identification procedure should be performed considering all the modeled random variables.

Figure 3 shows the 95% envelope (that is to say the confidence region constructed with a probability level of 0.95) for the rate-of-penetration and the rotational speed of the bit for a standard deviation σ = 1000 N, which means δ = σ/ is approximately 1x10^{-3}. The envelopes (the upper and lower envelopes of the confidence region) are calculated using the method of quantiles, Serfling (1980).

We are plotting two important variables: the rate-of-penetration (ROP) and the rotational speed at the bit (*ω*_{bit}). So, we analyze the influence of the random weight-on-hook in the system response. It can be seen that, for *σ* = 1000 N, the random response presents tie confidence intervals. Figure 4 shows the stochastic response of the torque and force on the bit.

It is noted that for *σ* = 1000 N, the response changes just a little, therefore, σ will be increased in the next analysis. In our analysis we can not increase σ too much, because the model used for the bit-rock interaction assumes a weight-on-bit *ƒ*_{bit} ~ -100 kN. Hence, the standard deviation σ of the *W* is increased in a way that the *ƒ*_{bit} has a maximum variation around 5%, that is to say that σ_{max} = 5000 N and, therefore, δ_{max} ~ 0.05 (0.5% variation), which is a constraint to our analysis. But, as it will be seen, a small variation on *W* may cause a big variation in the system response. Figures 5 and 6 show the system response for σ = 3000 N (δ ~ 0.003).

We want to see how uncertainties in the weight-on-hook affect the performance of the system; hence, Fig. 7 shows the evolution of the dispersion of the response for four dynamic responses: ROP, rotational speed of the bit, torque-on-bit, and force-on-bit. The dispersion of the response is calculated taking the square root of the variance divided by the value of the mean response for each time instant (it is the instant coefficient of variation). It can be noticed that the mean coefficient of variation of the force-on-bit (FOB) (see Fig. 7(d)) is about 3%, which is much more than 0.3%, which is the coefficient of variation of the random weight-on-hook *W*. This might be explained due to the fact that the absolute value of *W* is much higher than the absolute value of FOB; therefore, a small percentage variation of *W* has a great effect in the FOB, in terms of percent. What is more critical is the fact that the dispersion is even higher for the ROP and the rotational speed of the bit (Figs. 7(a) and (b)).

Figure 8 shows the system response for σ = 5000 N(δ ~ 0.005).

As expected, as δ increases the envelope of the response gets wider. Figure 9 shows the dispersion of the response for σ = 5000 N of the *W*.

It is noted that, even for a small variation of *W* (~ 0.5%), there is a big dispersion in the response. See, for instance, the rate-of-penetration: the mean coefficient of variation is 4.3%, which is more than eight times greater than the coefficient of variation of *W*. It gets worse if we take the maximum coefficient of variation, which is 16%. It means that if the *W* has a coefficient of variation of half percent, the variation in the ROP may achieve sixteen percent and the coefficient of variation of the rotational speed of the bit may achieve twenty six percent!

**Concluding Remarks **

A stochastic model of the drill-string dynamics has been analyzed. The weight-on-hook has been modeled as a random variable with probability density functions constructed using the Maximum Entropy Principle. It has been shown that the system response is sensitive to dispersion on the weight-on-hook. There are many sources of uncertainties in this problem; hence, more stochastic analysis should be done to identify the uncertainties that most affect the performance of the system.

**Acknowledgements **

The authors acknowledge the financial support of the Brazilian agencies CNPQ, CAPES, and FAPERJ.

**References **

Aldrich, J., 1997, "R.A. Fisher and the making of maximum likelihood 1912-1922", *Statistical Science*, Vol. 12, No. 3, pp. 162-176. [ Links ]

Jaynes, E., 1957a, "Information theory and statistical mechanics", *The Physical Review*, Vol. 106, No. 4, pp. 1620-630. [ Links ]

Jaynes, E., 1957b, "Information theory and statistical mechanics II", *The Physical Review*, Vol. 108, pp. 171-190. [ Links ]

Kapur, J.N., Kesavan, H.K., 1992, "Entropy Optimization Principles with Applications", Academic Press, Inc., USA. [ Links ]

Khulief, Y.A., AL-Naser, H., 2005, "Finite element dynamic analysis of drillstrings", *Finite Elements in Analysis and Design*, Vol. 41, pp. 1270-1288. [ Links ]

Kotsonis, S.J., Spanos, P.D., 1997, "Chaotic and random whirling motion of drillstrings", *Journal of Energy Resources Technology*, Vol. 119, No. 4, pp. 217-222. [ Links ]

Paidoussis, M.P., Luu, T.P., Prabhakar, S., 2008, "Dynamics of a long tubular cantilever conveying fluid downwards, which then flows upwards around the cantilever as a confined annular flow", *Journal of Fluids and Structures*, Vol. 24, pp. 111-128. [ Links ]

Piovan, M.T., Sampaio, R., 2009, "Modelos continuos de sondas de perforación para la industria petrolera: Análisis de enfoques y su discretización", *Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería*, Vol. 35, No. 3, pp. 259-277. [ Links ]

Ritto, T.G., Sampaio, R., Cataldo, E., 2008, "Timoshenko beam with uncertainty on the boundary conditions", *Journal of the Brazilian Society of Mechanical Sciences and Engineering*, Vol. 30, No. 4, pp. 295-303. [ Links ]

Ritto, T.G., Soize, C., Sampaio, R., 2009, "Nonlinear dynamics of a drillstring with uncertain model of the bit-rock interaction", *International Journal of Non-Linear Mechanics*, Vol. 44, No. 8, pp. 865-876. [ Links ]

Ritto, T.G., Soize, C., Sampaio, R., 2010, "Robust optimization of the rate of penetration of a drill-string using a stochastic nonlinear dynamical model", *Computational Mechanics*, article in press, DOI: 10.1007/s00466-009-0462-8. [ Links ]

Rubinstein, R.Y., 2007, "Simulation and the Monte Carlo Method", 2nd Edition, Series in Probability and Statistics, John Wiley and Sons, New Jersey, USA. [ Links ]

Sampaio, R., Piovan, M.T., Lozano, G.V., 2007, "Coupled axial/torsional vibrations of drilling-strings by mean of nonlinear model", *Mechanics Research Comunications*, Vol. 34, No. 5-6, pp. 497-502. [ Links ]

Sampaio, R., Ritto, T.G., 2008, "Short course on dynamics of flexible structures - deterministic and stochastic analysis", Seminar on Uncertainty Quantification and Stochastic Modeling, PUC-Rio, September. [ Links ]

Serfling, R.J., 1980, "Approximation Theorems of Mathematical Statistics", John Wiley and Sons, USA. [ Links ]

Shannon, C.E., 1948, "A mathematical theory of communication", *Bell System Tech. J.*, Vol. 27, pp. 379-423 and 623-659. [ Links ]

Soize, C., 2000, "A nonparametric model of random uncertainties for reduced matrix models in structural dynamics", *Probabilistic Engineering Mechanics*, Vol. 15, pp. 277-294. [ Links ]

Spall, J.C., 2005, "Introduction to Stochastic Search and Optimization", John Wiley and Sons, Hoboken, NJ, USA. [ Links ]

Spanos, P., Politis, N., Esteva, M., Payne, M., 2009, "Drillstring vibrations", *Advanced Drilling and Well Technology (Society of Petroleum Engineers)*, pp. 117-156. [ Links ]

Spanos, P.D., Payne, M., Secora, C., 1997, "Bottom-hole assembly modeling and dynamic response determination", *Journal of Energy Resources Technology*, Vol. 119, No. 3, pp. 153-158. [ Links ]

Trindade, M.A., Wolter, C., Sampaio, R., 2005, "Karhunen-Loève decomposition of coupled axial/bending of beams subjected to impacts", *Journal of Sound and Vibration*, Vol. 279, pp. 1015-1036. [ Links ]

Tucker, R.W., Wang, C., 1999, "An integrated model for drill-string dynamics", *Journal of Sound and Vibration*, Vol. 224, No. 1, pp. 123-165. [ Links ]

Tucker, R.W., Wang, C., 2003, "Torsional vibration control and Cosserat dynamics of a drill-rig assembly", *Meccanica*, Vol. 38, No. 1, pp. 143-159. [ Links ]

Paper accepted February, 2010.

Technical Editor: Domingos Alves Rade

**Appendix A - Data used in the simulations **

= | 1.06 x 10^{6} N (mean weight-on-hook), | |

w_{1} | = | 0.49 x 10^{6} N (weight-on-hook related to the buckling limit), |

w_{2} | = | 1.16 x 10^{6} N (weight of the column), |

L_{dp} | = | 1400 m (length of the drill pipe), |

L_{dc} | = | 200 m (length of the drill collar), |

D_{odp} | = | 0.127 m (outside diameter of the drill pipe), |

D_{odc} | = | 0.2286 m (outside diameter of the drill collar), |

D _{idp} | = | 0.095 m (inside diameter of the drill pipe), |

D_{idc} | = | 0.0762 m (inside diameter of the drill collar), |

E | = | 210 GPa (elasticity modulus of the drill string material), |

ρ | = | 7850 kg/m^{3} (density of the drill string material), |

Ω_{x} | = | 100 RPM (constant speed at the top), |

ρ_{f} | = | 1200 kg/m^{3} (density of the fluid), |

C_{f} | = | 1.25 x 10-2 (fluid viscous damping coefficient), |

ƒ_{c} | = | 100 kN (initial reaction force at the bit), |

g | = | 9.81 m/s^{2} (gravity acceleration), |

a_{1} | = | 3.429 x 10^{-3} m/s (constant of the bit-rock interaction model), |

a_{2} | = | 5.672 x 10_{-8} m/(N.s) (constant of the bit-rock interaction model), |

a_{3} | = | 1.374 x 10^{-4} m/rd (constant of the bit-rock interaction model), |

a_{4} | = | 9.537 x 10^{-6} N.rd (constant of the bit-rock interaction model), |

a_{5} | = | 1.475e x 10^{3} N.m (constant of the bit-rock interaction model), |

e | = | 2 rd/s (regularization parameter). |

The damping matrix is constructed using the relation [*C*] = *a*[*M*] + *b*([*K*] + [*K _{g}*(

**u**

_{S})]) with

*a*= 0.01 and

*b*= 0.0003.