. This paper presents a model to describe the nonlinear dynamics of a drillstring in horizontal conﬁguration, which is intended to correctly predict the three-dimensional dynamics of this complex structure. This model uses a beam theory, with eﬀe cts of rotatory inertia and shear deformation, which is capable of reproducing the large displac ements that the beam undergoes. Also, it considers the eﬀects of torsional friction and normal shock due to the transversal impacts between the rotating beam and the borehole wall, as w ell as, the force and the torque induced by the bit-rock interaction. This is done as a ﬁr st eﬀort to solve a robust optimization problem, which seeks to maximize the rate of penetration of the drillstring into the soil, to reduce the drilling process costs. Numerical simulat ions reported in this work shown that the developed computational model is able to quantitatively well describe the dynamical behavior of a horizontal drillstring, once its reproduces some phenom ena observed in real drilling systems, such as bit-bounce, stick-slip, and transverse impacts . Mecánica Computacional Vol XXXIII, págs. 1517-1527 (artículo completo) Graciela Bertolino, Mariano Cantero, Mario Storti y Federico Teruel (Eds.) San Carlos de Bariloche, 23-26 Setiembre 2014 Copyright © 2014 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar 1 INTRODUCTION A drillstring is a device, used to drill oil wells, which pres ents an extremely com- plex three-dimensional nonlinear dynamics. The dynamical system associated with this physical system involves the nonlinear coupling between th ree diﬀerent mechanisms of vi- bration (longitudinal, transverse, and torsional), as wel l as lateral and frontal shocks, due to drill-pipes/borehole and drill-bit/soil and impacts re spectively ( Spanos et al. ,2003). Traditionally, a drillstring conﬁguration is vertical, bu t directional or even horizontal con- ﬁgurations, where the boreholes are drilled following a non -vertical way, are also possible. Once oil drilling a topic of great relevance in the context of engineering, the dynamics of a vertical drillstring has been studied in several works ( Chevallier ,2000;Ritto et al. , 2009,2010;Chatjigeorgiou ,2013;Liu et al. ,2013;Depouhon and Detournay ,2014). How- ever, although of most of the oil wells today be drilled with c olumns using non-vertical conﬁgurations, very few papers in the open literature model s drillstring in directional conﬁgurations ( Sahebkar et al. ,2011;Hu et al. ,2012;Ritto et al. ,2013). Aiming to ﬁll the gap in the scientiﬁc literature on horizonta l drillstring dynamics, this work presents the modeling of a drillstring in a horizon tal conﬁguration. This model takesintoaccountthethree-dimensionaldynamicsofthest ructure, aswellasthetransver- sal/torsional eﬀects of shock, which the structure is subje ct due to the impact with the borehole wall. Also, the model considers the bit-rock intera ction eﬀects, and the weight of the drilling ﬂuid. This rest of this paper is organized as follows. The mathemat ical modeling of the nonlinear dynamics appears in section 2. Then, in section 3, the results of numerical simulations are presented and discussed. Finally, in secti on4, the main conclusions are emphasized and some directions for future work outlined. 2 MATHEMATICAL MODELING 2.1 Mechanical system of interest The mechanical system of interest in this work is sketched in Figure1. It consists of a horizontal rigid pipe (illustrated as the pair of stationa ry rigid walls), perpendicular to gravity acceleration g, which contains in its interior a deformable tube under rotat ion (rotating beam), subjected to three-dimensional displace ments. This deformable tube has a length L, cross section area A, and is made of a material with mass density ρ, elastic modulus E, and Poisson ratio ν. It loses energy through a mechanism of viscous dissipation, proportional to the mass operator, with dampi ng coeﬃcient c. Inside the tube there is a ﬂuid without viscosity, with mass density ρf. Concerning the boundary conditions, the rotating beam is blocked for transversal di splacements in both extremes; blocked to transversal rotations on the left extreme; and, o n the left extreme, has a constant angular velocity around x equal to Ω, and an imposed longitudinal velocity V0. 2.2 Beam theory The beam theory adopted takes into account the rotatory iner tia and shear deforma- tion of the beam cross section. Also, as the beam is conﬁned wit hin the borehole, it is reasonable to assume that it is undergoing small rotations i n the transverse directions.A. CUNHA JR, C. SOIZE, R. SAMPAIO 1518 Copyright © 2014 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar x Lyz Rext Rint Figure 1: Schematic representation of the rotating beam which models t he horizontal drillstring. By another hand, large displacements are observed in x,y, andz. Therefore, the analysis that follows uses a beam theory which assumes large rotation in x, and large displacements the three spatial directions, which couples the longitudinal, transverse and torsional vibrations ( Bonet and Wood ,2008). Regarding the kinematic hypothesis adopted for the beam the ory, it is assumed that the three-dimensional displacement of a beam point, occupy ing the position ( x,y,z) at the instant of time t, can be written as ux(x,y,z,t) =u−yθz+zθy, (1) uy(x,y,z,t) =v+y(cosθx−1)−zsinθx, uz(x,y,z,t) =w+z(cosθx−1)+ysinθx, where letters u,v, andware used to denote the displacements of a beam neutral ﬁber point inx,y, andzdirections, respectively, while θx,θy, andθzrepresent rotations of the beam around the x,y, andzaxes respectively. Note that these quantities depend on the positionxand the time t. 2.3 Friction and shock eﬀects This rotating beam is also able to generate normal shocks and torsional friction in random areas of the rigid tube, which are respectively descr ibed by the Hunt and Cross- ley shock model Hunt and Crossley (1975), and the standard Coulomb friction model. Therefore, the force of normal shock is given by FFS=−kFS1δFS−kFS2δ3 FS−cFS|δ|3˙δFS, (2) and the Coulomb frictional torque by TFS=−µFSFFSRbhsgn/parenleftBig ˙θx/parenrightBig . (3) In the above equations, kFS1,kFS2andcFSare constants of the shock model, while µFSis a friction coeﬃcient, Rbhis the borehole radius, and sgn(·) the sign function. The ˙ is an abbreviation for time derivative, and the paramet erδFS=r−gap, where r=√ v2+w2, is dubbed indentation , and is a measure of penetration in the wall of a beam cross section, such as illustrated in Figure 2.Mecánica Computacional Vol XXXIII, págs. 1517-1527 (2014) 1519 Copyright © 2014 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar gap rgap r r−gap≤0 r−gap>0 Figure 2: Illustration of the indentation parameter in a situation withou t impact (left) or with impact (right). 2.4 Bit-rock interaction eﬀects At the right extreme of the rotating beam act a force and a torq ue, which emulate the eﬀects of interaction between the drill-bit and the soil. Th ey are respectively given by FBR=  ΓBR/parenleftBig exp/parenleftbig −αBR˙u(L,·)/parenrightbig −1/parenrightBig ,for ˙u(L,·)>0 (4) 0, for ˙u(L,·)≤0 and TBR=−µBRFBRξBR/parenleftBig ˙θx/parenrightBig , (5) where Γ BRis the bit-rock limit force; αBRis the rate of change of bit-rock force; µBR bit-rock friction coeﬃcient; and ξBRis a regularization function, which takes into account the dimension of length, to the Eq.( 5) gives a torque. The expression for the bit-rock interaction models above were, respectively, proposed by Ritto et al. (2013) andKhulief et al.(2007). 2.5 Variational formulation of the nonlinear dynamics Using a modiﬁed version of the extended Hamilton’s principle, t o include the eﬀects of dissipation, one can write the weak form of the nonlinear e quation of motion of the mechanical system as M/parenleftBig ψ,¨U/parenrightBig +C/parenleftBig ψ,˙U/parenrightBig +K(ψ,U) =FNL/parenleftBig ψ,U,˙U,¨U/parenrightBig , (6) whereMrepresents the mass operator, Cis the damping operator, Kis the stiﬀness operator, and FNLis the nonlinear force operator. Also, the ﬁeld variables and the ir weight functions are lumped in the vectors ﬁelds U=/parenleftbig u,v,w,θ x,θy,θz/parenrightbig , andψ=/parenleftbig ψu,ψv,ψw,ψθx,ψθy,ψθz/parenrightbig .A. CUNHA JR, C. SOIZE, R. SAMPAIO 1520 Copyright © 2014 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar The above operators are respectively deﬁned by M/parenleftBig ψ,¨U/parenrightBig =/integraldisplayL x=0ρA(ψu¨u+ψv¨v+ψw¨w)dx+ (7) /integraldisplayL x=0ρfAf(ψv¨v+ψw¨w)dx+ /integraldisplayL x=0ρI4/parenleftBig 2ψθx¨θx+ψθy¨θy+ψθz¨θz/parenrightBig dx, C/parenleftBig ψ,˙U/parenrightBig =/integraldisplayL x=0cρA(ψu˙u+ψv˙v+ψw˙w)dx+ (8) /integraldisplayL x=0cρI4/parenleftBig 2ψθx˙θx+ψθy˙θy+ψθz˙θz/parenrightBig dx, K(ψ,U) =/integraldisplayL x=0EAψ′ uu′dx+ (9) /integraldisplayL x=0EI4/parenleftBig ψ′ θyθ′ y+ψ′ θzθ′ z/parenrightBig dx+ /integraldisplayL x=02κsGI4ψ′ θxθ′ xdx+ /integraldisplayL x=0κsGA/parenleftBig/parenleftbig ψθy+ψ′ w/parenrightbig/parenleftbig θy+w′/parenrightbig +/parenleftbig ψθz−ψ′ v/parenrightbig/parenleftbig θz−v′/parenrightbig/parenrightBig dx, and FNL/parenleftBig ψ,U,˙U,¨U/parenrightBig =FKE/parenleftBig ψ,U,˙U,¨U/parenrightBig +FSE(ψ,U) + (10) FFS(ψ,U)+FBR/parenleftBig ψ,˙U/parenrightBig +FG(ψ), where FKE=−/integraldisplayL x=02ρI4ψθx/parenleftBig θy¨θz+˙θy˙θz/parenrightBig dx (11) +/integraldisplayL x=02ρI4ψθy/parenleftBig θy˙θ2 z+˙θx˙θz/parenrightBig dx −/integraldisplayL x=02ρI4ψθz/parenleftBig θy¨θx+θ2 y¨θz/parenrightBig dx −/integraldisplayL x=02ρI4ψθz/parenleftBig ˙θx˙θy+2θy˙θy˙θz/parenrightBig dx is a nonlinear force due to inertial eﬀects;Mecánica Computacional Vol XXXIII, págs. 1517-1527 (2014) 1521 Copyright © 2014 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar FSE=/integraldisplayL x=0/parenleftbig ψθxΓ1+ψθyΓ2+ψθzΓ3/parenrightbig dx+ (12) /integraldisplayL x=0/parenleftbig ψ′ uΓ4+ψ′ vΓ5+ψ′ wΓ6/parenrightbig dx+ /integraldisplayL x=0/parenleftBig ψ′ θxΓ7+ψ′ θyΓ8+ψ′ θzΓ9/parenrightBig dx, is a nonlinear force due to geometric nonlinearity; FFS=Nnodes/summationdisplay m=1/parenleftbig FFS(vψv+wψw)/r+TFSψθx/parenrightbig/vextendsingle/vextendsingle/vextendsingle x=xm, (13) is a nonlinear force due to the eﬀects of friction and shock; FBR=FBRψu/vextendsingle/vextendsingle/vextendsingle x=L+TBRψθx/vextendsingle/vextendsingle/vextendsingle x=L, (14) is a nonlinear force due to the bit-rock interaction; and FG=−/integraldisplayL x=0/parenleftbig ρA+ρfAf/parenrightbig gψwdx, (15) is a linear force due to the gravity. The nonlinear functions Γn, withn= 1,···,9, in Eq.(12) are very complex and, for sake of space limitation, are not p resented here. See Cunha Jr (2015) for details. The weak form of the initial conditions reads M/parenleftbig ψ,U(0)/parenrightbig =M(ψ,U0), (16) and M/parenleftBig ψ,˙U(0)/parenrightBig =M/parenleftBig ψ,˙U0/parenrightBig , (17) whereU0and˙U0, respectively, denote the initial displacement, and the ini tial velocity ﬁelds. The model presented above is an adaptation, for the case of ho rizontal drillstrings, with some variations in the friction and shock treatment, of the m odel proposed by Ritto et al. (2009) to describe the nonlinear dynamics of vertical drillstrin gs. 2.6 Discretization of the model equations The Eqs.( 6), (16) and (17) are discretized by means of the standard ﬁnite element method ( Hughes,2000), using an interdependent interpolation scheme ( Reddy,1997), which adopts aﬃne functions for the axial displacement/tor sional rotation, and Hermite cubic polynomials for the transverse displacements/rotat ions.A. CUNHA JR, C. SOIZE, R. SAMPAIO 1522 Copyright © 2014 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar Therefore, one arrives in the following initial value probl em [M]¨Q(t)+[C]˙Q(t)+[K]Q(t) =F/parenleftBig Q(t),˙Q(t),¨Q(t)/parenrightBig , (18) and [M]Q(0) =Q0,and [M]˙Q(0) =˙Q0, (19) whereQ(t) is the nodal displacement vector (translations and rotati ons),˙Q(t) is the nodal velocity vector, ¨Q(t) is the nodal acceleration vector, [ M] is the mass matrix, [ C] is the damping matrix, [ K] is the stiﬀness matrix, and Fis a nonlinear force vector, which contains contributions of an inertial force and a force of ge ometric stiﬀness. The geometric boundary conditions are included as constrai nts, via the method of Lagrange multipliers. Nominally, they are the velocity of tra nslation,V0, and the velocity of rotation, Ω, which are imposed at the left end of the beam. 2.7 Reduction of the nonlinear dynamics To reduce the computational cost of the simulations, the ini tial value problem of Eqs.(18) and (19) is projected in a vector space of dimension Nred, spanned by the linear modes associatedto the conservative part of the underlying linear dynamical system. This results in the reduced initial value problem given by [M]¨q(t)+[C]˙q(t)+[K]q(t) =f/parenleftbig q(t),˙q(t),¨q(t)/parenrightbig , (20) and q(0) =q0,and ˙q(0) =˙q0, (21) which is integrated using the Newmark method ( Newmark ,1959), and the nonlinear system of algebraic equations, resulting from the time disc retization, is solved by a ﬁxed point iteration. 3 RESULTS AND DISCUSSION In order to simulate the nonlinear dynamics of the mechanica l system, the physical parameters presented in the Table 1are adopted, as well as the length L= 35m, the ro- tationalandaxialvelocitiesinx, respectivelygivenbyΩ = 2πrad/s, andV0= 1/720m/s. For the geometry discretization, 105 ﬁnite elements are use d. This results in FEM model with 636 degrees of freedom. In the reduced order model, 51 DO F are considered. The dynamics is investigated for a“temporal window”of 90 s, with a nominal time step ∆t= 69ms, which is reﬁned whenever necessary to capture the eﬀects of shock. For the initial conditions, the static equilibrium conﬁguration o f the beam is adopted. The drill-bit longitudinal displacement and velocity, can be seen in Figure 3. It is noted that, during the interval of analysis, the column pres ents an advance in the forward directionwithsmallaxialoscillations. Theseaxialoscil lations, whicharemorepronounced in the velocity curve, correspond to the vibration mechanis m known as bit-bounce , where the drill-bit loses contact with the soil and then hits the ro ck abruptly. This phenomenon is widely observed in real systems ( Spanos et al. ,2003).Mecánica Computacional Vol XXXIII, págs. 1517-1527 (2014) 1523 Copyright © 2014 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar Table 1: Physical parameters of the mechanical system that are used in th e simulation. parameter value unit ρ 7900kg/m3 ρf 1200kg/m3 E 203GPa ν 0.3 — Rint 40mm Rext 60mm Rbh 70mm c 0.03 — g 9.81m/s2 kFS11×1010N/m kFS21×1016N/m3 cFS 1×109(N/m3)/(m/s) µFS 1×10−5— ΓBR 250kN αBR 180 1/(m/s) µBR400×10−4— 0 20 40 60 80−20020406080100120140 time (s) displacement (mm) drill−bit longitudinal displacement 0 20 40 60 80−3−2−101234 time (s) velocity (m/s) drill−bit longitudinal velocity Figure 3: Illustration of the drill-bit displacement (left) and of th e drill-bit velocity (right). The drill-bit rotation and angular velocity, can be seen in F igure4. What it is observed now is a almost monotonic rotation. However, when one looks to the angular velocity, it is possible to see packages of ﬂuctuations with amplitude va riations that can reach up to four orders of magnitude. This indicates that the drill-bit undergoes a blockage due to the torsional friction, and then it is released subtly, so that i ts velocity is sharply increased, in astick-slip phenomenon type. This is also seen experimentally ( Spanos et al. ,2003) in real drilling systems. The evolution of the radial displacement, for x= 20, of the beam cross-section can be seen in the Figure 5. Analyzing this ﬁgure it is clear that transverse impacts betwe en the drillstring and the borehole wall occur during the drilling process, which is also reported experimentally ( Spanos et al. ,2003).A. CUNHA JR, C. SOIZE, R. SAMPAIO 1524 Copyright © 2014 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar 0 20 40 60 80−20020406080100 time (s) rotation (rev) drill−bit rotation 0 20 40 60 80−6−4−20246x 104 time (s) angular velocity (rpm) drill−bit angular velocity Figure 4: Illustration of the drill-bit rotation (left) and of the drill -bit angular velocity (right). 0 20 40 60 8050556065707580 time (s) displacement (mm) beam cross−section at x = 20.0 m radial displacement borehole wall Figure 5: Illustration of the beam radial displacement for x= 20m. 4 CONCLUDING REMARKS A model was developed in this work to describe the nonlinear d ynamics of horizontal drillstrings. The model uses a beam theory, with eﬀects of ro tatory inertia and shear deformation, which is capable of reproducing the large disp lacements that the beam un- dergoes. This model also considers the eﬀects of friction an d shock due to the transversal impacts between the beam and the borehole wall, as well as, th e force and the torque induced by the bit-rock interaction. Numerical simulations reported in this work shown that the dev eloped computational model is able to quantitatively well describe the dynamical behavior of a horizontal drill- string, once its reproduces some phenomena observed in real drilling systems, such as bit-bounce, stick-slip, and transverse impacts.Mecánica Computacional Vol XXXIII, págs. 1517-1527 (2014) 1525 Copyright © 2014 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar In a future work, the authors intend to develop a stochastic m odeling of the nonlinear dynamics of horizontal drillstrings, in order to quantify t he uncertainties associated with this problem, which are due to the variability of its paramet ers (Schu¨eller,2007), and/or epistemic in nature, i.e., result of the ignorance about the physics of the problem ( Soize, 2013). Also, in a next step, they want to solve an robust optimizati on problem, which seeks to maximize the rate of penetration of the column into t he soil (Ritto et al. ,2010). ACKNOWLEDGMENTS The authors are indebted to the Brazilian agencies CNPq, CAPES , and FAPERJ, and theFrenchagencyCOFECUBfortheﬁnancialsupportgiventoth isresearch. Theywould also like to acknowledge professors Anas Batou and Thiago Ritt o for valuable discussions on this work.