This work intends to analyze the nonlinear stochastic dynamics of drillstrings in horizontal con g- uration. For this purpose, it considers a beam theory, with e ects of rotatory inertia and shear deformation, which is capable of reproducing the large displacements that the beam undergoes. The friction and shock e ects, due to beam/borehole wall transversal impacts, as well as the force and torque induced by bit-rock interaction, are also considered in the model. Uncertainties of bit- rock interaction model are taken into account using a parametric probabilistic approach. Numerical simula- tions have shown that the mechanical system of inter- est has a very rich nonlinear stochastic dynamics, which generate phenomena such as bit-bounce, stick-slip, and transverse impacts. A study aiming to maximize the drilling process eciency, varying drillstring velocities of translation and rotation is presented. Also, the work presents the de nition and solution of two optimizations problems, one deterministic and one robust, where the objective is to maximize drillstring rate of penetration into the soil respecting its structural limits. A. Cunha Jr (corresponding author) Universidade do Estado do Rio de Janeiro, Instituto de Matem atica e Estat stica, Departamento de Matem atica Aplicada, Rua S~ ao Francisco Xavier, 524, Pav. Jo~ ao Lyra, Bl. B, Sala 6032, Rio de Janeiro, 20550-900, Brasil E-mail: americo@ime.uerj.br A. Cunha Jr, C. Soize Universit e Paris-Est, Laboratoire Mod elisation et Simula- tion Multi Echelle, MSME UMR 8208 CNRS, 5, Boulevard Descartes 77454, Marne-la-Vall ee, France E-mail: christian.soize@univ-paris-est.fr A. Cunha Jr, R. Sampaio PUC{Rio, Departamento de Engenharia Mec^ anica, Rua M. de S~ ao Vicente, 225 - Rio de Janeiro, 22451-900, Brasil E-mail: rsampaio@puc-rio.brKeywords nonlinear dynamics horizontal drill- stringuncertainty quanti cation parametric probabilistic approach robust optimization 1 Introduction High energy demands of 21st century make that fossil fuels, like oil and shale gas, still have great importance in the energy matrix of several countries. Prospection of these fossil fuels demands the creation of exploratory wells. Traditionally, an exploratory well con guration is vertical, but directional or even horizontal con gura- tions, where the boreholes are drilled following a non- vertical way, are also possible [62]. An illustration of the di erent types of con gurations which an exploratory well can take is presented in Figure 1. reservoirvertical well directional well Fig. 1 Schematic representation of two exploratory wells. The left well con guration is vertical while the right one is directional. 2 A. Cunha Jr et al. The equipment used to drill the soil until the reser- voir level is called drillstring . This device is a long col- umn, composed of a sequence of connected drill-pipes and auxiliary equipment. Furthermore, within the col- umn ows drilling mud, which is used to cool the drilling system and to remove drilling cuttings from the bore- hole. The bottom part of this column is called bottom hole assembly (BHA) and consists of a pipe of greater thickness, named drill-colar , and a tool used to stick the rock, the drill-bit [20]. The BHA presents stabiliz- ers throughout its length, whose function is to maintain structural integrity of the borehole before cementation process. A schematic representation of a typical vertical drillstring and its components is presented in Figure 2, but a column in horizontal con guration essentially has the same structure. drill pipe drill colar drill bitstabilizer BHA Fig. 2 Schematic representation of a typical drillstring. Since the axial dimension of a drillstring is orders of magnitude larger than the characteristic dimension of its cross section area, the column is a long exible structure with a very complex exural dynamics. Fur- thermore, during drilling process, the drillstring is also subjected to other two mechanisms of vibration (longi- tudinal and torsional), which interact nonlinearly with the exural mechanism, resulting in a further compli- cated dynamics [60]. The coupling between these three mechanisms of vibration, which imposes severe com- plications on the drillstring dynamics modeling, comes from the action of several agents, such as: structure self weight (for a vertical column); tensile and compressive loads due to the weight on bit (WOB) and soil reac- tion force; dry friction and impacts with borehole wall; bit-rock interaction forces; internal ow pressure; forces induced by internal ow instabilities; etc [60]. The dynamics of a drillstring is not a new subject in the technical/scienti c literature. Works on this sub- ject, covering experimental analysis, numerical and/oranalytical modeling, can be seen since the 1960s. Most of the numerical works developed between 1960s and 1990s, have used lumped parameters approach to gain insight about drillstrings dynamical behavior. On the other hand, the analytical works focused on simple dis- tributed parameters models. Little has been done us- ing nite element-based approaches until the beginning of 1990s. A comprehensive literature survey of the re- search work produced until 2000 can be found in [10] and [60]. In recent studies, lumped parameters approach have been used, for example, to seek con gurations which re- duce stick-slip occurrence during drillstring operation [52]; to identify suitable values for drilling system op- erational parameters [32]; to analyze the coupling be- tween axial and torsional vibrations and its stability [19, 17, 34, 15]. On the other hand, approaches based on distributed parameters models have been used to: in- vestigate drillstring failure mechanisms [27]; better un- derstand transversal impacts between the column and borehole wall [61]; study the e ects induced by the non- linear coupling between longitudinal and torsional dy- namics of the drillstring [47]; describe the column dy- namic behavior taking into account the coupling be- tween the three mechanisms of vibration [41, 39]; inves- tigate the chaotic regime which drillstring transverse vibration mechanism is subjected [9]. Despite the fact that directional drilling has been used in practical engineering for a few decades, and most of the exploratory wells drilled today be direc- tional in con guration, all the works mentioned above model vertical drillstrings only. To the best of authors' knowledge, there are very few papers in open literature which models drillstring in directional con gurations [46, 24, 44]. All of these works use a distributed pa- rameters approach, but while [46, 44] only address the drillstring longitudinal dynamics, [24] uses generalized Euler-Bernoulli beam theory to describe the drillstring three-dimensional dynamics in a sloped directional well. In [46], the authors study a sloped con guration for the borehole and uses a perturbation technique to discretize the model equations. Conversely, model equations are discretized by nite element in [44]. In addition to the diculties inherent to the non- linear dynamics, drillstrings are subjected to random- ness on their geometrical dimensions, physical proper- ties, external forcing, etc. The lack of knowledge on these parameters, known as system-parameter uncer- tainty , is a source of inaccuracies in drillstring modeling, which may, in an extreme case, completely compromise the model predictability [49, 50]. Furthermore, during the modeling process, hypotheses about the drillstring physical behavior are made. These considerations may Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 3 be or not be in agreement with reality and should in- troduce additional inaccuracies in the model, known as model uncertainty induced by modeling errors [56, 57]. This source of uncertainty is essentially due to the use of simpli ed computational model for describing the phe- nomenon of interest and, usually, is the largest source of inaccuracy in computational model responses [56, 57]. Therefore, for a better understanding of the drill- string dynamics, these uncertainties must be modeled and quanti ed. In terms of quantifying these uncer- tainties for vertical drillstrings, the reader can see [59], where external forces are modeled as random objects and the method of statistical linearization is used along with the Monte Carlo (MC) method to treat the stochas- tic equations of the model. Other works in this line in- clude: [41, 39], where system-parameter and model un- certainties are considered using a nonparametric prob- abilistic approach; and [43, 40], which use a standard parametric probabilistic approach to take into account the uncertainties of the system parameters. Regarding the works that model directional con gurations, only [44] considers the uncertainties, which, in this case, are related to the friction e ects due to drillstring/borehole wall contact. From what is observed above, considering only the theoretical point of view, the study of drillstring nonlin- ear dynamics is already a rich subject. However, a good understanding of its dynamics also has signi cant im- portance in applications. For instance, it is fundamental to predict the fatigue life of the structure [33] and the drill-bit wear [68]; to analyze the structural integrity of an exploratory well [14]; to optimize the drill-bit rate of penetration (ROP) of into the soil [42], and the last is essential to reduce cost of production of an exploratory well. In this sense, this study aims to analyze the three- dimensional nonlinear dynamics of a drillstring in hor- izontal con guration, taking into account the system- parameter uncertainties. Through this study it is ex- pected to gain a better understanding of drillstring physics and, thus, to improve drilling process eciency, and maximize the column ROP accordingly. All results pre- sented here were developed in the thesis of [12]. The rest of this work is organized as follows. Sec- tion 2 presents the mechanical system of interest in this work, its parametrization and modeling from the physical point of view. Mathematical formulation of initial/boundary value problem that describes the me- chanical system behavior, as well as the conservative dynamics associated, is shown in section 3. The compu- tational modeling of the problem, which involves model equations discretization, reduction of the discretized dy- namics, the algorithms for numerical integration andsolution of nonlinear system of algebraic equations, can be seen in section 4. The probabilistic modeling of un- certainties is presented in section 5. Results of numeri- cal simulations are presented and discussed in section 6. Finally, in section 7, the main conclusions are empha- sized, and some paths to future works are pointed out. 2 Physical model for the problem 2.1 De nition of the mechanical system The mechanical system of interest in this work, which is schematically represented in Figure 3, consists of a hor- izontal rigid pipe, perpendicular to gravity, which con- tains in its interior a deformable tube under rotation. This deformable tube is subjected to three-dimensional displacements, which induces longitudinal, lateral, and torsional vibrations of the structure. These mechanisms of vibration are able to generate slips and shocks in ran- dom areas of the rigid tube. Also, the contact between drill-bit, at the right extreme of the tube, with soil gen- erates nonlinear forces and torques on drillstring right extreme, which may completely block the structure ad- vance over the well. Fig. 3 Schematic representation of the mechanical system under analysis. 2.2 Nonlinear dynamical system parameterization For purposes of modeling, the only part of column con- sidered is the BHA. So, any variation of diameter along the column is ignored. In this way, the bottom part of the deformable tube described, in section 2.1, is mod- eled as a rotating beam in horizontal con guration, whose transverse displacement ( yandz) at both ends is blocked, as well as transverse rotations on the left extreme. It looks like the left end of the system is a sta- bilizer and the right one a support. This beam is free to rotate around the xaxis, and to move longitudinally. The rigid pipe is treated as a stationary cylindrical rigid wall in horizontal con guration. As the beam is con ned within the borehole, it is reasonable to assume that it undergoes small rotations 4 A. Cunha Jr et al. in transverse directions. On the other hand, large dis- placements are observed in x,y, andz, as well as large rotations around the x-axis. Therefore, the analysis that follows uses a beam theory which assumes large rotation inx, large displacements in the three spatial directions, and small deformations [5]. Seeking not to make the mathematical model exces- sively complex, this work will not model the uid ow inside the beam, nor the dissipation e ects induced by the ow on the system dynamics. Due to the horizontal con guration, the beam is un- der action of the gravitational eld, which induces an acceleration g. The beam is made of an isotropic mate- rial with mass density , elastic modulus E, and Pois- son's ratio. It has length Land annular cross section, with internal radius Rintand external radius Rext. An illustration of beam geometric model is presented in Figure 4. It is important to note that this model also ignores the mass of drill-bit and its geometric shape. x Lyz Rext Rint Fig. 4 Schematic representation of the beam geometry used to model the deformable tube under rotation, and the inertial system of coordinates used. Using the cartesian coordinate system ( x;y;z ), de- ned by orthonormal basis fex;ey;ezg, xed in the in- ertial frame of reference R, and shown in Figure 4, one can describe the beam undeformed con guration by Bb=n (x;y;z )2R3 0xL;(y;z)2Sbo ;(1) where the beam cross section undeformed con guration is described by Sb=n (y;z)2R2 R2 inty2+z2R2 exto : (2) With the undeformed cross section con guration char- acterized, one can de ne and compute the cross-sectional area, A=ZZ Sbdydz = R2 extR2 int ; (3) the second moment of area around the yaxisIyy=ZZ Sbz2dydz =I4; (4) the second moment of area around the zaxis Izz=ZZ Sby2dydz =I4; (5) the polar moment of area Ixx=ZZ Sb y2+z2 dydz = 2I4; (6) the fourth moment of area around the zaxis Izzzz=ZZ Sby4dydz = 3I6; (7) and the fourth product of area Iyyzz=ZZ Sby2z2dydz =I6; (8) where I4= 4 R4 extR4 int ; (9) and I6= 24 R6 extR6 int : (10) In this work other three coordinate systems (all of then with the same origin as the ( x;y;z ) coordinate sys- tem) are also used, each one xed in a non-inertial frame of referenceRn, wheren= 1;2;3, and de ned by an orthonormal basis of vectors of the form fexn;eyn;ezng. These systems of coordinates are related by a se- quence of elementary rotations, such as follows Rx! R 1y! R 2z! R 3; (x;y;z ) (x1;y1;z1) (x2;y2;z2) (x3;y3;z3)(11) wherexis the rotation around the xaxis,yis the rotation around the yaxis, andzis the rotation around thezaxis. These rotations follow the right hand rule. Thus, with respect to the non-inertial frame of ref- erence, the instantaneous angular velocity of rotating beam is written as !=_xex+_yey1+_zez2; (12) Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 5 where _x,_y, and _zdenote the rate of rotation around thex,y, andzdirections, respectively. From now on, upper dot will be used as an abbreviation for time derivative, i.e., _@ @t. Referencing vector !to the inertial frame of refer- ence, and using the assumption of small rotations in transversal directions, one obtains !=0 B@_x+_zy _ycosx_zsinx _ysinx+_zcosx1 CA: (13) Regarding the kinematic hypothesis adopted for beam theory, it is assumed that the three-dimensional dis- placement of a beam point, occupying position ( x;y;z ) at instant of time t, can be written as ux(x;y;z;t ) =uyz+zy; (14) uy(x;y;z;t ) =v+y(cosx1)zsinx; uz(x;y;z;t ) =w+z(cosx1) +ysinx; whereux,uy, anduzrespectively denote the displace- ment of a beam point in x,y, andzdirections. More- over,u,v, andware the displacements of a beam neu- tral ber point in x,y, andzdirections, respectively. Finally, it is possible to de ne the vectors r=0 B@x y z1 CA;v=0 B@_u _v _w1 CA;and _=0 B@_x _y _z1 CA; (15) which, respectively, represent the position of a beam point, the velocity of a neutral ber point, and the rate of rotation of a neutral ber point. Note that the kinematic hypothesis of Eq.(14) is ex- pressed in terms of three spatial coordinates ( x,y, and z) and six eld variables ( u,v,w,x,y, andz). It is important to mention that, as the analysis as- sumed small rotations in yandz, this kinematic hy- pothesis presents nonlinearities, expressed by trigono- metric functions, only in x. Besides that, since a beam theory is employed, the eld variables in Eq.(14) depend only on the spatial coordinate xand timet. Therefore, although the kinematic hypothesis of Eq.(14) is three- dimensional, the mathematical model used to describe the beam nonlinear dynamics is one-dimensional. 2.3 Modeling of friction and shock e ects When a drillstring deforms laterally, there may occur a mechanical contact between the rotating beam andthe borehole wall, such as illustrated in Figure 5. This mechanical contact, which generally take place via a strong impact, gives rise to friction and shock e ects [21, 63, 31]. AA AA Fig. 5 Schematic representation of the situation where there is a mechanical contact between a drillstring and the borehole wall. The modeling of friction and shock phenomena is made in terms of a geometric parameter dubbed inden- tation , which is de ned as FS=rgap; (16) wherer=p v2+w2is the neutral ber lateral dis- placement, and gapdenotes the spacing between the undeformed beam and the borehole wall. One has that FS>0 in case of an impact, or FS0 otherwise, as can be seen in Figure 6. Note that the indentation cor- responds to a measure of penetration in the wall of a beam cross section [21]. gap rgap r δFS=r−gap≤0 δFS=r−gap>0 Fig. 6 Illustration of the indentation parameter in a situa- tion without impact (left) or with impact (right). When an impact occurs, a normal force of the form Fn FS=kFS1FSkFS23 FScFSjj3_FS; (17) wherekFS1,kFS2andcFSare constants of the shock model, begins to act on the beam cross section. In this nonlinear shock model, proposed by Hunt and Crossley [26], the rst (linear spring) and the second (nonlinear spring) terms describe the elastic deformation during an impact, while the third term (nonlinear damper) takes into account the loss of energy during an impact. 6 A. Cunha Jr et al. Once the column is rotating and moving axially, an impact also induces a frictional force in axial direction, Fa FS, and a torsional friction torque, TFS. Both are mod- eled by Coulomb friction law [11], so that the force is given by Fa FS=FSFn FSsgn( _u); (18) where the torque is described by TFS=FSFn FSRbhsgn _x ; (19) beingFSthe friction coecient, sgn() the sign func- tion, and the radius of the borehole is Rbh=Rext+gap. In order to nd all points of contact between the beam and the borehole wall, it is necessary to discover all values of xwhereFS>0. This is usually done by solving an optimization problem with constraints [65, 64]. Although the strategy of detection based on the op- timization problem is robust in terms of accuracy, it is extremely complex in terms of implementation and computational cost. For this reason, this work uses an approach that introduces the forces of Eqs.(17) and (18), and the torque of Eq.(19), as e orts concentrated on the nodes of nite element mesh, de ned in the sec- tion 4.1. This procedure sacri ces accuracy, but simpli- es the friction and shock model implementation. 2.4 Modeling of bit-rock interaction e ects During the drilling process, in response to drillstring rotational advance, a force and a torque of reaction be- gin to act on the drill-bit, giving rise to the so-called bit-rock interaction e ects [16, 18]. In this work, the model proposed by [44] is consid- ered to describe the bit-rock interaction force FBR=8 < :BR e BR_ubit1 for _ubit>0, (20) 0 for _ ubit0, whereBRis the bit-rock limit force, BRis the rate of change of bit-rock force, and _ ubit= _u(L;). The graph ofFBRis illustrated in Figure 7. Also, for the bit-rock interaction torque, it is adopted the regularized Coulomb model used by [28], which is expressed as TBR=BRFBRRbhBR(!bit); (21) ˙ubitFBR −ΓBRFig. 7 Illustration of the function used to describe the reac- tion force on drill-bit, due to bit-rock interaction e ects. where!bit=_x(L;), the bit-rock friction coecient is BR, and BR(!bit) = tanh (!bit) +2!bit 1 +!2 bit; (22) is a regularization function. The graph of BRis illus- trated in Figure 8. ωbitξBR Fig. 8 Illustration of the smooth function used to regularize the reaction torque on drill-bit, due to bit-rock interaction e ects. 2.5 Kinetic energy The kinetic energy of rotating beam is given by T=1 2ZZZ Bbvvdxdydz + (23) 1 2ZZZ Bb!(rrIr r)!dxdydz; where rst triple integral corresponds to beam transla- tional kinetic energy, and the second one is associated to beam rotational kinetic energy. In this equation, I denotes the identity tensor, the symbol represents the standard inner product between two Euclidean vectors, Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 7 and the symbol is used to designate the tensor prod- uct. Developing the vector operations indicated in Eq.(23), using (1) and (2) to de ne the limits of integration, us- ing the de nitions of A,Iyy,Izz, andIxx, and mak- ing the other calculations one can show that Eq.(23) is equivalent to T=1 2ZL x=0A _u2+ _v2+ _w2 dx+ (24) 1 2ZL x=02I4 _x+_zy2 dx+ 1 2ZL x=0I4 _ycosx_zsinx2 dx+ 1 2ZL x=0I4 _ysinx+_zcosx2 dx: 2.6 Strain energy The analysis of the beam assumes that it is subjected to large displacements, and small deformations. In this way, its strain energy is given by V=1 2ZZZ Bb:dxdydz; (25) wheredenotes the Green-Lagrangian strain tensor, is the second Piola-Kirchho stress tensor, and the symbol :represents the double inner product between two tensors. It is further considered that the beam is made of an isotropic material, such that stress and strain are related by the following constitutive equation (Hooke's law) = 2G+tr()I; (26) where tr() represents the trace operator, Gis material shear modulus, and is used to designate the material rst Lam e parameter. In terms of the elastic modulus Eand the Poisson's ratio , these elastic parameters can be written as G=E 2 (1 +); and=E (1 +)(12):(27) According to the beam theory used in this work, there is no tension in any cross section of the beam that is perpendicular to xaxis, i.e.,yy= 0,zz= 0, yz= 0, andzy= 0. When this hypothesis is com- bined with the three-dimensional Hooke's law, repre- sented by Eq.(26), one can conclude that xx=Exx,xy= 2Gxy, andxz= 2Gxz, which is an one- dimensional version of Hooke's law. Combining this one-dimensional Hooke's law with the stress tensor symmetry, one can express the double contraction between strain and stress tensors, within the integral in Eq.(25), as a quadratic form :=E2 xx+ 4G2 xy+ 4G2 xz; (28) which is modi ed, by the introduction of shearing factor s, as :=E2 xx+ 4sG2 xy+ 4sG2 xz: (29) This modi cation aims to take into account the e ect of shear deformation in the beam cross section area, which is neglected when one uses the one-dimensional Hooke's law. Hence, after replace Eq.(29) in Eq.(25), one nally obtains V=1 2ZZZ Bb E2 xx+ 4sG2 xy+ 4sG2 xz dxdydz: (30) As the analysis is using large displacements, one has xx=1 2@ux @x+@ux @x + (31) 1 2@ux @x@ux @x+@uy @x@uy @x+@uz @x@uz @x ; xy=1 2@uy @x+@ux @y + (32) 1 2@ux @x@ux @y+@uy @x@uy @y+@uz @x@uz @y ; and xz=1 2@uz @x+@ux @z + (33) 1 2@ux @x@ux @z+@uy @x@uy @z+@uz @x@uz @z ; where the quadratic terms on the right hand side of the above equations are associated to beam model geomet- ric nonlinearity. Substituting the kinematic hypothesis of Eq.(14) in Eqs.(31) to (33), and then calculating the partial deriva- tives, one concludes that the deformations are respec- tively given by 8 A. Cunha Jr et al. xx=u0y0 z+z0 y+u0 z0 yy0 z yz0 y0 z+ 0 x yw0zv0 cosx yv0+zw0 sinx + 1 2 u02+v02+w02+y202 z+z202 y+ y2+z2 02 x ;(34) xy=1 2 v0cosx+w0sinxz0 x + (35) 1 2z y0 zz0 yu01 ; and xz=1 2 w0cosxv0sinx+y0 x + (36) 1 2y y0 z+z0 y+u0+ 1 ; where prime is used as an abbreviation for space deriva- tive, i.e., 0@ @x. 2.7 Energy dissipation function It is assumed that the beam loses energy through a mechanism of viscous dissipation, with a (dimension- less) damping constant c. In this way, there is an energy dissipation function (per unit of length) associated to the system, which is given by D=1 2ZZ Sbcvvdydz + (37) 1 2ZZ Sbc_(rrIr r)_dydz; where the rst term is a dissipation potential due to the translational movement, and the second term represents a dissipation potential due to the movement of rotation. Making a development almost similar to the one per- formed to obtain Eq.(24), it can be shown that D=1 2cA _u2+ _v2+ _w2 + (38) 1 2cI 4 2_2 x+_2 y+_2 z : 2.8 External forces work The work done by the external forces acting on the beam is given by W=ZL x=0Agwdx +WFS+WBR; (39)where the rst term is due to gravity, the second one is associated to friction and shock e ects, and the last term accounts the work done by the force/torque that comes from bit-rock interaction. Note that, due to non-holonomic nature of the forces and torques that comes from the e ects of friction/shock, and bit-rock interaction, it is not possible to write ex- plicit formulas for WFSandWFS[30]. However, it is known that the virtual work of WFS, denoted by WFS, is written as WFS=NnodesX m=1 Fa FSu+Fn FS(vv+ww )=r+TFSx x=xm(40) wherexmare the global coordinates of nite element nodes,Nnodes is the number of nodes in nite element mesh, andu,v,w, andxrespectively denote the variations of u,v,w, andx. On the other hand, the virtual work of WBR, denoted byWBR, reads as WBR=FBRu x=L+TBRx x=L: (41) 3 Mathematical model for the problem 3.1 Equation of motion of the nonlinear dynamics A modi ed version of Hamilton's extended principle [30] is employed to derive the equations which describe the mechanical system nonlinear dynamics, so that the rst variation is expressed as Ztf t=t0(TV+W)dtZtf t=t0ZL x=0U@D @_Udxdt = 0;(42) where the rst term corresponds to dynamics conser- vative part, and the second one is associated to energy dissipation. Also, Uis a vector eld which lumps the eld variables, the initial and nal instants of obser- vation are respectively denoted by t0andtf, and the symbolrepresents the variation operator [45]. The development of Eq.(42) results in the following weak equation of motion M ;U +C ;_U +K( ;U) =F ;U;_U;U ;(43) valid for any chosen in a \suitable" space of weight functions, where the eld variables and their correspond- ing weight functions are represented by the vector elds U= u;v;w;x;y;z , and = u; v; w; x; y; z . Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 9 Furthermore, M ;U =ZL x=0A( uu+ vv+ ww)dx+ (44) ZL x=0I4 2 xx+ yy+ zz dx; represents the mass operator, C ;_U =ZL x=0cA ( u_u+ v_v+ w_w)dx+ (45) ZL x=0cI 4 2 x_x+ y_y+ z_z dx; is the damping operator, K( ;U) =ZL x=0EA 0 uu0dx+ (46) ZL x=0EI4 0 y0 y+ 0 z0 z dx+ ZL x=02sGI4 0 x0 xdx+ ZL x=0sGA y+ 0 w y+w0 dx+ ZL x=0sGA z 0 v zv0 dx; is the sti ness operator, and F ;U;_U;U =FKE ;U;_U;U + (47) FSE( ;U) +FFS( ;U) + FBR ;_U +FG( ); is the force operator, which is divided into ve parts. A nonlinear force due to inertial e ects FKE=ZL x=02I4 x yz+_y_z dx (48) +ZL x=02I4 y y_2 z+_x_z dx ZL x=02I4 z yx+2 yz dx ZL x=02I4 z _x_y+ 2y_y_z dx; a nonlinear force due to geometric nonlinearityFSE=ZL x=0 x1+ y2+ z3 dx+ (49) ZL x=0 0 u4+ 0 v5+ 0 w6 dx+ ZL x=0 0 x7+ 0 y8+ 0 z9 dx; a nonlinear force due to friction and shock e ects FFS=NnodesX m=1 Fa FS u+Fn FS(v v+w w)=r+TFS x x=xm(50) a nonlinear force due to bit-rock interaction FBR=FBR u x=L+TBR x x=L; (51) and a linear force due to gravity FG=ZL x=0Ag wdx: (52) The nonlinear functions n, withn= 1;;9, in Eq.(49) are very complex and, for sake of space limita- tion, are not presented in this section. But they can be seen in Appendix A. The model presented above is an adaptation, for horizontal drillstrings, of the model proposed by [41, 39] to describe the nonlinear dynamics of vertical drill- strings. To be more precise, in the reference problem gravity is parallel to drillstring main axis, while in this work, it is perpendicular to the structure primal direc- tion. Therefore, the former problem primarily addresses the dynamics of a column, while the new problem deals with the dynamics of a beam. Also, the original problem treated the nonlinear dynamics around a pre-stressed equilibrium con guration, while the new problem does not consider the dynamics around any particular con- guration. It is worth mentioning that changes made in the modeling of friction and shock e ects are signi cant. For instance, a nonlinear shock model that also takes into account the dissipation of energy during an impact is introduced, in contrast to the reference work, that only consider the linear elastic deformation e ects. In addition, the boundary conditions are di erent, as well as the bit-rock interaction model. On the other hand, for the sake of simplicity, the uid structure interac- tion e ects, considered in reference work are neglected in this study. 10 A. Cunha Jr et al. 3.2 Initial conditions With regard to mechanical system initial state, it is assumed that the beam presents neither displacement nor rotations, i.e., u(x;0) = 0,v(x;0) = 0,w(x;0) = 0, x(x;0) = 0,y(x;0) = 0, and z(x;0) = 0. These eld variables, except for uandx, also have initial velocities and rate of rotations equal to zero, i.e. _ v(x;0) = 0, _w(x;0) = 0, _y(x;0) = 0, and _z(x;0) = 0. It is also assumed that, initially, the beam moves horizontally with a constant axial velocity V0, and ro- tates around the xaxis with a constant angular velocity . Thereby, one has that _ u(x;0) =V0, and _x(x;0) = . Projecting the initial conditions in a\suitable"space of weight functions, weak forms for them are obtained, respectively, given by M ;U(0) =M( ;U0); (53) and M ;_U(0) =M ;_U0 ; (54) whereU0= (0;0;0;0;0;0) and _U0= (V0;0;0; ;0;0). In formal terms, the weak formulation of the initial{ boundary value problem that describes the mechanical system nonlinear dynamics consists in nd a vector eld U, \suciently regular", which satis es the weak equa- tion of motion given by Eq.(43) for all \suitable" , as well as the weak form of initial conditions, given by Eqs.(53), and (54) [25]. 3.3 Associated linear conservative dynamics Consider the linear homogeneous equation given by M ;U +K( ;U) = 0; (55) obtained from Eq.(43) when one discards the damping, and force operators, and which is valid for all in the space of weight functions. Suppose that Eq.(55) has a solution of the form U=ei!t , where!is a natural frequency (in rad/s), is the associated normal mode, and i=p 1 is the imaginary unit. Replacing this expression of Uin Eq.(55) and using the linearity of operators M, andK, one gets  !2M( ; ) +K( ; ) ei!t= 0; (56) which is equivalent to!2M( ; ) +K( ; ) = 0; (57) a generalized eigenvalue problem. Since operatorMis positive-de nite, and operator Kis positive semi-de nite, the generalized eigenvalue problem above has a denumerable number of solutions. The solutions of this eigenvalue problem have the form (!2 n; n), where!nis then-th natural frequency and nis then-th normal mode [23]. Also, the symmetry of operators M, andKimplies the following orthogonality relations M( n; m) =nm;andK( n; m) =!2 nnm;(58) wherenmrepresents the Kronecker delta symbol. See [23] for more details. The generalized eigenvalue problem of Eq.(57), as well as the properties of (58), will be useful to construct a reduced order model for discretized dynamical system which approximates the solution of the weak initial{ boundary value problem of Eqs.(43), (53), and (54). 4 Computational model for the problem 4.1 Discretization of the nonlinear dynamics To proceed with the discretization of the weak initial{ boundary value problem of Eqs.(43), (53), and (54), which describes the rotating beam nonlinear dynamics, it is used the standard nite element method (FEM) [25], where the spaces of basis and weight functions are constructed by the same ( nite dimensional) class of functions. In this procedure, the beam geometry is discretized by a FEM mesh with Nelem nite elements. Each one of these elements is composed by two nodes, and each one of these nodes has six degrees of freedom associated, one for each eld variable in the beam model described in section 3.1. Thus, the number of degrees of freedom associated with FEM model is Ndofs = 6(Nelem+ 1). An illustration of FEM mesh/element can be seen in Figure 9. Concerning the shape functions, it is adopted an in- terdependent interpolation scheme which avoids shear- locking e ect [38]. This scheme uses, for transverse dis- placements/rotations, Hermite cubic polynomials, and, for the elds of axial displacement/torsional rotation, ane functions [2]. Thus, each eld variable of the physical model is ap- proximated by a linear combination of basis functions, in such way that Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 11 u1v1w1 θx1θy1θz1 u2v2w2 θx2θy2θz2 Fig. 9 Illustration of FEM mesh/element used to discretize the beam geometry. u(x;t)NdofsX m=1Qm(t)Nm(x); (59) x(x;t)NdofsX m=1Qm(t)Nm(x); v(x;t)NdofsX m=1Qm(t)H(1) m(x); w(x;t)NdofsX m=1Qm(t)H(1) m(x); y(x;t)NdofsX m=1Qm(t)H(2) m(x); z(x;t)NdofsX m=1Qm(t)H(2) m(x); whereNm(x),H(1) m(x), andH(2) m(x) are the (position dependent) shape functions, and Qm(t) are the (time dependent) coecients of the linear combination. In physical terms, each one of these temporal coecients represents a degree of freedom of the FEM model. The discretization results is the NdofsNdofsnon- linear system of ordinary di erential equations given by [M]Q(t) + [C]_Q(t) + [K]Q(t) =F Q;_Q;Q ;(60) whereQ(t) is the nodal displacement vector (transla- tions and rotations), _Q(t) is the nodal velocity vector, andQ(t) is the nodal acceleration vector. The other ob- jects in Eq.(60) are the mass matrix [ M], the damping matrix [C], the sti ness matrix [ K], and the force vector F. A discretization procedure similar to one presented above is applied to the initial conditions of Eqs.(53) and (54), which results in linear systems of algebraic equations given by[M]Q(0) =Q0; and [M]_Q(0) = _Q0: (61) 4.2 Reduction of nite element model In order to reduce the dimension of nite element model developed in section 4.1, it is considered a nite dimen- sional version of the generalized eigenvalue problem pre- sented in section 3.3, which is de ned by [K] n=!2 n[M] n: (62) Due to the properties of M, andKoperators, dis- cussed in section 3.3, and inherited by nite dimensional operators [M] and [K], the above eigenvalue problem hasNdofssolutions. But Eq.(62) is solved only for n= 1;2;;Nred, where the reduced model dimension Nred is an integer chosen such that Nred Ndofs. The procedure that follows consists in project the nonlinear dynamics, de ned by the initial value problem of Eqs.(60) and (61), into the vector space spanned by f 1; 2;; Nredg. For this purpose, de ne the NdofsNredprojection matrix by [] =2 64j j j 1 2 Nred j j j3 75; (63) make in Eqs.(60) and (61) the change of basis de ned by Q(t) = []q(t); (64) and then pre-multiply the resulting equations by matrix []T, where superscriptTrepresents the transposition operation. This development results in the reduced initial value problem given by [M]q(t)+[C]_q(t)+[K]q(t) =f q(t);_q(t);q(t) ;(65) and q(0) =q0; and _q(0) = _q0; (66) whereq(t) is the reduced displacement vector, _q(t) is the reduced velocity vector, q(t) is the reduced acceler- ation vector. The reduced matrices of mass, damping, 12 A. Cunha Jr et al. and sti ness, as well as the reduced vectors of force, ini- tial displacement, and initial velocity are, respectively, de ned by [ M] = []T[M] [], [C] = []T[C] [], [K] = []T[K] [],f= []TF []q(t);[]_q(t);[]q(t) ,q0= []TQ0,_q0= []T_Q0. These matrices are NredNred, while these vectors are Nred1. Furthermore, due to the orthogonality properties de ned by Eq.(58), that are inherited by the operators in nite dimension, these matrices are diagonal. Thus, although the initial value problem of Eqs.(65) and (66) is apparently similar to the one de ned by Eqs.(60) and (61), the former has a structure that makes it much more ecient in terms of computational cost, and so, it will be used to analyze the nonlinear dynam- ics under study. 4.3 Integration of discretized nonlinear dynamics In order to solve the initial value problem of Eqs.(65) and (66), it is employed the Newmark method [35], which de nes the following implicit integration scheme _qn+1=_qn+ (1 )tqn+ tqn+1; (67) qn+1=qn+t_qn+1 2  t2qn+ t2qn+1;(68) whereqn,_qnandqnare approximations to q(tn),_q(tn) andq(tn), respectively, and tn=nt is an instant in a temporal mesh de ned over the interval [ t0;tf], with an uniform time step t. The parameters and are associated with accuracy and stability of the nu- merical scheme [25], and for the simulations reported in this work they are assumed as = 1=2 + , and = 1=4 1=2 + 2, with = 15=1000. Handling up properly Eqs.(67) and (68), and the dis- crete version of Eq.(65), one arrives in a nonlinear sys- tem of algebraic equations, with unknown vector qn+1, which is represented by ^[K]qn+1=^fn+1(qn+1); (69) where^[K] is the e ective sti ness matrix, and ^fn+1is the (nonlinear) e ective force vector. 4.4 Incorporation of boundary conditions As can be seen in Figure 4, the mechanical system has the following boundary conditions: (i) left extreme with no transversal displacement, nor transversal rotation;(ii) right extreme with no transversal displacement. It is also assumed that the left end has: (iii) constant axial and rotational velocities in x, respectively equal to V0 and . Hence, for x= 0, it is true that u(0;t) =V0t, v(0;t) = 0,w(0;t) = 0,x(0;t) = t,y(0;t) = 0, andz(0;t) = 0. On the other hand, for x=L, one has v(L;t) = 0, andw(L;t) = 0. The variational formulation presented in section 3.1, was made for a free-free beam, i.e. the above geometric boundary conditions were not considered. For this rea- son, they are included in the formulation as constraints, using the Lagrange multipliers method [25]. The details of this procedure are presented below. Observe that the boundary conditions can be rewrit- ten in matrix form as [B]Q(t) =h(t); (70) where the constraint matrix [ B] is 8Ndofs and has almost all entries equal to zero. The exceptions are [B]ii= 1 fori=f1;;6g, [B]7(Ndofs5)= 1, and [B]8(Ndofs4)= 1. The constraint vector is given by h(t) =0 BBBBBBBBBBB@u(0;t) v(0;t) w(0;t) x(0;t) y(0;t) z(0;t) v(L;t) w(L;t)1 CCCCCCCCCCCA: (71) Making the change of basis de ned by Eq.(64), one can rewrite Eq.(70) as [B]q(t) =h(t); (72) where the 8Nredreduced constraint matrix is de ned by [B] = [B] []. The discretization of Eq.(72) results in [B]qn+1=hn+1; (73) wherehn+1is an approximation to h(tn+1). This equa- tion de nes the constraint that must be satis ed by the variational problem \approximate solution". In what follows it is helpful to think that Eq.(69) comes from the minimization of an energy functional qn+17!F(qn+1), which is the weak form of this non- linear system of algebraic equations. Then, one de nes the Lagrangian as Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 13 L(qn+1;n+1) =F(qn+1) +T n+1 [B]qn+1hn+1 ;(74) being the (time-dependent) Lagrange multipliers vector of the form n+1=0 BBBBBBBBBBB@1(tn+1) 2(tn+1) 3(tn+1) 4(tn+1) 5(tn+1) 6(tn+1) 7(tn+1) 8(tn+1)1 CCCCCCCCCCCA: (75) Invoking the Lagrangian stationarity condition one arrives in the following ( Nred+ 8)(Nred+ 8) system of nonlinear algebraic equations " ^[K] [B]T [B] [0]# qn+1 n+1! = ^fn+1 hn+1! ; (76) where [0] is a 8 8 null matrix. The unknowns are qn+1andn+1, and must be solved for each instant of time in the temporal mesh, in order to construct an approximation to the mechanical system dynamic response. The solution of the nonlinear system of algebraic equations, de ned by Eq.(76), is carried out rst ob- taining and solving a discrete Poisson equation for n+1 [22], and then using the rst line of (76) to obtain qn+1. To solve these equations, a procedure of xed point iter- ation is used in combination with a process of successive over relaxation [66]. 5 Probabilistic modeling of system-parameter uncertainties The mathematical model used to describe the physical behavior of the mechanical system is an abstraction of reality, and its use does not consider some aspects of the problem physics. Regarding the system modeling, either the beam theory used to describe the structure dynamics [41], as the friction and shock model used [26] are fairly established physical models, who have gone through several experimental tests to prove their validity, and have been used for many years in similar situations. On the other hand, the bit-rock interaction model adopted in this work, until now was used only in a purely numeric context [44], without any experi- mental validation. Thus, it is natural to conclude thatbit-rock interaction law is the weakness of the model proposed in this work. In this sense, this work will focus on modeling and quantifying the uncertainties that are introduced in the mechanical system by bit-rock interaction model. For convenience, it was chosen to use a parametric proba- bilistic approach [56], where only the uncertainties of system parameters are considered, and the maximum entropy principle is employed to construct the proba- bility distributions. 5.1 Probabilistic framework LetXbe a real-valued random variable, de ned on a probability space ( ;;P), for which the probability distribution PX(dx) onRadmits a density x7!pX(x) with respect to dx. The support of the probability den- sity function (PDF) pXwill be denoted by SuppXR. The mathematical expectation of Xis de ned by E[X] =Z SuppXxpX(x)dx; (77) and any realization of random variable Xwill be de- noted by X() for2. LetmX=E[X] be the mean value,2 X=Eh (XmX)2i be the variance, and X= q 2 Xbe the standard deviation of X. The Shannon en- tropy ofpXis de ned by S(pX) =E lnpX(X) . 5.2 Probabilistic model for bit-rock interface law Recalling that bit-rock interaction force and torque are, respectively, given by Eqs.(20) and (21), the reader can see that this bit-rock interface law is characterized by three parameters, namely, BR,BR, andBR. The con- struction of the probabilistic model for each one param- eter of these parameters, which are respectively mod- eled by random variables BR,BR, andBR, is pre- sented below. 5.3 Distribution of force rate of change As the rate of change BRis positive, it is reasonable to assume Supp BR=]0;1[. Therefore, the PDF of BRis a nonnegative function p BR, such that Z+1 =0p BR( )d = 1: (78) 14 A. Cunha Jr et al. It is also convenient to assume that the mean value of BRis a known positive number, denoted by m BR, i.e., E[ BR] =m BR>0: (79) One also need to require that E ln ( BR) =q BR;jq BRj<+1; (80) which ensures, as can be seen in [53, 54, 55], that the inverse of BRis second order random variable. This condition is necessary to guarantee that the stochastic dynamical system associated to this random variable is of second order, i.e., it has nite variance. Employing the principle of maximum entropy one need to maxi- mize the entropy function S(p BR), respecting the con- straints imposed by (78), (79) and (80). The desired PDF corresponds to the gamma distri- bution and is given by p BR( ) =1]0;1[( )1 m BR 1 2 BR!1=2 BR 1 (1=2 BR) m BR1=2 BR1 exp  2 BRm BR! ; (81) where the symbol 1]0;1[( ) denotes the indicator func- tion of the interval ]0 ;1[, 0 BR= BR=m BR< 1=p 2 is a type of dispersion parameter, and (z) =Z+1 y=0yz1eydy; (82) is the gamma function. 5.4 Distribution of limit force The parameter BRis also positive, in a way that SuppBR=]0;1[, and consequently Z+1 =0pBR( )d = 1: (83) The hypothesis that the mean is a known positive numbermBRis also done, i.e., E[BR] =mBR>0; (84)as well as that the technical condition, required for the stochastic dynamical system associated be of second or- der, is ful lled, i.e. E ln (BR) =qBR;jqBRj<+1: (85) In a similar way to the procedure presented in sec- tion 5.3, it can be shown that PDF of maximum entropy is also gamma distributed, and given by pBR( ) =1]0;1[( )1 mBR 1 2 BR!1=2 BR 1 (1=2 BR) mBR1=2 BR1 exp  2 BRmBR! : (86) 5.5 Distribution of friction coecient With respect to the parameter BR, one know it is non- negative and bounded above by the unity. Thus, one can safely assume that SuppBR= [0;1], so that the normalization condition read as Z1 =0pBR()d= 1: (87) The following two conditions are also imposed E ln (BR) =q1 BR;jq1 BRj<+1; (88) E ln (1BR) =q2 BR;jq2 BRj<+1; (89) representing a weak decay of the PDF of BRin 0+and 1respectively [53, 54, 55]. Evoking again the principle of maximum entropy considering now as known infor- mation the constraints de ned by (87), (88), and (89) one has that the desired PDF is given by pBR() =1[0;1]()(a+b) (a)(b)a1(1)b1;(90) which corresponds to the beta distribution The parameters aandbare associated with the shape of the probability distribution, and can be re- lated withmBRandBRby a=mBR 2BR 1 mBR2 BR1! ; (91) Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 15 and b=mBR 2BR 1 mBR2 BR1! 1 mBR1! :(92) 5.6 Stochastic nonlinear dynamical system Due to the randomness of parameters BR,BR, and BR, the mechanical system physical behavior is now described, for all in, by the stochastic nonlinear dynamical system de ned by [M] q(t;) + [C] _q(t;) + [K]q(t;) =f(q;_q;q);(93) q(0;) =q0; and _q(0;) =_q0; a:s: (94) whereq(t) is the random reduced displacement vector, _q(t) is the random reduced velocity vector, and  q(t) is the random reduced acceleration vector, and fis the random reduced nonlinear force vector. The methodology used to calculate the propagation of uncertainties through this stochastic dynamical sys- tem is Monte Carlo (MC) method [29], employing a strategy of parallelization described in [13]. 6 Numerical experiments and discussions In order to simulate the mechanical system nonlinear dynamics, the physical parameters presented in Table 1 are adopted, as well as L= 100m, the rotational and axial velocities in x, respectively given by = 2 rad=s , andV0= 1=180m=s. The values of these parameters do not correspond exactly to actual values used in a real drillstring, but are of the same order of magnitude. For this con guration, the beam geometry is discretized by 500 nite elements, and the interval of integration [ t0;tf] = [0;10]sis considered. For the friction and shock model constants, are con- sidered the values shown in Table 2, which have order of magnitude typical of a borehole wall made of steel [67].The low value for FSis justi ed by the fact that, in a real system, there is a uid between the column and borehole wall, which carries a substantial reduction in the torsional friction. The bit-rock interaction model constants can be seen in Table 3, and were estimated following a similar strat- egy as that shown in [44].Table 1 Physical parameters of the mechanical system that are used in the simulation. parameter value unit  7900 kg=m3 g 9:81m=s2  0:3 | s 6/7 | c 0:01 | E 203109Pa Rbh 95103m Rint 50103m Rext 80103m Table 2 Parameters of the friction and shock model that are used in the simulation. parameter value unit kFS1 11010N=m kFS2 11016N=m3 cFS 1106(N=m3)=(m=s) FS 0:25 | Table 3 Parameters of the bit-rock interaction model that are used in the simulation. parameter value unit BR 30103N BR 400 1 =(m=s) BR 0:4 | 6.1 \Validation" of the computational model In order to validate to computational model developed in this work, the most widely accepted approach would be through comparison with experimental data [36]. Unfortunately, it is very dicult to obtain such experi- mental measurements for the type of mechanical system analyzed in this work. Accordingly, as in many ana- lyzes of dynamic systems with complex behavior (e.g. aerospace, automotive, etc), this work used numerical simulation without experimental analysis. However, it is clear that such numerical simulations do have meaning, and are only useful, if prior validation of models, for- mulations, and numerical simulation software were ob- tained. In the cases treated, \validations" of the model, formulation, and software were obtained, step by step, by comparing the results obtained with the software with the following reference results: (i) linear static; (ii) nonlinear static with contact; (iii) linear dynamic nonrotating without contact; (iv) linear dynamic rotat- ing without contact; (v) nonlinear nonrotating dynamic with contacts (shocks); and (vi) nonlinear rotating dy- namic with contacts (shocks). All tests were performed and the analyzes were used to \validate" the model, the formulation, and the results of numerical simulations. 16 A. Cunha Jr et al. 6.2 Modal analysis of the mechanical system In this section, the modal content of the mechanical sys- tem is investigated. This investigation aims to identify the natural frequencies of the system, and, especially, to check the in uence of slenderness ratio , de ned as the ratio between beam length and external diameter, in the natural frequencies distribution. Therefore, the dimensionless frequency band for the problem is assumed as being B= [0;4], with the dimen- sionless frequency de ned by f=fL cL; (95) wherefis the dimensional frequency (Hz), and cL=p E= is the longitudinal wave velocity. Once it was de ned in terms of a dimensionless frequency, the band of analysis does not change when the beam length is varied. Also, the reader can check that this band is rep- resentative for the mechanical system dynamics, once the beam rotates at 2  rad=s , which means that it is excited at 1 Hz. In Figure 10 one can see the distribution of the ex- ural modes as a function of dimensionless frequency, for several values of slenderness ratio. Clearly it is observed that the exural modes are denser in the low frequency range. Further, when the slenderness ratio increases, the modal density in the low frequencies range tend to increase. A completely di erent behavior is observed for the torsional and longitudinal (traction-compression) modes of vibration, as can be seen in Figures 11 and 12, re- spectively. One can note that, with respect to these two modes of vibration, the modal distribution is almost uniform with respect to dimensionless frequency, and invariant to changes in the slenderness ratio. It may also be noted from Figures 10 to 12 that lowest natural frequencies are associated with exural mechanism. This is because beam exural sti ness is much smaller than torsional sti ness, which, in turn, is less than axial sti ness. In other words, it is much easier to bend the beam than twisting it. However, twists the beam is easier than buckling it. The dimensionless frequency band adopted in the analysis corresponds to a maximum dimensional fre- quency offmax= 4cL=L. In this way, a nominal time step oft= (2fmax)1is adopted for time integration. This time step is automatically re ned by the algorithm of integration, whenever necessary, to capture the shock e ects.6.3 Construction of the reduced model In the reduced model construction, are taken into ac- count the mechanical system rigid body modes, as well as the modes of bending, torsion and traction-compression. The construction strategy consists of including: (i) the two rigid body modes (translation and rotation); (ii) all exural modes such that 0 <f5L=cL; (iii) all tor- sional modes such that 0 <f4; (iv) all longitudinal modes such that 0 <f4. In this way, the total number of modes used in FEM model is a function of beam length. In Table 4 the reader can see a comparison, for di erent values of L, of the full FEM model dimension and the correspond- ing reduced order model dimension. Note that using the above strategy, the reduced model dimension is always much smaller than the full model dimension. Table 4 Dimension of the FEM model as a function of beam length. beam length full model reduced model (m) DoFs DoFs 50 306 37 100 3006 49 150 4506 60 6.4 Calculation of static equilibrium con guration Before the beginning of drilling operation, the drill- string is inserted into the borehole, without axial veloc- ity and rotation imposed. Due to gravitational e ects, the column de ects until it reaches a static equilibrium con guration. This con guration can be calculated by temporal integration of the dynamical system de ned by Eqs.(65) and (66), assuming zero initial conditions, i.e., = 0rad=s , andV0= 0m=s. In this way, after a short transient, the system reaches static equilibrium and remains in this con guration inde nitely. An illustration of this equilibrium con guration, for a 100mlong column is presented in Figure 13. In this illustration, one can see the mechanical system sec- tioned by the plane y= 0m, as well as by the planes x=f0;50;100gm. A visual inspection clearly indi- cates that this equilibrium is stable. Moreover, as this equilibrium con guration is the initial state of a real system, it will be used as initial condition in all other simulations reported bellow. An animation which illustrates the beam static equi- librium calculation can be seen in Online Resource 1. Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 17 0 1 2 3 4020406080 dimensionless frequency number of modes flexural modes density (a) slenderness = 312.5 0 1 2 3 4020406080 dimensionless frequency number of modes flexural modes density (b) slenderness = 625.0 0 1 2 3 4020406080 dimensionless frequency number of modes flexural modes density (c) slenderness = 937.5 Fig. 10 Distribution of the exural modes as function of dimensionless frequency, for several values of slenderness ratio. 0 1 2 3 40246810 dimensionless frequency number of modes torsional modes density (a) slenderness = 312.5 0 1 2 3 40246810 dimensionless frequency number of modes torsional modes density (b) slenderness = 625.0 0 1 2 3 40246810 dimensionless frequency number of modes torsional modes density (c) slenderness = 937.5 Fig. 11 Distribution of the torsional modes as function of dimensionless frequency, for several values of slenderness ratio. 0 1 2 3 40246810 dimensionless frequency number of modes longitudinal modes density (a) slenderness = 312.5 0 1 2 3 40246810 dimensionless frequency number of modes longitudinal modes density (b) slenderness = 625.0 0 1 2 3 40246810 dimensionless frequency number of modes longitudinal modes density (c) slenderness = 937.5 Fig. 12 Distribution of the longitudinal modes as function of dimensionless frequency, for several values of slenderness ratio. 6.5 Drill-bit nonlinear dynamic behavior The drill-bit longitudinal displacement and velocity, can be seen in Figure 14. For practical reasons, some scal- ing factors were introduced in the units of measure of these quantities. They allow one to read the displace- ment in \millimeter", and the velocity in \meters per hour". Accordingly, it is noted that, during the interval of analysis, the column presents an advance in forward direction with very small axial oscillations for displace- ment. The axial oscillations for velocity curve are morepronounced, and correspond to the vibration mecha- nism known as bit-bounce , where the drill-bit loses con- tact with the soil and then hits the rock abruptly. This phenomenon, which is widely observed in real systems [60], presents itself discreetly in the case analyzed. Note that the velocity exhibits a mean value of 19 :36 \me- ters per hour", close to the velocity V0= 20 \meters per hour", which is imposed on the beam left end. Also, throughout the \temporal window" analyzed, one can observe packages where the drill-bit velocity presents 18 A. Cunha Jr et al. 0 10 20 30 40 50 60 70 80 90 100−80−5005080beam deflexion in z position (m)deflexion (mm) time = 0.000 sec −80 −50 0 50 80−80−5005080x = 0 m y (mm)z (mm) −80 −50 0 50 80−80−5005080x = 50 m y (mm)z (mm) −80 −50 0 50 80−80−5005080x = 100 m y (mm)z (mm) Fig. 13 Illustration of static equilibrium con guration of a horizontal drillstring with 100 mlength. large uctuations, which can reach up to 40 times the mean value. The drill-bit rotation and angular velocity, can be seen in Figure 15. Now the scale factors allow one to read rotation in \revolution", and angular velocity in \revolution per minute". Thus, what is observed is a almost monotonic rotation. However, when one looks to the angular velocity, it is possible to see packages of uctuations with amplitude variations that can reach up to an order of magnitude. This indicates that drill- bit undergoes a blockage due to torsional friction, and then it is released subtly, so that its velocity is sharply increased, in a stick-slip phenomenon type. This is also seen experimentally [60] in real drilling systems, and a serious consequence of this blockage is the reduction of drilling process eciency. 6.6 Transverse nonlinear dynamics of the beam Observing the beam cross section at x= 50m, for which transversal displacement (top) and velocity (bot- tom) are shown in Figure 16, one can see an asymmetry of the displacement, with respect to the plane z= 0m. This is due to gravity, which favors the beam to move below this plane. Furthermore, one can note that this signal is composed of \packages", which has a recurring oscillatory pattern. As will be seen in section 6.7, these packages present a strong correlation with the number of impacts which the mechanical system is subjected. The evolution of the beam cross-section radial dis- placement, for x= 50m, can be seen in Figure 17, which shows that several transverse impacts occur be- tween drillstring and borehole wall during the drilling 0 2 4 6 8 10−100102030405060 time (s) displacement ( × 10−3 m) drill−bit longitudinal displacement 0 2 4 6 8 10−1000−50005001000 time (s) velocity ( × 1/3600 m/s) drill−bit longitudinal velocityFig. 14 Illustration of drill-bit displacement (top) and drill- bit velocity (bottom). Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 19 0 2 4 6 8 10−20246810 time (s) rotation ( × 2 π rad) drill−bit rotation 0 2 4 6 8 10−600−400−2000200400600 time (s) angular velocity ( × 2 π/60 rad/s) drill−bit angular velocity Fig. 15 Illustration of drill-bit rotation (top) and drill-bit angular velocity (bottom). process. This fact is also reported experimentally [60], and is an important cause of damage for both, well and drillstring. Note that, after an impact, the oscillations ampli- tudes decreases until subtly increase sharply, giving rise to a new impact, and then the entire process repeats again. 6.7 In uence of transverse impacts on the nonlinear dynamics In Figure 18 it is shown the map t2R7!number of shocks 2 N, which associates for any instant tthe number of im- pacts su ered by the mechanical system. The\packages of uctuation"observed in Figures 14 to 16 correspond to transitory periods of the dynamical system, and are highly correlated with the process of collision between beam and borehole wall. This asser- tion can be veri ed if the reader compares the graphs of Figures 14 to 16 with the graph of Figure 18, which 0 2 4 6 8 10−20−15−10−50510 time (s) displacement ( × 10−3 m) transversal displacement in z at x = 50.0 m 0 2 4 6 8 10−4−3−2−101234 time (s) velocity (m/s) transversal velocity in z at x = 50.0 mFig. 16 Illustration of transversal displacement (top) and velocity in z (bottom) when x= 50m. shows the existence of \shock packages". The existence of a correlation is clearly evident. Whenever there is a shock, the system\loses it mem- ory" about previous dynamic behavior, and undergoes a new transient period until reach a steady state again. This behavior is repeated 11 times in the \temporal window" analyzed. Regarding the distribution of impacts along the beam, the mapx2[0;L]7!number of shocks 2N, which associates for any position xthe number of impacts suf- fered by the mechanical system, is shown in Figure 19. It is clear that impacts do not occur near the beam ends. This is natural due to the restrictions of move- ment imposed by the boundary conditions. The impacts between drillstring and borehole wall generate nonlinear elastic deformations in the beam, but without residual deformation e ects. In this contact also occurs energy dissipation, due to the normal shock, and the torsional friction induced by beam rotation. These mechanical contacts also activate exural modes 20 A. Cunha Jr et al. 0 2 4 6 8 105060708090100 time (s) displacement ( × 10−3 m) beam cross−section at x = 50.0 m radial displacement borehole wall Fig. 17 Illustration of beam radial displacement for x= 50m. 0 2 4 6 8100510152025 time (s) number of shocks shock evolution Fig. 18 Illustration of the number of impacts su ered by the mechanical system as function of time. 0 20 40 60 80 1000510152025 position (m) number of shocks shock distribution Fig. 19 Illustration of the number of impacts su ered by the mechanical system as function of position.of vibration associated to high natural frequencies, so that the mechanical system assumes complex spatial con gurations, as can be seen, for several instants, in Figure 20. It is also very clear from Figure 20 that, the me- chanical contacts between beam and borehole wall, do not occur all the time among discrete points, they can also be seen along continuous line segments. For a qualitative illustration of the nonlinear dy- namics, the reader can see Online Resource 2. 6.8 Spectral analysis of the nonlinear dynamics All signals presented above, that are associated with the mechanical system response, have stochastic char- acteristics. Thereby, for a good understanding of them, it is necessary to analyze their spectral content through the power spectral density (PSD) function [37]. The PSDs that are presented in this section (ma- genta line) were estimated using the periodogram method [37], and the smooth curves (blue line) appearing were obtained by a ltering process, using a Savitzky-Golay lter [48]. The PSDs are measured in dB/Hz, where the intensity of reference is adopted as being equal to one. An illustration of PSD functions of drill-bit velocity and angular velocity is show in Figure 21. One can note that, for velocity, the two peaks of highest amplitude correspond to the frequencies 84 :55Hz, and 115:20Hz, respectively. These frequencies are very close to the exural frequencies 83 :65Hz, and 114:27Hz, so that drill-bit axial dynamics is controlled by the transversal mechanisms of vibration. Furthermore, with respect to angular velocity, it is noted a peak standing out in rela- tion to the others. This peak is associated with 7 :92Hz frequency, which is very close to the exural frequency 7:63Hz. In Figure 22 the reader can see an illustration of PSD functions of beam transversal velocity in z and angular velocity around x when x= 50m. The two peaks of highest amplitude, for velocity in z, correspond to the frequencies 143 :20Hz, and 172:50Hz, respec- tively. These frequencies are close to the exural fre- quencies 144 :94Hz, and 174:07Hz, which indicates that lateral vibrations in z, when x= 50m, are in- duced by the transversal vibration mechanism. In ad- dition, in what concerns angular velocity around x, the two peaks of largest amplitude are associated to the fre- quencies 6:93Hz, and 107:10Hz, respectively close to the exural frequencies 6 :59Hz, and 106:18Hz. According to Figure 23, exion is the primary mech- anism of vibration that causes beam and borehole wall impacts, since the highest peak of the PSD shown in this Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 21 0 10 20 30 40 50 60 70 80 90 100−80−5005080beam deflexion in z position (m)deflexion (mm) time = 2.145 sec −80−50 0 5080−80−5005080x = 0 m y (mm)z (mm) −80−50 0 5080−80−5005080x = 50 m y (mm)z (mm) −80−50 0 5080−80−5005080x = 100 m y (mm)z (mm) (a)t= 2:145s 0 10 20 30 40 50 60 70 80 90 100−80−5005080beam deflexion in z position (m)deflexion (mm) time = 4.932 sec −80−50 0 5080−80−5005080x = 0 m y (mm)z (mm) −80−50 0 5080−80−5005080x = 50 m y (mm)z (mm) −80−50 0 5080−80−5005080x = 100 m y (mm)z (mm) (b)t= 4:932s 0 10 20 30 40 50 60 70 80 90 100−80−5005080beam deflexion in z position (m)deflexion (mm) time = 6.214 sec −80−50 0 5080−80−5005080x = 0 m y (mm)z (mm) −80−50 0 5080−80−5005080x = 50 m y (mm)z (mm) −80−50 0 5080−80−5005080x = 100 m y (mm)z (mm) (c)t= 6:214s Fig. 20 Illustration of the mechanical system, for several instants, sectioned by the planes y= 0m, andx=f0;50;100gm. 22 A. Cunha Jr et al. 0 50 100 150 200− 80− 70− 60− 50− 40− 30 frequency (Hz)power spectral density (dB/Hz)drill− bit velocity 84.6 Hz (flexural)115.2 Hz (flexural) 0 50 100 150 200− 15− 10− 50510152025 frequency (Hz)power spectral density (dB/Hz)drill− bit angular velocity (flexural) 7.9 Hz Fig. 21 Illustration of power spectral density functions of drill-bit velocity (top) and angular velocity (bottom). gure is associated with frequency 57 :42Hz, which is close to the exural frequency 57 :77Hz. This result is consistent with the intuition and with what is reported in the literature of drillstrings in vertical con guration [60]. 6.9 Analysis of drilling process eciency The drilling process eciency is de ned as E=Rtf t0Poutdt Rtf t0Pindt; (96) wherePoutis the useful (output) power used in the drilling process, and Pinis the total (input) power in- jected in the system. The output power is due to drill- bit movements of translation and rotation so that Pout= _u+ bit(FBR)++!+ bit(TBR)+; (97) 0 50 100 150 200504540353025 frequency (Hz)power spectral density (dB/Hz)beam velocity in z at x = 50 m 143.2 Hz (flexural)172.5 Hz (flexural) 0 50 100 150 200− 20− 100102030 frequency (Hz)power spectral density (dB/Hz)beamangular velocity at x =50 m (flexural) 6.9 Hz(flexural) 107.1 HzFig. 22 Illustration of power spectral density functions of beam transversal velocity in z (top) and angular velocity around x (bottom) when x= 50m. 0 50 100 150 200403530252015 frequency (Hz)power spectral density (dB/Hz)shock evolution(flexural) 57.4 Hz Fig. 23 Illustration of power spectral density function of number of shocks per unit of time. Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 23 velocity ( × 1/3600 m/s)rotation ( × 2//60 rad/s)efficiency (%) 10 15 20 25 3045505560 2468101214 Fig. 24 Illustration of eciency function contour plot, for an \operating window" de ned by 1 =360m=sV01=120m=s and 3=2rad=s 2 rad=s . The maximum is indicated with a blue cross. where the upper script+means the function positive part. The input power is de ned as Pin= _u(0;t)+(1)++_x(0;t)+(4)+; (98) where the rst and the fourth Lagrange multipliers, re- spectively, represent the drilling force and torque on the beam origin. The reason for considering, in the above de nitions, only functions positive part is that negative powers do not contribute to the drilling process. One can observe the contour map of E, for an \oper- ating window"de ned by 1 =360m=sV01=120m=s and 3=2rad=s 2 rad=s , in Figure 24. Note that, by operating window of a drillstring, one means the subset of R2that provides acceptable values for the pair ( ; V 0). In order to facilitate the results interpre- tation, some scaling factors were introduced in the units of measure. They allow one to read velocity in \meters per hour" and rotation in \rotation per minute". Accordingly, it can be noted in Figure 24 that the optimum operating condition is obtained at the point (V0; ) = (1=144m=s; 5=3rad=s ), which is indicated with a blue cross in the graph. This point corresponds to an eciency of approximately 16%. Suboptimal op- eration conditions occur in the vicinity of this point, and some points near the \operating window" bound- ary show lower eciency. 6.10 Optimization of drillstring rate of penetration In order to optimize the drilling process of an oil well in horizontal con guration, it is necessary to maximize the drillstring ROP into the soil.The instantaneous rate of penetration is given by the function _ ubit(t), de ned for all instants of analy- sis. Meanwhile, only contributes to the column advance, the positive part of this function _ u+ bit(t). In addition, as objective function, it is more convenient to consider a scalar function. Thus, the temporal mean of _ u+ bit(t) is adopted as rate of penetration, and, consequently, ob- jective function of the optimization problem rop( ; V 0) =1 tft0Ztf t=t0_u+ bit(t)dt: (99) Furthermore, respect the material structural limits is indispensable to avoid failures in drillstring during the drilling process. For this reason, von Mises criterion of failure is considered, where it is established that, for all pairs ( ; V 0) in the \operating window", one has UTSmax 0xL t0ttf VM(V0; ; x; t ) 0; (100) where UTSis the material ultimate tensile strength, and VMis the von Mises equivalent stress. Regarding the rate of penetration analysis, \operat- ing window" is de ned by the inequations 1 =360m=s V01=90m=s and 3=2rad=s 7=3rad=s , andUTS= 650106Pa. The contour map of constraint (100), is shown in Figure 25. From the way (100) is written, the Mises criterion is not satis ed when the function is negative, which occurs in a \small neighborhood" of the upper left corner of the rectangle that de nes the \operating window". It is noted that all other points respect the material structural limits. In this way, then, the \oper- ating window" admissible region consists of all points that satisfy the constraint. In Figure 26 the reader can see the contour map ofrop function. Taking into account only points in the admissible region, the maximum of ropoccurs at (V0; ) = (7=720m=s; 2 rad=s ), which is indicated on the graph with a blue cross. This point corresponds to a mean rate of penetration, during the time interval analyzed, approximately equal to 90 \meters per hour". It is worth remembering that the de nition of rop uses temporal mean of _ ubit(t) positive part. In such a way, it is not surprising to nd the maximum value of ropmuch higher than the corresponding velocity, V0 imposed on the column left end. This occurs because, by taking only the function positive part, the rate of penetration value increases. To see how signi cant is the inclusion of _ ubit(t) pos- itive part in the de nition of rop, the reader can see Figure 27. This gure shows the same information as 24 A. Cunha Jr et al. velocity ( × 1/3600 m/s)rotation ( × 2//60 rad/s)UTS − von Mises stress ( × 106 Pa) 10 15 20 25 30 35 40455055606570 0100200300400500 Fig. 25 Illustration of maximum von Mises stress contour plot, for an \operating window" de ned by 1 =360m=sV0 1=90m=sand 3=2rad=s 7=3rad=s . velocity ( × 1/3600 m/s)rotation ( × 2//60 rad/s)rop (× 1/3600 m/s) 10 15 20 25 30 35 40455055606570 20304050607080 Fig. 26 Illustration of rate of penetration function contour plot, for an \operating window" de ned by 1 =360m=sV0 1=90m=s and 3=2rad=s 7=3rad=s . The maximum is indicated with a blue cross. Figure 26, i.e., the contour map of ropfunction, but now considering _ ubit(t) instead of _ u+ bit(t) in the de ni- tion of rop. Note that, in comparison with the contour map of Figure 26, lower values for the function levels are observed, and these values now are closer to V0. Furthermore, the topology of contour lines change, so that no local extreme point can be seen isolated. This example shows the importance of considering _ u+ bit(t) in the de nition of rop. 6.11 Probabilistic analysis of the dynamics For the probabilistic analysis of the dynamic system a parametric approach is used, where the random pa- rameters distributions are constructed according to the velocity ( × 1/3600 m/s)rotation ( × 2//60 rad/s)rop (× 1/3600 m/s) 10 15 20 25 30 35 40455055606570 101520253035Fig. 27 Illustration of the contour plot of the rate of pene- tration function, with an alternative de nition, for an \oper- ating window" de ned by 1 =360m=sV01=90m=s and 3=2rad=s 7=3rad=s . procedure presented in section 5. In this case, the ran- dom variables of interest are characterized by the mean valuesm BR= 400 1=m=s ,mBR= 30103N, and mBR= 0:4, and the dispersion factors  BR= 0:5%, BR= 1%, and BR= 0:5%. To compute the parameters uncertainties propaga- tion through the model, MC method is employed. To analyze the convergence of MC simulations, it is taken into consideration the map ns2N7!convMC(ns)2R, being conv MC(ns) =0 @1 nsnsX n=1Ztf t=t0 q(t;n) 2dt1 A1=2 ;(101) wherensis the number of MC realizations, and kkde- notes the standard Euclidean norm. This metric allows one to evaluate the approximation q(t;n) convergence in mean-square sense. For further details the reader is encouraged to see [55]. The evolution of conv (ns) as function of nscan be seen in Figure 28. Note that for ns= 1024 the metric value has reached a steady value. In this sense, if some- thing is not stated otherwise, all the stochastic simula- tions that follows in this work use ns= 1024. An illustration of the mean value (blue line), and a con dence band (grey shadow), wherein a realization of the stochastic dynamic system has 95% of probability of being contained, for drill-bit longitudinal displace- ment and velocity is shown in Figure 29. For sake of reference, the deterministic model, which the numer- ical results were presented earlier, is also shown and called nominal model (red line). It is observed that for Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 25 0200 400 600 800100050.9550.95250.95450.95650.95850.96 number of MC realizations convergence metric study of Monte Carlo convergence Fig. 28 This gure illustrates the convergence metric of MC simulation as function of the number of realizations. the displacement, mean value and nominal model are very similar. Meanwhile, for the velocity, mean value presents oscillations that are correlated with the nomi- nal model, but with very di erent amplitudes. Regard- ing the con dence band, there are signi cant ampli- tudes in the instants that corresponds to the uctua- tion packages, and negligible amplitudes in the other moments. Fixing time in t= 10s, it is possible to analyze the behavior of drill-bit longitudinal velocity through its normalized PDF, which is presented in Figure 30. In this context normalized means a distribution of proba- bility with zero mean and unit standard deviation. It is observed an unimodal behavior, with maximum value occurring in a neighborhood of the mean value, with small dispersion around this position. In Figure 31, the reader can see the nominal model, the mean value, and the 95% probability envelope of drill-bit rotation and angular velocity. A good agree- ment between rotation nominal model and mean value is observed, and the con dence band around it is neg- ligible. On the other hand, with respect to the angular velocity, it is possible to see discrepancies in the nominal model and mean value amplitudes. These di erences oc- cur in the instants when the system is subject to shocks, as in the case of drill-bit longitudinal velocity. The band of uncertainty shows that the dispersion around mean value increases with time, due to accumulation of un- certainties, but also in reason of the impacts, once its amplitude increases a lot near the instants where the mean value presents large uctuations, i.e., the instants which are correlated with the beam and borehole wall impacts. Fort= 10s, the reader can see the normalized PDF of drill-bit angular velocity in Figure 32. It is noted 0 2 4 6 810−100102030405060 time (s) displacement ( × 10−3 m) drill−bit longitudinal displacement nominal mean value 95% prob. 0 2 4 6 810−1000−500050010001500 time (s) velocity ( × 1/3600 m/s) drill−bit longitudinal velocity nominal mean value 95% prob.Fig. 29 Illustration of the nominal model (red line), the mean value (blue line), and the 95% probability envelope (grey shadow) for the drill-bit longitudinal displacement (top) and velocity (bottom). −40 −30 −20 −10 0 1000.20.40.60.811.21.4 velocity (normalized) probability density function drill−bit longitudinal velocity PDF Fig. 30 Illustration of the normalized probability density function of drill-bit longitudinal velocity. 26 A. Cunha Jr et al. 0 2 4 6 810−2024681012 time (s) rotation ( × 2π rad/s) drill−bit rotation nominal mean value 95% prob. 0 2 4 6 810−600−400−2000200400600800 time (s) angular velocity ( × 2π/60 rad/s) drill−bit angular velocity nominal mean value 95% prob. Fig. 31 Illustration of the nominal model (red line), the mean value (blue line), and the 95% probability envelope (grey shadow) for the drill-bit rotation (top) and angular ve- locity (bottom). again an unimodal behavior, with maximum again near the mean value. But now a large dispersion around the mean can be seen. Moreover, in Figure 33 it is shown the nominal model, the mean value, and the 95% probability envelope of beam transversal displacement and velocity in z at x= 50m. Here the mean values of both, velocity and dis- placement, present correlation with the nominal mod- els. Indeed, both present discrepancies in the oscillation amplitudes, especially the velocity, that are more pro- nounced, as before, in the instants wherein the system is subject to impacts. The con dence bands present mean- ingful amplitudes, what evidentiates a certain level of dispersion around the means, which are more signi - cant, as expected, at the instants of impact. The PDF of drilling process eciency function is shown in Figure 34. One can observe a unimodal distri- bution with the maximum around 16% and wide disper- −10 −5 0 5 1000.10.20.30.40.5 angular velocity (normalized) probability density function drill−bit angular velocity PDFFig. 32 Illustration of the normalized probability density function of drill-bit angular velocity. 0 2 4 6 810−20−15−10−50510 time (s) displacement ( × 10−3 m) transversal displacement in z at x = 50.0 m nominal mean value 95% prob. 0 2 4 6 810−4−3−2−101234 time (s) velocity (m/s) transversal velocity in z at x = 50.0 m nominal mean value 95% prob. Fig. 33 Illustration of the nominal model (red line), the mean value (blue line), and the 95% probability envelope (grey shadow) for the beam transversal displacement (top) and velocity in z (bottom) at x= 50m. Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 27 0 20 40 60 80 100012345 efficiency (%) probability density function efficiency PDF Fig. 34 Illustration of the probability density function of the drilling process eciency. 0 100 200 300 4000100200300400500600 rate of penetration ( × 1/3600 m/s) probability density function rate of penetration PDF Fig. 35 Illustration of the probability density function of the rate of penetration function. sion between 0 and 40%, declining rapidly to negligible values outside this range. Finally, in Figure 35 one can see the PDF of drill- string rate of penetration function. One notes an uni- modal behavior in a narrow range between 20 and 50 \meters per hour", with the maximum around 30 \me- ters per hour". 6.12 Robust optimization of drillstring rate of penetration To improve the level of con dence of drilling process op- timization, uncertainties intrinsic to the problem should be taken into account. This leads to a robust optimiza- tion problem, i.e, optimization under uncertainty where the range of random parameters are known, but not necessarily their distribution [4, 51, 7, 6, 8, 58, 3].Taking into account the uncertainties, through the parametric approach presented in section 5, drill-bit ve- locity becomes the stochastic process Ubit(t;), so that the random rate of penetration is de ned by ROP(V0; ; ) =1 tft0Ztf t=t0_U+ bit(t;)dt: (102) In the robust optimization problem, who plays the role of objective function is the expected value of the random variable ROP(V0; ; ), i.e.,E ROP(V0; ; ) . Regarding the restriction imposed by the von Mises criteria, now the equivalent stress is the random eld VM(V0; ; x; t; ), so that the inequality is written as UTSmax 0xL t0ttf VM(V0; ; x; t; ) 0: (103) However, the robust optimization problem consid- ers as restriction a probability of the event de ned by inequality (103), P8 >< >:UTSmax 0xL t0ttf VM(V0; ; x; t; ) 09 >= >;1Prisk, (104) where 0< Prisk<1 is the risk percentage acceptable to the problem. A robust optimization problem very similar to this one, in the context of a vertical drillstring dynamics, is considered in [42]. To solve this robust optimization problem it is em- ployed a trial strategy which discretizes the \operating window" in a structured grid of points and then eval- uates the objective function E ROP(V0; ; ) and the probabilistic constraint (104) in these points. Accordingly, it is considered the same \operating window" used in the deterministic optimization prob- lem solved above, i.e., 1 =360m=sV01=90m=s and 3=2rad=s 7=3rad=s , in addition to UTS= 650106PaandPrisk= 10%. Each MC simu- lation in this case used 128 realizations to compute the propagation of uncertainties. Concerning the simulation results, the probabilis- tic constraint (104) is respected in all grid points that discretize the \operating window". Thus, the admissi- ble region of robust optimization problem is equal to the \operating window". In what follows, the contour map of function E ROP(V0; ; ) can be see in Fig- ure 36. Note that the maximum, which is indicated on the graph with a blue cross, occurs at the point 28 A. Cunha Jr et al. (V0; ) = (1=90m=s; 7=3rad=s ). This point is located in the boundary of the admissible region, in the upper right corner, and corresponds to a expected value of the mean rate of penetration, during the time interval analyzed, approximately equal to 58 \meters per hour". velocity ( × 1/3600 m/s)rotation ( × 2//60 rad/s)mean ROP ( × 1/3600 m/s) 10 15 20 25 30 35 40455055606570 2025303540455055 Fig. 36 Illustration of the contour plot of the mean rate of penetration function, for an \operating window" de ned by 1=360m=sV01=90m=s and 3=2rad=s  7=3rad=s . The maximum is indicated with a blue cross in the upper right corner. This result says that, in the\operating window"con- sidered here, increasing drillstring rotational and trans- lational velocities provides the most robust strategy to maximize its ROP into the soil. This is in some ways an intuitive result, but is at odds with the result of the deterministic optimization problem, which provides an- other strategy to achieve optimum operating condition. 7 Concluding remarks A model was developed in this work to describe the non- linear dynamics of horizontal drillstrings. The model uses a beam theory, with e ects of rotatory inertia and shear deformation, which is capable of reproducing large displacements that the beam undergoes. This model also considers the friction and shock e ects due to transver- sal impacts, as well as, the force and torque induced by bit-rock interaction. It should be emphasized that the paper not only proposes a mechanical analysis of this complex dynamic system, but, through this analysis, also provides a for- mulation for this type of system, and a methodology, that can be reused to other rotating slender mechanical systems subjected to friction and shock e ects. Numerical simulations showed that the mechanical system of interest has a very rich nonlinear dynam-ics, which reproduces complex phenomena such as bit- bounce, stick-slip, and transverse impacts. The study also indicated that large velocity uctuations observed in the phenomena of bit-bounce and stick-slip are corre- lated with the transverse impacts, i.e., with the number of shocks per unit time which the system is subjected. Also, the mechanical impacts cause the beam to as- sume complex spatial con gurations, which are formed by exural modes associated to high natural frequen- cies. A study aiming to maximize drilling process e- ciency, varying drillstring velocities of translation and rotation was presented. The optimization strategy used a trial approach to seek for a local maximum, which was located within \operating window" and corresponds to an eciency of approximately 16%. The probabilistic analysis of the nonlinear dynam- ics showed that, with respect to the velocities, nominal model and mean value of the stochastic model di er sig- ni cantly. Furthermore, at the instants which the sys- tem was subjected to mechanical impacts, it was pos- sible to see a more pronounced dispersion around the mean value. Regarding the velocities probability distri- butions, it was noticed a unimodal behavior essentially. Two optimizations problems, one deterministic and one robust, where the objective was to maximize drill- string rate of penetration into the soil, respecting its structural limits, were formulated and solved. The solu- tions of these problems provided two di erent strategies to optimize the ROP. Finally, it sounds stressing that the mathematical model used in this work has not gone through any pro- cess of experimental validation. This is because exper- imental data for this type of system is dicult to be obtained, and to construct an experimental apparatus in real scale is virtually impossible. An interesting pro- posal for future work would be the construction of an experimental test rig, in reduced scale, that emulates the main aspects of a real drillstring. The model used in this study could be validated, following, for instance, the methodology presented in [1], with the aid of ex- perimental measurements taken from this reduced ap- paratus. A Geometric nonlinearly force coecients This appendix presents the coecients which appears in ge- ometric nonlinearity force of Eq.(49). For the sake of sav- ing space, in the following lines it is used the abbreviations: Sx= sinx, andCx= cosx. Computational modeling of the nonlinear stochastic dynamics of horizontal drillstrings 29 1=EI41 +u0 v00 y+w00 z Sx0 x+ (105) EI41 +u0 v00 zw00 y Cx0 x+ ksGA1 +u0  zv0yw0 Sx ksGA1 +u0  yv0+zw0 Cx; 2=ksGI4 y 02 y+02 z 0 x0 z + (106) ksGA w0+u0y2 +u0  ksGA1 +u0  v0Sxw0Cx ; 3=ksGI4 z 02 y+02 z +0 x0 y + (107) ksGA v0+u0z2 +u0  ksGA1 +u0  w0Sx+v0Cx ; 4=EA1 21 +u0 v02+w02 +1 2u023 +u0 + (108) EI4 Sx v00 zw00 y Cx v00 y+w00 z 0 x+ EI41 +u0 02 x+3 2 02 y+02 z + ksGA Cx yw0zv0 Sx yv0+zw0 + ksGA1 +u0 2 y+2 z ; 5=EA u0+1 2 u02+v02+w02 v0+ (109) EI4 202 x+1 2 02 y+02 z v0+ EI41 +u0 0 zSx0 yCx 0 x+ ksGA1 +u0  zySxzCx ; 6=EA u0+1 2 u02+v02+w02 w0+ (110) EI4 202 x+1 2 02 y+02 z w0+ EI41 +u0 0 ySx0 zCx 0 x+ ksGA1 +u0  y+yCxzSx ; 7=EI4 u02+ 2 u0+v02+w02 0 x+ (111) EI41 +u0 v00 zw00 y Sx EI41 +u0 v00 y+w00 z Cx+ EI6 402 x+ 2 02 y+02 z 0 x+ ksGA z0 yy0 z ;8=EI4 3u0+1 2 3u02+v02+w02 0 y+ (112) EI41 +u0  w0Sxv0Cx 0 x+ EI6 202 x+3 2 02 y+02 z 0 y+ ksGI4 z0 x+0 y 2 y+2 z ; and 9=EI4 3u0+1 2 3u02+v02+w02 0 z+ (113) EI41 +u0  v0Sxw0Cx 0 x+ EI6 202 x+3 2 02 y+02 z 0 z+ ksGI4 y0 x+0 z 2 y+2 z : Acknowledgements The authors are indebted to Brazilian agencies CNPq, CAPES, and FAPERJ, and French agency COFECUB for the nancial support given to this research. The rst author is grateful for the institutional support re- ceived from PUC-Rio and Universit e Paris-Est to carry out this work.