Academic Journal of Science and Technology ISSN: 2771-3032 | Vol. 13, No. 3, 2024 281 Nonlinear Dynamic Response Characteristics of Deep‐ water Riserless Drilling‐string System Under Combined Excitation of Platform Motion and Current Yi Zhang1, Jun Liu1, 2, *, Shuang Liang1, Lin Chen1, Yili Chen1 1School of Mechatronic Engineering, Southwest Petroleum University, Chengdu, 610500, China 2School of Mechanical Engineering, Chengdu University, Chengdu, 610106, China *Corresponding author, Jun Liu Abstract: The study develops a model for axial-lateral-torsional (ALT) bidirectional coupling nonlinear vibrations in riserless drill-strings, using energy and finite element methods. The model incorporates environmental factors such as platform motion, sea currents, and wave loads, as well as operational interactions including wellbore contact and the drill-rock interface. Validated with field test results and data from four wells in the South China Sea, the findings emphasize the importance of heave compensation devices in offsetting the effects of sea motions on drilling operations. The analyses indicate that the weight on bit is mainly impacted by heave movements; however, slow drift motions can exacerbate the nonlinearity of bit rotation. Managing the amplitude and frequency of heave motions is crucial to prevent resonances that can arise between the drill-string and the bit, ensuring improved safety and efficiency in deep-water drilling processes. Keywords: Deepwater drilling; Riserless Drilling; Drill-string; Coupling vibration. 1. Introduction The drilling technology with riserless mud recovery(RMR, Fig.1) is promising and disruptive deepwater drill technology (Yang et al., 2019). It controls wellbore pressure through a subsea pump, creating a dual-gradient pressure effect that expands the safe density window between formation pore pressure and fracture pressure (Thakar et al., 2018). At the same time, it eliminates the conventional riser, significantly reduces the platform load, reduces the casing level, increase the casing running depth, shortens the well construction cycle and extends the water depth limit of drilling operations, which can solve many problems in deepwater drilling technology (Xu et al., 2011). Figure 1. Schematic diagram of deepwater riserless drilling 282 The drill-string system is a core piece of equipment for the riserless drilling process. Its stability not only affects the drilling efficiency, but also concerns the safety of the entire drilling process. Existing studies have shown that harmful vibrations, such as Stick-slip vibration, drill jump, and eddy motion, are the main causes of premature failure of drill-string and bits, deviation of borehole trajectory, reduction of drilling speed, and reduction of measurement while drilling accuracy (Guzek et al., 2015; Kessai et al., 2022). Therefore, systematic, comprehensive and accurate dynamic analysis of drill-string system can better understand and master the working condition of drill-string, which is of great significance to improve drilling efficiency and reduce drilling cost. Depending on the mode of vibration, the drill-string vibrations can be divided into axial, lateral, and torsional vibrations. Each individual vibration mode has a significant effect on the drill-string (Ren et al., 2018). Limited by cognitive level and computational conditions, early studies focused on the analysis of uncoupled unidirectional vibrations in vertical wells. In spite of the limitations of the method and the means, it is important to reveal to some extent the law of vibration in each direction of the drill-string (Ghasemloonia et al., 2015). With the development and application of vibration-while-drilling instruments and the revelation of drill-string vibration mechanism, it is proved that drill-string vibration is a whole, and axial, lateral, and torsional vibration exist and interact with each other at the same time. Any single aspect of vibration research is incomplete (Li et al., 2017). With the development of computer technology, some scholars began to focus on the axial-lateral-torsional (ALT) coupled vibrations of the drill-string. Dykstra (1996) was one of the earliest scholars to study the ALT coupled vibration of drill- string. He used the energy method to develop a coupled vibrational model of the drill-string. In his paper, he gave a detailed modeling process, which provided a reference for later scholars. Baumgar (2000) derived a ATL nonlinear differential equation of drill-string to predict the effects of mud on bit pressure, drill-string buckling, and bit stick-slip vibration. Tucker and Wang (1999) established the ALT coupled vibration model of drill-string considering the nonlinear strain and realized the coupled vibration of three vibration forms.This model was used to discuss the stability of axisymmetric drill-string configurations in vertical boreholes under both coupled ATL perturbations as well as general non-perturbative coupled vibrational states under extreme conditions of lateral whirl. de Moraes and Savi (2019) established a four-degree-of-freedom mass ALT coupled vibration model of drill-string, and discussed drilling skip, stick-slip vibration, eddy motion and their combined effects. Some scholars (Khoshroo and Eftekhari, 2022; Liu and Gao, 2017; Rajabali et al., 2020) studied the vibration behavior of local horizontal drill-string by using ALT coupled vibration model, providing theoretical support for the combination design of horizontal drill-string, especially BHA tools. Some scholars (Cai et al., 2022; Liu et al., 2022; Li et al., 2020) have begun to use whole-hole ALT coupled vibration model to study drill-string vibration. However, these studies mainly focus on drill-string dynamics in traditional land drilling, while there are relatively few studies on drill-string dynamics in offshore drilling (Qin et al., 2022). Compared with conventional land drilling drill-string, the biggest difference of RMR drill-string is that its top is connected to an offshore drilling platform or drilling ship, and the drill-string is directly exposed to seawater, resulting in a more complex dynamic environment (Su et al., 2013). Under the action of wind and waves, the floating drilling platform will produce three types of motion, including heaving, swaying and slow drift, of which the heaving and slow drift have obvious influence on drilling operations. In the process of drilling, the platform heave motion drives the drill-string and downhole drill tool up and down, which makes the pressure of the bottom hole bit fluctuate and even the bit skip, which not only affects the drilling efficiency, but also may lead to fatigue fracture of the drill-string and premature failure of the bit. At the same time, the changing axial tension of the drill-string will affect the lateral stiffness, and the lateral displacement of the drill-string will fluctuate constantly under the action of seawater load. The lateral slow drift motion of the drilling platform will cause a large change in the lateral displacement of the drill-string in the seawater segment, which not only changes the bending angle and bending moment of the drill-string in the seawater segment, but also fluctuates the axial displacement of the drill-string under the action of the axial and lateral coupling (Xing and Liang, 2015), which further influences the pressure of the bit. The conventional theory and understanding of drill-string dynamics cannot meet the requirements of analysis and design, thus special research on RMR drill-string dynamics must be carried out. However, there is a lack of research on riserless drill-stringing. Wada et al, (2018) analyzed the drilling data obtained by the scientific drilling ship Chikyu with a 8000m drill-string, and found that the axial vibration frequency of the drill pipe is higher than the frequency of the heave motion. Xia et al, (2019) investigated various loads in riserless drilling and pointed out that the variable load was significantly reduced after the drilling platform was equipped with riserless drilling technology, up to more than 50%. Su et al, (2013) established a lateral vibration model of seawater section riserless drill-string under the combined excitation of wave current and platform motion, and analyzed the effects of platform drift, drill-string size, and drill-string length on the dynamic stress in the drill-string. Robello (2013) established a nonlinear dynamic model of drill-string in deepwater highly deviated well and investigated the influence of mud depth, platform drift, borehole inclination angle and azimuth angle on the drilling string vibration response and pointed out that large vibration of drilling string can be avoided by changing excitation parameters. In view of the lack of studies on drill-string dynamics in deep-water riserless drilling, in this paper, we derive ALT bidirectionally coupled nonlinear drill-string elements suitable for RMR drill-string using energy method and finite element method, and propose a solution suitable for the bidirectionally coupled drill-string elements. The ALT coupled nonlinear vibration model of the RMR drill-string in deep water is established by considering the drilling platform heave and slow drift motion, wind current, sea current and wave load, wall contact and drill-bit rock interaction. We analyze the impact of the platform heave and slow drift motion on the drilling operation and the vibration of the drill- string at different heave and slow drift motion amplitudes and frequencies. 283 2. ALT Coupled Nonlinear Vibration Model of Drill-string The drill-string continuous system is discretized by finite element method. As shown in Fig.2, the drill-string is discrete into n drill-string units, where the platform (drill ship) is connected to the upper end of the drill-string by a heave compensation device, and the drill bit is connected to the lower node of the nth drill-string unit. Then the equation for the drill-string dynamics is:  MD+CD+KD= F (1) Where D , D , D and F are the generalized acceleration array, generalized velocity array, generalized displacement array, and load array of the whole drill-string system, respectively. M , C K and are the mass matrix, damping matrix, and stiffness matrix, respectively, of the whole drill-string system. Element 1 Element 2 Element n-1 Element n Drill string Bit · · · Tk Tc Tm Platform Figure 2. Drill-string discrete element system In this paper, we use two nodal twelve degrees of freedom drill-string elements. As shown in Fig. 3, each node contains one axial displacement w , two lateral displacements u and v , one torsional angle z , and two transverse rotation angles x and y . In total, a drill-string unit contains two nodes with 12 degrees of freedom. 1 x j z u1 v1 w1 v2 θx1 θy1 θz1 u2 w2 θx2 θy2 θz2 y 2 Figure 3. Drill-string element with two nodes and 12 degrees of freedom The potential energy 1 2 12( , ,..., )V d d d of drill-string element is a function of the generalized displacement of the element, and V is assumed to be zero at the equilibrium position of the drill-string element. Expanding V in Taylor series and retaining the second-order small quantities of the generalized shift, one obtains: 12 12 1 1 1 2 ij i j i j V k d d     (2) where ijk is the stiffness coefficient and is expressed as follows: 2 ( , 1, 2,...,12)ij i j V k i j d d      (3) Similarly, the kinetic energy 1 2 12( , ,..., )T d d d , expressed by the generalized displacement of the drill-string element, and the mass coefficient ijm can be obtained. 12 12 1 1 2 1 2 ( , 1, 2,...,12) ij i j i j ij i j T m d d T m i j d d                 (4) By introducing the drill-string element generalized displacement array ( )idd , the stiffness matrix ( )ijkk , and the mass matrix ( )ijmm , the potential energy V and kinetic energy T of drill-string element can be expressed as: T T 1 2 1 2 V T       d kd d md  (5) Therefore, the element stiffness matrix and the element mass matrix can be obtained by obtaining the expressions of potential energy and kinetic energy of the drill-string element. Fig. 4 (a) Three-dimensional Cartesian coordinate system was established. The center of the upper end of the drill-string is the origin of the coordinate, and the axis is oriented vertically downward, consistent with the direction of gravity. The structure and forces of the RMR drill-string are extremely complex. To facilitate the derivation of the kinetic and potential energy equations for the drill-string element, the following assumptions are made: (1) The cross section of the drill-string is a ring, and each position can have any geometric size and material properties, but each segment remains constant. (2) The effect of the drill-string joints is neglected. (3) The upper end of the drill-string rotates at a constant speed. (4) The contact between drill-string and borehole wall is regarded as elastic collision. 284 O1 x y  ,r （ ） ,r  （ ） y z O1 u v x w O2 O2 (a) (b) Figure 4. (a) 3D coordinates of drill-string system; (b) Sketch of drill-string cross section Figure 4 (b) is the schematic diagram of the torsion deformation of the drill-string, and 2O is the section centroid. 1O is the initial position of the drill-string axis; ,z z t（ ） is the torsional angular displacement of the cross section, rad. Then the relative torsional angular displacement of any profile is z t   (6) Where,  is the rotation speed of the platform turntable, rad/s. 2.1. Kinetic energy The displacements of the drill-string elements include translational and torsional displacements. As shown in Fig. 4 (a), ( , )u z t and ( , )v z t represent the transverse displacements along the x and y axes, respectively, while the axial displacements along the z axis are represented by ( , )w z t . Thus, the translational displacement of the drill- string element can be obtained. cr uI vJ wK      (7) Equation (7) is applied to obtain the first partial derivative of time to obtain the translational velocity of the element: c cv r uI vJ wK           (8) Thus, the translational kinetic energy ktT of the drill- string element can be expressed as follows:  2 2 2 2 0 0 1 1 d d 2 2 e el l kt cT Av z A u v w z         (9) where,  is the density of the drill-string, A is the cross-sectional area of the drill-string, el is the length of the drill-string element. The rotational kinetic energy of the drill-string element is divided into three parts, namely the rotational kinetic energy around the x , y and z axes. The drill-string has a circular cross-section, so the moment of inertia with respect to the x and y axis is the same, but the moment of inertia with respect to the z axis is different. The drill-string element rotational velocity is denoted as: x y zI J K         (10) The rotation angle of the drill-string around the x and y axes is related to the transverse displacements u and v as follows. tan , tanx y u v z z        (11) On the other hand, in the drill-string, which has a extremely high slenderness ratio, the angle is small, so: tan , tanx x y y u v z z            (12) The angular velocities x and y of the drill-string element around the x and y axes can be obtained by computing the first-order partial derivatives of the angular velocity with respect to time. Similarly, the angular velocity z of the drill-string element around the z axis can be obtained by computing the first order partial derivative of the twist angle ( , )z z t with respect to time , ,x y z zu v          (13) where,    ' z    , similarly,     2 2 '' z    . By substituting equation (13) into equation (10), the rotational speed of drill-string element can be obtained: zu I v J K          (14) Thus, the rotational kinetic energy of the drill-string element can be expressed as  2 2 2 0 0 0 0 1 1 0 0 d d 2 2 0 0 e e x l l kr x y z y z z I T I z Iu Iv I z I                                              (15) where, I is the moment of inertia of drill-string section, I is the polar moment of inertia of the section Therefore, the kinetic energy of the drill-string element is    2 2 2 2 2 2 0 0 1 1 d d 2 2 e el l kt kr zT T T A u v w z Iu Iv I z                 (16) 2.2. Potential energy In order to derive the potential energy of the drill-string element, the nonlinear geometric equations for the drill-string element must first be given. In the derivation, the drill-string motion is treated as the combined motion of translation and rotation. As shown in Fig. 4 (b), the displacement of any point 285 in the cross section is equal to the displacement of the section centroid 2O plus the displacement of the torsional deformation  , so the displacements of any point  ,r  on the section along the x and y axes are:     cos cos cos sin sin sin cos sin x y u u r r u x x y u v r r v y x y                                     (17) When  is small, cos 1, sin    , plugging into the above equation to sort: x yu u y u v x    ， (18) Then the displacement ( xu , yu , zu ) is the displacement field function of the coordinate system ( x , y , z ), whose expression is as follows: , ,x y z x y u u y u v x u w xu yu w xu yv                (19) Following Green's strain formula, the strain in the displacement field can be expressed as (Yin, 2019):  1 1 , , , , 2 2 i j k k ij u u u u i j k x y z i i i j                    (20) where  is the nonlinearity or coupling factor and =1 is the actual nonlinear geometric equation of the element. When =0 , equation (20) degenerates into a linear geometric equation. For drill-string with ultra-long and fine ratio, the stress of the drill-string can be assumed to be dominated by axial stress. 0,x y xy     0,xy x y z       (21) where,  is the Poisson's ratio. The strain can be obtained according to Einstein's summation rule. 22 2 2 2 2 1 2 1 ( ) ( ) ( ) 2 yxz z z z z uuu u z z z z w xu yv u y v x w xu yv                                                  (22a)   1 1 2 2 1 1 = ( ) ( ) 2 2 y yx x xz z z xz z z z u uu u uu u u z x z x z x z x y v x w xu yv u                                          (22b)   1 1 2 2 1 1 ( ) ( ) 2 2 y y yx xz z z yz z z z u u uu uu u u z y z y z y z y x u y w xu yv v                                             (22c) It can be obtained from the element potential energy formula. 2 2 0 1 d ( )d d 2 el A V E G A z         (23) where, E is the Modulus of elasticity, G is the Modulus of shear,  is the Shear strain,    222 = 2 2xz yz   . It has also been argued that the following relations exist in the integration process:   3 3 2 2 2 2 2 2 d 0 d 0 d 0 d 0 d 0 d 0 d 0 d d d A A A A A A A A A A x A y A x A y A xy A xy A x y A x A I y A I x y A I                                    (24) In combination with Equations (22)-(24), the potential energy V expression of the drill-string element is obtained by ignoring the higher-order minor terms:       2 2 2 3 2 0 0 2 2 2 2 2 0 d 2 2 d 2 2 2 2 2 3 3 d 2 e e e l l z z z l z z z EA GI V w w u w v w z u v u v z EI u v w u v u v w u w v z                                                           (25) 2.3. Drill-string element shape function In the above derivation of the expressions of drill-string kinetic energy and potential energy, the continuous displacement function is used, while the expressions of drill- string element kinetic energy and potential energy in Equation (5) are expressed by the generalized displacement of drill- string element nodes. Therefore, the continuous displacements in Equations (16) and (25) need to be replaced by nodal displacements, which can be done by introducing a shape function. As can be seen from the derivation, the drill- string displacement consists of two lateral displacements u and v , one axial displacement w , two lateral rotation angles x and y , and one torsional angle z . According to Equation (12), the two lateral angles are not independent degrees of freedom and can be obtained from the lateral displacement to the axial partial derivative. Thus, the cubic Hermite interpolation function is used to express the transverse displacement of the drill-string in the x and y directions, and the Lagrangian interpolation function is used to express the axial and torsional displacements: T T 1 2 T T 3 4 , , z u v w        φ d φ d φ d φ d (26) Where: 1 1 2 2T 1 1 1 1 2 2 2 2 2 3 2 3 2 3 3 2 T 1 2 3 2 2 3 2 2 3 2 3 T 2 2 3 2 3 2 2 3 21 , ,0,0,0,0, , ,0,0,0,0 3 2 20,0,1 , ,0,0,0 z z e e e e e e e e e e e e u v u vu v w u v w z z z z z z z z z z z zz l l l l l l l l z z z zz l l l l                            d φ φ 2 3 3 2 2 3 2 T 3 T 4 3 2,0, , ,0,0 0,0,0,0,1 ,0,0,0,0,0, ,0 0,0,0,0,0,1 ,0,0,0,0,0, e e e e e e e e z z z z l l l l z z l l z z l l                           φ φ (27) 286 The shape functions of the drill-string element: 1φ , 2φ , 3φ and 4φ are obtained above, which can be put into equations (16) and (25) to obtain the kinetic and potential energy of the drill-string element expressed by the generalized displacement of the drill-string element node. 2.4. Dynamic equation for the drill-string element Let  iff be the conservative generalized force corresponding to the generalized displacement of the drill- string element and L T V  be the Lagrangian function. General form of the Lagrangian of the second kind: d ( 1, 2, ...,12) d i ii L L f i t dd          (28) Combined with Equations (2) and (4), the dynamics equation of the drill-string element can be obtained: 12 1 ( ) ( 1,2,...,12)ij j ij j i j m d k d f i      (29) Write matrix form as: e ie m d k d f (30) The mass matrix em and stiffness matrix iek of the drill-string element: 1 2 3 4 T T T T 1 2 3 4 4 5 5 6 6 7 7 8 e ie i i i i i i                m m m m m k k k k k k k k k k k k k+ + (31) where i indicates different coupling vibration modes. When i=1, the transverse and torsional vibration of the drill-strings will cause axial coupled vibration. i=2 denotes that the axial tension affects the lateral and torsional stiffness. Some scholars (Xing and Liang, 2015) have studied the case of i=1, but more scholars (Xie et al., 2019) are concerned about the case of i=2. The linear mass and stiffness matrix of drill-string element is as follows:       1 0 2 0 3 0 4 0 1 0 2 0 3 0 d d d 2 d d d 2 d e e e e e e e l T T l T T l T l T l T T l T l T A I z A I z A z I z EI z EA z GI z                                          1 1 1 1 2 2 2 2 3 3 4 4 1 1 2 2 3 3 4 4 m φ φ φ φ m φ φ φ φ m φ φ m φ φ k φ φ φ φ k φ φ k φ φ (32a) where, 1m and 2m are the lateral mass matrices, 3m is the axial mass matrix, 4m is the torsional mass matrix, The element mass matrix is the same as that in the uncoupled case, because the coupling term is not considered in the kinetic energy function of Equation (16). And 1k is the lateral linear stiffness matrix, 2k is the axial linear stiffness matrix, 3k is the torsional linear matrix. The nonlinear stiffness matrix of the drill-string element is expressed as follows:             T T T 41 30 T T 42 2 1 0 0 T T T 51 0 T T 52 2 1 0 0 T T 61 0 62 2 1 d 2 d d 2 3 d 2 3 d d 2 d e e e e e e e l l l l l l l EA z EA w w z z l EI z EI w w z z l EI z EI w w l                                           1 1 2 2 1 1 2 2 3 1 1 2 2 1 1 2 2 3 4 4 4 k φ d φ φ φ φ k φ φ φ φ k φ d φ φ φ φ k φ φ φ φ k φ d φ φ k φ     T 0 T 7 3 32 1 0 T T T 8 0 d d ( ) d e e e l l l z EA zw w l E G I z                                    4 4 1 2 1 2 φ k φ φ k φ d φ φ φ φ (32b) where, 4ik and 5ik are the axial and lateral coupling stiffness matrices, 6ik is the axial and torsional coupling stiffness matrix, 7k is the nonlinear axial stiffness matrix, 8k is the lateral and torsional coupling stiffness matrix. The dynamic equation (30) of drill-string element considering viscous damping of drilling fluid is rewritten as: e ie ie  m d c d k d f  (33) where, iec is the element damping matrix, The viscous damping of drilling fluid is commonly treated as Rayleigh damping, which is represented by a linear combination of stiffness matrix and mass matrix: ie e ie  c m k (34) where,  is the mass proportional coefficient,  is the stiffness proportional coefficient, which are calculated as follows: 0 0 2 2 = , =i j i j i j          (35) where, 0 is the target damping ratio. In this paper, 0 =0.03 is selected. i is the first-order natural frequency. To avoid excessive damping of the high-frequency vibrations, 287 the natural frequency of order 10 (Moharrami et al., 2021). For linear systems, the natural frequencies of the drill- string are different in axial, lateral, and torsional ways, so the damping matrices of axial, lateral, and torsional ways are also different, which can be expressed as: 1 2 1 3 2 4 3 u u u w w w                   c m m k c m k c m k （ ） (36a) For the coupling part, the damping matrix is expressed as:       T T 4 4 5 5 T 6 6 7 T 8 8 uwi u w i i i i w i w i i ww w w u u                           c k k k k c k k c k c k k (36b) Thus, the damping matrix of the element can be expressed as: ie u w uwi w i ww u        c c c c c c c c (37) It can be noted that the element stiffness matrix and damping matrix are related to the element generalized displacement d at every moment, so they are time-varying. 2.5. Solution scheme The drill-string element is assembled to obtain the dynamic equation (1) of the overall discrete form of the drill-string. Since the stiffness matrix and damping matrix of the drill- string element have two different forms, equation (1) is rewritten as: ( ) ( ) ( )i it t t  MD C D K D F  (38) Where, C and K are the overall mass matrix and damping matrix of the drill-string, respectively. Equation (38) is a set of second-order ordinary differential equations, which can be solved by standard methods for ordinary differential equations. However, in the finite element dynamic response analysis, when the order of the matrix is high, the efficiency of the standard algorithm is low, so Newmark-β is adopted in this paper to solve the discrete equations (38). Due to the large slenderness ratio of the drill-string, the influence of the axial tension of the drill-string on the lateral and torsional stiffness can not be ignored under the action of gravity. For the traditional land drill-string, due to the existence of borehole wall constraints, there will be no large lateral displacement, so the impact of lateral displacement on the axial direction can be ignored. However, for the riserless drill-string, the current load and the platform movement force the drill-string to produce a large lateral displacement in the seawater section, and the axial vibration caused by the lateral motion can not be ignored. Whether case 1i  or case 2i  , it can only satisfy the analysis of unidirectional coupled vibration. For a deep-sea riserless drill-string, the influence of axial tension on transverse/torsional stiffness and the axial vibration caused by transverse displacement should not be ignored, therefore, it needs a two-way coupled vibration analysis. In this paper, a two-way coupling scheme is used to solve the equations (22), and the two one-way coupling systems are solved the independently and in turn to achieve the purpose of two-way coupling vibration analysis. The solution idea is given as follows: Ⅰ: Substituting the displacement  1it D at the previous time into the unidirectional coupling matrices 1K , 1C , 2K and 2C at the current time. Ⅱ: Calculating the displacement field at the current time 1( )itD , 1( )itD , 1( )itD , 2 ( )itD , 2 ( )itD and 2 ( )itD . Ⅲ: The coupled vibration displacement ( )itD , ( )itD , ( )itD of the system at the current time is obtained from the lateral and torsional displacement of 1( )itD , 1( )itD , 1( )itD and the axial displacement of 2 ( )itD , 2 ( )itD , 2 ( )itD . 2.6. Boundary condition (1) Upper end of drill-string The upper end of the riserless drill0string is connected to the drill platform by a heave compensator, and moves horizontally with the latter. Therefore, the corresponding boundary conditions can be described as: (0, ) , (0, ) (0, ) cos( ), (0, ) 0 (0, ) sin( ), (0, ) 0 z TF t F t u t S u t v t S v t                (39) where TF is the top tension from the heave compensator and its detailed description is can be found in appendix A1, is the horizontal displacement of the drilling platform, is the angle between horizontal displacement orientation and the x- axis, is the rotation speed of the rotary table. Considering that the slow drift motion of the drilling ship mainly consists of two parts: the average offset of the drilling platform and the long-term slow drift motion simulating the response of the positioning system of the drilling ship, the following formula is used to describe the horizontal displacement of the drilling platform (Wang et al., 2015). 0 2 ( ) sinL L L S t S S t T           (40) where 0S is the average offset of the drilling ship, LS is the amplitude of slow drift motion, LT is the period of slow drift motion, and L is the phase difference. The heaving motion of a drilling platform is formed by the action of ocean waves and may be approximated by a simple harmonic motion, whose amplitude is typically considerably 288 smaller than the wave amplitude and depends on the wave height and the structure size of the platform, with a period comparable to that of the wave. Here, it is assumed that the heaving motion of the drilling platform is similar to the simple harmonic motion. ( ) sin( )h sw t h t (41) where sw is the heave motion of the drilling platform, sh is the amplitude,  is wave circular frequency. (2) Mudline Wellhead The lateral displacement of the drill-string at the mudline wellhead is constrained, so the boundary conditions here can be written as: ( , ) 0, ( , ) 0 ( , ) 0, ( , ) 0 u H t v H t u H t v H t       (42) where H is the water depth at mudline wellhead. (3) Ocean loads The ocean load on the drill-string is described by Morison equation (Morison et al., 1950):  2 1( , ) ( ) 4 2 H H M W b s D W H h b s h b s Df z t C a a C D v v v v v v        (43) where hv is the horizontal velocity of the water point caused by the sea current, bv is the horizontal velocity of the water point caused by the wave, ba is the horizontal acceleration of the water point caused by the wave, sv is the horizontal velocity of the drill-string, sa is the horizontal acceleration of the drill-string. hv , bv and ba are determined in Appendix A2. (4) Drill-string wellbore contact in the formation When the radial displacement of the drill-string reaches the annular space between the drill-string and the wellbore, the drill-string will be constrained by the wellbore, and the two will contact each other. In this paper, Hertz contact theory is used to describe this boundary problem (Christoforou and Yigit, 1997; Zhao et al., 2016). As shown in Figure , the radial contact force rF , the tangential friction force tF and the friction torque cT can be expressed as: y r cT tF rF x  Figure 5. Drill-string-wellbore contact diagram   3 2 2 2 0 2 sign 2 c H c H r h c H t b r p c t d D d D F K r r d D r F F d T F                          (44) where 2 (=1.9 )h p p Eb b K D D is the collision rigidity coefficient, r is the radial displacement of the drill-string, = 2 p pD d b  is the wall thickness of the drill-string, cd is the inner diameter of the borehole wall and  is the relative rotation speed of the borehole wall of the drill-string. (5) Bit-rock interaction at bottom hole The boundary at the lower end of the drill-string, i.e. the bit-rock interaction, is described by the friction torque received by the drill bit (Moharrami et al., 2021): 0 0 0 0 , (stick mode), sgn( ), (transition from stick to slip), T sgn( ), (slip mode). ( ) st b d st st b d st b sl b b d s b b sl b b b T if V and T T T T if V and T T T if V T R W T R W                           (45) where dV is critical rotational speed, bW is weight on the drill bit, ( )b  is the exponential decay friction coefficient of the drill. For the determination of )b （ and bW , refer to Appendix A3. 2.7. Establishment of initial working conditions In actual drilling work, the riserless drill-string is subjected to the tension of the top hook, the pressure of the bottom bit and the gravity in the axial direction at the same time before the turntable starts to rotate, and the seawater section is subjected to the transverse current load, so there is a certain initial displacement in the axial and transverse direction. Taking initial drillship displacement, top tension, bit pressure and current load and ignoring inertial force and damping force of equation (38), the balance equation of drill-string under initial static force is obtained as follows: 0 0i t t K D F (46) The displacement field of the drill-string at the initial time can be obtained by solving equation (46). It is worth noting that from the expression of equation (32b), it can be seen that the stiffness matrix iK changes with the change of the numerical results of displacement 0tD , and the stiffness 289 matrix is changed in each solution result so that the equation does not hold. In this paper, Equation (46) is solved by several iterations, the calculation result is brought into the update system stiffness matrix iK , and the displacement 0tD is obtained by solving again. When the maximum error of the two displacement values is less than the set error precision, the iteration is stopped to obtain the initial system displacement. The error accuracy is taken to be 1e-4. The process of riserless drill-string displacement field initialization and bidirectional coupling solution is shown in Fig. 6. Input parameters and division of drilling units Calculate initial load array initial displacement array Calculate mass matrix 0tF 0tD Calculate stiffness matrix 1 2K K， 1 0 0 i t t K D F 2 0 0 i t t K D F Obtained the lateral and torsional displacement Obtained the axial displacement 0 i tD 1 4 0 0 1ei i t t     D D1i i  No Yes M Calculate stiffness matrix Calculate damping matrix Calculate load array 1 2K K， 1 2C C， 1 1  MD C D K D F  F 2 2  MD C D K D F  Obtained the lateral and torsional acceleration velocity and displacement Obtained the axial acceleration velocity and displacement ( ) ( ) ( )it it itD D D 、 、 1it it  Newmark-β Initial displacement  field calculation Bidirectional coupled  solution Figure 6. The process of riserless drill-string displacement field initialization and bidirectional coupling solution Based on the above models and MATLAB software, we developed the computational program: coupling vibration of riserless drill-string for ALT coupling nonlinear vibration analysis of deepwater riserless drill-string. Next, we will test the effectiveness and adaptability of the model and code. 2.8. Model validation At present, the research on the dynamics of riserless drill- string is still limited, and the experimental and field test study on vibration have not been reported. Therefore, this paper uses Kyllingstad and Nessjøen (2009) and Teng et al, (2017)'s drill-string field test results, as well as four example well parameters (LD-A, LD-B, LD-C and LD-D) in Ledong, South China Sea to analyze the correctness and validity of the model. (1) Correctness verification Tab.1 and Tab.2 show the drill-string field parameters of Kyllingstad and Nessjøen (2009) and Teng et al, (2017) respectively. Fig.7 shows comparisons between the numerical simulations of the model and the field test results. The red curve in Fig.7(a) is from the data of Kyllingstad field test within 4200~4240 s. It can be seen that the maximum angular velocity of the bit obtained from the model in this paper and the two field tests are consistent in both trend and magnitude. In Fig.7(a), the numerical simulation and field test results show that the maximum rotational speed of the bit is 21.96 rad/s and 21.67 rad/s, respectively, and the frequency of stick- slip vibration of the bit is 0.2 Hz and 0.19Hz, respectively, and the maximum speed error is only 1%. In Fig. 7(b), the maximum rotational speed of field test data is 34.54 rad/s, which is obviously higher than the peak speed in other measured data. Therefore, the maximum rotational speed of 34.51rad/s is accidental and can not reflect the general law of stick-slip vibration of the bit. The maximum speed should be removed in the data comparison, and the maximum speed of 31.68 rad/s is obtained. Then, by comparing the drilling bit rotational speed amplitude, the accuracy of model is 88.96%. In addition, the stick-slip vibration frequencies are 0.102 Hz and 0.108 Hz, respectively. The comparison of the numerical results of the proposed model with two field tests proves the correctness of the proposed model. Table 1. Model calculation parameters Item Value Item Value Inner/Outer diameter of the drill-string 1 (m) 0.1214/0.1397 Inner diameter of the HBA (m) 0.0635 Length of the drill-string 1 (m) 6200 Outer diameter of the HBA (m) 0.2032 Inner/Outer diameter of the drill-string 2 (m) 0.1086/0.127 Length of the HBA (m) 100 Length of the drill-string 2 (m) 800 Rotational velocity (rad/s) 9.42 290 Table 2. Model calculation parameters Item Value Item Value Inner/Outer diameter of the drill-string (m) 0.1086/0.127 Inner diameter of the HBA (m) 0.08573 Length of the drill-string (m) 6100 Outer diameter of the HBA (m) 0.2032 Rotational velocity (rad/s) 12.56 Length of the HBA (m) 100 Figure 7. Comparison of rotational speed between the proposed model and the field test (2) Validity verification The stick-slip vibration of the drill-string is a deleterious movement in the drilling process that seriously affects the effective drilling of the bit. That is, the longer the viscous phase in the stick-slip vibration, the smaller the ROP. Then, the ROP can be increased by improving the stick-slip vibration of the bit (Vromen et al., 2019). Tab.3 show the parameters of four wells(LD-A, LD-B, LD-C, LD-D) in Ledong block, South China Sea. The Viscous time and slippage time of the drill-string were calculated(as show in red and blue, respectively, in Fig.8), which were compared with the measured ROP data(as show in black curve of Fig.8). Compared with the field measured ROP data, it is found that the variation trend of viscous time is negatively correlated with the variation trend of ROP, while the slippage time is positively correlated, indicating that with the increase of viscous time of the bit, the time of holding drill increases, which obstructs the rock breaking efficiency and reduces the ROP. This shows that the proposed model can effectively predict the ROP variation trend of the drill-string by calculating the stick-time, which proves the validity of the model. Table 3. Model calculation parameters Item Value Item Value Well Depth (m) 4000/4050/ 4430/4260 Length of heavy weight drill pipe (HWDP) (m) 132 Inner diameter of the drill-string (m) 0.1186 Inner diameter of HWDP (m) 0.0825 Outer diameter of the drill-string (m) 0.1397 Outer diameter of HWDP (m) 0.1397 Length of the HBA (m) 88.49~159.28 WOB (kN) 30~180 Inner diameter of the HBA (m) 0.08573/0.0572 Ground rotary speed (rad/s) 5.23~13.39 Outer diameter of the HBA (m) 0.2032/0.1651 Simulation time (s) 40 3. Dynamic Analysis Tab. 4 shows the geometric model and material parameters of the drill-strings. Tab. 5 shows the load parameters in the seawater section, where the current velocity and wind speed are set to 0.93m/s and 15.43m/s respectively under the normal drilling conditions in the South China Sea (Mao et al., 2016). Based on field experience (Wang and Liang, 2013), the resistance coefficient is set to 1.2 from the sea surface to 150m depth, and set to 0.7 from 150m depth to the mudline, and the inertial force coefficient is set to 2.0 in the whole seawater section. Tab. 6 shows the parameters of wellbore contact load and bit-rock interaction. In order to observe the effect of the lateral movement on the drill-string response, it is assumed that the slow drift motion of the platform and the direction of the current velocity are both in the x -axis direction. 291 Figure 8. The relation curve between stick-slip vibration time and mechanical speed of example well Table 4. Parameters of drill-string system Property Value Unit Outer diameter of the drill-string 0.1270 m Inner diameter of the drill-string 0.1086 m Length of the drill-string 2500 m Outer diameter of the drill collar 0.2286 m inner diameter of the drill collar 0.0762 m Length of the drill collar 150 m Density of material 7850 kg/m3 Young's modulus of material 210 GPa Shear modulus of material 79.6 GPa Bit radius 0.1555 m Rotary speed 5 rad/s Piston cross-sectional area 0.24 m2 Initial pressure 3.25 MPa Total volume of compressed air 10.5 m3 Specific heat ratio of air. 1.2 Zero dimension Table 5. Load parameters of seawater section Property Value Unit Depth of seawater 1200 m Density of seawater 1025 kg/m3 Wave height 6 m Wave length 10.2 m Wave period 11.2 s Maximum measured wind speed 15.43 m/s Maximum surface current speed 0.93 m/s Friction depth 200 m Wind speed factor 0.03 Zero dimension Drag coefficient CD 1.2(between water surface and 150m) 0.7(between 150m and seabed) Zero dimension Inertia cofficient CM 2.0(at all depths) Zero dimension 292 Table 6. Load parameters of stratigraphic section Property Value Unit Static friction coefficient 0.8 Zero dimension Kinetic (Coulomb) friction coefficient 0.5 Zero dimension Decay coefficient 0.9 Zero dimension Threshold velocity 0.000001 m/s drill-string wall friction coefficient 0.3 Zero dimension Bit pressure 80 kN 3.1. Comparison of heave compensation with and without heave compensation The slow drift motion of the platform has an average drift of 15 m, an amplitude of 10 m, and a period of 80 s. Due to the initialization of the static method, the drill-string was in a static state in the seawater section in the initial state, so the phase angle of the slow drift motion of the platform is set to / 2 , making the lateral velocity of the platform in the initial state also zero. The plateau heave has an amplitude of 0.3m and a period of 16 s. Fig. 9 shows a comparison of the time-domain response of the drill-string with and without the heave compensation device. Fig. 9 (a) shows the axial displacement curve of the top of the drill-string, which is consistent with the heave motion of the platform without the heave compensation device. When the heave compensation device is present, the lateral slow drift motion of the platform plays a major role. Under the lateral drag of the platform, the axial displacement of the drill-string first increases and then decreases, and the influence of the heave motion is effectively weakened. Fig. 9 (b) shows the response curve of WOB. It can be found that WOB fluctuates violently and even bit jumping (WOB=0) occurs when there is no heave compensation device. It can also be seen that the heave motion has a large impact on the WOB, while the WOB with heave compensation device is relatively stable. Fig. 9 (c) shows the lateral displacement response at the midpoint of the drill-string in the seawater section (depth 600 m), which is significantly different in two cases. In addition to the lateral movement of the platform, the lateral response of the drill-string without heave compensation device fluctuates obviously under the influence of heave motion. This is due to the fact that in the absence of heave compensation, the axial tension of the drill-string fluctuates dramatically, resulting in a change in the lateral stiffness, which fluctuates in response to the seawater load. Fig. 9 (d) shows the curve of bit speed. It can be found that without heave compensation device, the bit speed fluctuates more violently, and its maximum rotation speed even reaches 6 times of the platform turntable speed, and obvious bit inversion occurs, which is easy to cause damage to the bit. Therefore, for riserless drilling, heave compensation devices can not only compensate the platform heave motion, reduce bit pressure fluctuations, avoid bit jumping and reversal phenomena, but also reduce the lateral vibration of the drill- string in seawater. Figure 9. Time-history responses with and without heave compensation: (a) axial displacement at the top of the drill-string, (b) WOB, (c) lateral displacement at the midpoint of the drill-string in the seawater section, and (d) bit speed 293 3.2. Consider platform heave motion and drift motion The mean drift of the platform's slow drift motion is 20 m, the amplitude is 15 m, the period is 80 s, the phase angle is / 2 . The amplitude of the platform's heave motion is 1.5 m and the period is the same as the wave period of 11.2 s. Case 1 considers only slow drift motion, Case 2 considers only heave motion, and Case 3 considers both slow drift and heave motion. Fig. 10 shows the time-history response diagram of the drill-string under three conditions. Fig. 10 (a) and (b) show the axial displacement at the top of the drill-string and the bit pressure, respectively. It can be seen that the response results of Case 3 are approximately the superposition of the response results of Case 1 and Case 2, in which the heave movement has a significant influence on bit weight, while the slow drift movement mainly has an influence on the axial displacement at the top of the drill column. Therefore, when the platform has a pronounced slow drift motion, the effective travel of the heave compensator should be increased. Fig. 10 (c) shows the lateral displacement at the midpoint of the drill-string in the seawater segment. It can be found that in Case 2, the lateral stiffness of the drill-string in the seawater segment fluctuates due to changes in the axial tension, resulting in lateral displacement fluctuations, which is also the reason why the response result of Case 3 is slightly different from that of Case1. Fig. 10 (d) shows the bit rotation speed response. It can be seen that the response results of Case 3 and Case 1 are significantly different from those of Case 2, indicating that the bit-twist vibrations are greatly affected by the heaving motion of the platform, and slow-drift motion has some influence. Therefore, in order to simulate the actual drilling conditions as well as possible, it is necessary to analyze the response characteristics of the entire drill-string system considering the platform heave and slow drift motion. Figure 10. Response curves considering platform heave and drift separately and both: (a) axial displacement at the top of the drill-string, (b) WOB, (c) lateral displacement at the midpoint of the drill-string in the seawater section, and (d) bit speed 3.3. The influence of different slow drift motion amplitude For dynamically positioned drilling platforms, the period and amplitude of slow drift motion depends on the marine environment load, platform control system, position sensing system, and control algorithm (Qin et al., 2022). While the top of the drill-string is directly connected to the drilling platform, the long-term slow drift of the platform and the sea current load force the drill-string to generate large lateral displacements. Assuming that the platform is located at u=5m at the initial time, different slow drift amplitudes (10m, 15m and 20m) are taken to make the platform perform slow drift at 5-25m, 5-35m and 5-45m respectively, and the slow drift period is 80s. The amplitude of heave motion of the platform is 0.5 m, and its period is the same as the wave period of 11.2 s. Different from traditional land drilling, riserless drilling has no wall restrictions on the drilling string in seawater section. The lateral vibration in the seawater section is greatly 294 influenced by the slow drift movement of the platform and the combined excitation of the ocean current, while the influence on the drilling string and the drilling assembly in formation section is relatively small. Therefore, it is necessary to make a separate analysis of lateral vibration characteristics of drill- string in the seawater section. Fig. 11 shows the lateral time- history response diagram of drilling string in seawater segment. It can be seen that the drilling string exposed in seawater moves periodically in the lateral direction under the action of the platform's slow drift motion and ocean current. The larger the slow drift amplitude, the larger the transverse deformation, and the greater the bending moment effect of drilling string. Therefore, it is necessary to strictly control the slow drift motion of the rig during deep-water riserless drilling. In addition, it can be found that the drill-string motion slightly lags behind the platform, and the lag time increases with the increase of depth. Combined with the lateral displacement time-history diagram of the midpoint (depth 600 m) of the drill-string in the seawater section in Figure 11 (d), the platform returned to the initial position and the speed was zero at 80 s and 160 s, while the drill-string obviously lagged behind the platform, and the lag time increased with the increase of the slow drift amplitude. Figure 11 (e) is the lateral velocity curve of the midpoint of the drill- string in seawater section. One can find that the amplitude of the transverse velocity fluctuations increases with the amplitude of the slow drift motion. Figure 11. Lateral time-history response of drill-string in seawater section Fig. 12 shows the axial time-history response curve of the top of the drill-string. It can be seen from the time-history response of axial displacement at the top of drill column in Fig. 12 (a) that the maximum axial displacement increases with the increase of slow drift amplitude. This is because the drilling string in seawater segment is dragged by the platform, and the axial displacement is generated under the action of vertical and horizontal coupling. However, the slow drift amplitude of the platform is not proportional to the maximum axial displacement at the top of the drill-string. For a slow drift amplitude of 20 m, the maximum axial displacement is 3 times larger than for a slow drift amplitude of 10 m. Fig. 12(b) shows the difference between the axial displacement at the top of the drill-string and the platform, that is, the expansion of the heave compensator. It can be seen that the maximum expansion increases as the slow drift amplitude increases. Therefore, when the slow drift amplitude is large, the compensation cylinder stroke of the heave compensator should be increased. 295 Figure 12. Axial time-history response of the drill-string top: (a) displacement, and (b) relative displacement between the platform and the drill-string top Fig. 13 (a) shows the WOB time-history response diagram. It can be found that the amplitude of the fluctuation of the bit pressure increases with the amplitude of the slow drift. Fig. 13 (b) - (d) shows the response curve of WOB amplitude frequency. It can be seen that there are two main frequency components of the bit pressure, 0.0125 Hz and 0.0875 Hz, corresponding to the frequency of the slow drift and the heeling motion of the platform, respectively. Fig. 13 (b) - (d) shows the response curve of WOB amplitude frequency. It can be seen that there are two main frequency components of the bit pressure, 0.0125 Hz and 0.0875 Hz, corresponding to the frequencies of the slow drift and the heave motion of the platform, respectively. The amplitude of the corresponding frequency component is basically about 7.53kN because the heave motion is unchanged, while the amplitude of the corresponding frequency component of the slow drift motion increases with the increase of the slow drift amplitude, which is consistent with the axial displacement of the top of the drill- string. The amplitude is not proportional to the slow drift amplitude and the increment increases with the slow drift amplitude. It can be seen that the slow drift movement of the platform not only has a great impact on the drill-string in the seawater section, but also can make the downhole bit pressure fluctuate through the vertical and horizontal coupling, which further affects the vibration of the whole drill-string system. Figure 13. WOB time history curve and amplitude-frequency response curve Fig. 14 shows the torsional response of the bit, and Fig. 14 (a) - (c) shows the time-history response curves of friction torque, torsional displacement and rotational speed of the bit, respectively. It can be found that the bit changes back and forth between viscous and slip-off states, and the maximum rotational speed is more than 18 rad/s, which is nearly 4 times the rotational speed of the turntable. It can be found that the twist response of the bit varies for different values of the slow drift amplitude. Fig. 14 (d) - (f) shows the response curve of bit torsion amplitude and frequency, where the turntable rotary speed is 5 rad/s. It can be found that there is a dominant frequency component of 0.125 Hz, whose amplitude decreases with increasing slow drift amplitude. Fig. 14 (g) - (I) shows the phase diagram and the Pointcaré map of the bit torsion. It can be found that the bit rotation of the riser free drilling string is characterized by chaos, which is enhanced with the increase of the slow drift amplitude. For a value of the slow drift amplitude of 10 m, the phase diagram shows 296 multiple loops, while the other two do not have obvious fixed loop trajectories and have stronger nonlinear features. Figure 14. Bit torsional response: (a) friction torque, (b) displacement, (c) speed, (d)-(i) speed amplitude-frequency response, phase diagram and Poincaré map Fig. 15 shows the torsional response at the midpoint of the drill-string (depth 1250 m). it can be found that the response is more complex than that of the drill bit. From the amplitude frequency response curve of Figure 15 (b) - (d), it is found that there is an additional high-energy frequency component of 0.663 Hz. However, the amplitude of the 0.125 Hz 297 frequency component decreased greatly, and the frequency of Fig. 15 (b) and (d) decreased slightly to 0.121 Hz. It can be found from the torsive phase diagram and Poincaré map of the drill-string in Fig. 15 (e) - (g) that the torsive vibration also has chaotic characteristics and is enhanced with the increase of the slow drift amplitude. Figure 15. Torsion-response of drill-string at 1250m depth: (a) rotational speed, (b)-(g) amplitude-frequency response, phase diagram and Poincaré map 3.4. The influence of different slow drift motion frequencies The platform is assumed to be located at u = 5 m at the initial time with a slow drift amplitude of 15 m. The platform performs slow drift motion at 5-35m, respectively, and takes different slow drift periods of 40s, 80s and 120s. The platform heave motion is consistent with 2.3. Fig. 16 shows the lateral time-history response diagram of the drill-string in the seawater section. It can be found that the drill-string in the entire seawater section moves periodically along with the slow drift movement of the platform in the lateral direction. The phenomenon that the lateral movement of the drill-string lags behind the platform also exists. Comparing Fig. 16(a)-(c) and Fig. 16(d) at the midpoint of the drill*string in the seawater section (depth 600 m), it can be found that when the slow drift period is 40 s, the maximum lateral displacement of the drill-string at each depth is significantly smaller than the other two cases, which is due to the large movement speed of the platform, and the drill-string movement lags behind the platform and has not moved to the maximum displacement. It's dragged back by the platform. Fig. 16 (e) shows the curve of lateral velocity at the midpoint of the drill-string in the seawater section. It can be seen that the fluctuation amplitude of lateral velocity decreases with the increase of the period of slow drift movement. 298 Figure 16. Lateral time-history response of drill-string in seawater section Fig.17 (a) shows that the maximum axial displacement decreases with the increase of the slow drift period. Fig. 17 (b) shows the difference axial displacement between the top of the drill-string and the platform, which is also the expansion amount of the heave compensator. It can be found that the maximum amount of expansion increases as the slow drift period decreases. Therefore, when the period of slow drift motion is small and the slow drift velocity is large, the compensation cylinder stroke of the heave compensator should be increased. Figure 17. Axial time-history response of the drill-string top: (a) displacement, and (b) relative displacement between the platform and the drill-string top Fig. 18 (a) shows the WOB time-history response diagram. It can be found that the amplitude of the fluctuation of the bit pressure increases with the decrease of the slow drift period. Fig. 18 (b) - (d) shows the amplitude-frequency response curves of WOB. The amplitude of the corresponding frequency component is essentially about 7.53 kN due to the same heave motion, while the amplitude of the corresponding frequency component of the slow drift motion increases as the slow drift period decreases. Therefore, it is necessary to limit the slow drift motion period and the motion speed of the platform to reduce the down-hole bit pressure fluctuations. 299 Figure 18. WOB time history curve and amplitude-frequency response curve Fig. 19 shows the torsional response of the bit, and Fig. 19 (a) - (c) shows the time-history response curves of the friction torque, torsional displacement and rotational speed of the bit, respectively. It can be found that the bit changes back and forth between the viscous and slip states, and that the torsional response of the bit at different slow drift periods is somehow different. Fig. 19 (d) - (f) shows the response curve of bit torsion amplitude and frequency, where the turntable rotary speed is 5 rad/s. It can be found that the bit rotation velocities at different slow drift periods have the largest amplitudes around 0.12Hz, which are 4.3 rad/s, 5.42 rad/s and 5.1 rad/s, respectively. Fig. 19 (g) - (I) bit torsion phase diagram and Poincaré map. The absence of an obvious fixed toroidal trajectory in the phase diagram indicates that the bit rotation has a chaotic character and increases as the slow drift period decreases. 300 Figure 19. Bit torsional response: (a) friction torque, (b) displacement, (c) speed, (d)-(i) speed amplitude-frequency response, phase diagram and Poincaré map Fig. 20 shows the torsional response at the midpoint of the drill-string (1250 m depth). Compared with the drill bit, the response is more complex, because there is an extra high frequency component about 0.66Hz in the amplitude- frequency response curve of Fig. 20 (b) - (d). It can be seen from Fig. 20 (e) - (g) that the torsional vibration at the midpoint of the drill-string also has chaotic characteristics. Figure 20. Torsion-response of drill-string at 1250m depth: (a) rotational speed, (b)-(g) amplitude-frequency response, phase diagram and Poincaré map 301 3.5. The influence of different heave motion amplitude The amplitude of heave motion of the platform is 0.5m, 1m and 1.5m respectively, and its period is consistent with the wave period of 11.2s. The amplitude of the slow drift movement of the platform is 15m, the slow drift movement is performed at 5-35m, and the slow drift period is 80 s. Fig. 21 shows the time-history response of the lateral displacement at the midpoint of the drill-string in seawater segment. Compared with the lateral displacement curve without heave compensation device in Fig. 9 (c), the displacement fluctuation in Fig. 21 (a) is small, and the influence of different heave amplitude on the lateral displacement in seawater segment can be obviously observed from the velocity response curve. It can be found that as the amplitude of the heave increases, the amplitude of the axial tension of the drill-string increases, the amplitude of the transverse stiffness fluctuation increases, and the amplitude of the transverse velocity fluctuation of the drill-string in seawater also increases. Figure 21. Lateral time-history response of drill-string in seawater section Fig. 22 shows the axial time-history response curve of the top of the drill-string. It can be seen from the time-history response diagram of axial displacement at the top of the drill- string in Fig. 22 (a) that the amplitude of axial displacement fluctuation decreases with the decrease of heave amplitude. At the same time, it can be found that the slow drift movement of the platform is the main influencing factor for the axial displacement of the top of the drill-string, while the heave movement of the platform has relatively little influence because it is compensated by the heave compensator. Fig. 22 (b) shows the difference between the axial displacement at the top of the drill-string and the platform, that is, the extensibility of the heave compensator. It can be seen that the maximum extensibility increases with the increase of the heave amplitude. Therefore, when the heave amplitude is large, the compensation cylinder stroke of the heave compensator should be increased. Figure 22. Axial time-history response of the drill-string top: (a) displacement, and (b) relative displacement between the platform and the drill-string top Fig. 23 (a) is the response diagram of WOB time history. It can be found that platform heave movement has a significant impact on bit pressure, and the larger the heave amplitude, the larger the bit pressure fluctuation range. Fig. 23 (b) - (d) shows the amplitude-frequency response curve of bit weight. Since the slow drift motion is the same, there is little difference in the amplitudes of the corresponding frequency components, and the amplitude of the corresponding frequency component of the heave motion increases with the heave amplitude and is proportional. 302 Figure 23. WOB time history curve and amplitude-frequency response curve Figure 24. Bit torsional response: (a) friction torque, (b) displacement, (c) speed, (d)-(i) speed amplitude-frequency response, phase diagram and Poincaré map 303 Fig. 24 is the torsional response of the bit, and Fig. 24 (a) - (c) is the time-history response curve of the bit friction torque, torsional displacement and rotational speed, respectively. It can be found that for the heave amplitudes of 1.5m, the bit has a long time inversion torque, that is, the TOB is negative and the bit has a clear inversion phenomenon. Reversing the bit causes an abnormal working state of the bit, which can easily cause damage or even fall-off of the cutter, seriously affecting the life of the bit. In addition, the maximum drill speed is more than 4.6 times the rotary speed (23 rad/s), which will undoubtedly accelerate the wear of the drill bit. Combined with the bit torsional amplitude-frequency response curves in Fig. 24 (e) and (f), it can be found that the frequency component with the largest amplitude is 0.0875Hz, which is consistent with the heave frequency, that is, the drill-string torsional vibration resonates with the axial vibration. When the heave amplitude in Fig. 24 (e) and (f) is larger, the amplitude-frequency response has two more high-energy frequency components than when the heave amplitude in Fig. 24 (d) is smaller, which is also detrimental to the bit life. Fig. 24 (g) - (I) shows the twist phase diagram and Poincaré mapping of the bit. It can be found that when the heave amplitude reaches 1 m, the phase diagram rapidly expands into a ring and the Poincaré map points are relatively concentrated. For a heave amplitude of 1.5 m, the phase diagram continues to expand to the left and the minimum relative twist angle reaches -8 rad, that is, the bit twist angle displacement is much larger than the rotation angle, which is also responsible for the bit reversal. At the same time, the Poincaré map is more concentrated, indicating that the sites are undergoing approximately periodic rotations. Fig. 25 shows the torsional response at the midpoint of the drill-string (depth 1250 m); it can be found that the torsional vibration of the drill-string at the heave amplitude of 1 m and 1.5 m also resonates with the axial direction, and there is a larger frequency component of 0.0875 Hz. At the same time, from Fig. 25 (e) - (g), when the heave amplitude is larger, the phase trajectory expands to the left, which means that the drill-string has greater torsional deformation and is subjected to greater torque, which affects the drill-string life and drilling safety. As can be seen from the Poincaré map points, the drill- string rotates quasi-periodically for heave amplitudes of 1 m and 1.5 m, with stronger periodicity for larger heave amplitudes. Therefore, it is important to limit the heave amplitude of the platform during the drilling operation. Figure 25. Torsion-response of drill-string at 1250m depth: (a) rotational speed, (b)-(g) amplitude-frequency response, phase diagram and Poincaré map 3.6. The influence of different slow drift motion frequencies The platform heave has periods of 5.2 s, 8.2 s, and 11.2 s, respectively, with a heave amplitude of 1.5 m. The slow drift motion is consistent with 2.5 knots. Fig. 26 shows the time-history response of transverse displacement at the midpoint of drill-string in seawater section. When the amplitude of the heave is 1.5 m, there is little difference in the lateral displacement for different heave periods. From the velocity response curves, it is easier to observe the effect of different heave periods on the transverse motion of the drill-string in seawater. It can be found that the 304 transverse velocity fluctuates up and down with the heave motion at small amplitudes. Compared with Fig. 9 (c), the fluctuation range of the axial tension of the drill-string is effectively weakened under the action of the heave compensation device, so that the lateral stiffness of the drill- string in the seawater section changes less. Figure 26. Lateral time-history response of drill-string in seawater section Fig. 27 shows the axial time-history response curve of the top of the drill-string. As can be seen from the axial displacement of the top of the drill-string in Fig. 27 (a), the fluctuation amplitude of the axial displacement with a heave period of 5.2s is slightly larger, which is caused by the large heave motion frequency and heave speed, and the inertial force existing in the drill-string. It can be seen from Fig. 27 (b) that the heave frequency has little influence on the maximum shrinkage of the heave compensator. Figure 27. Axial time-history response of the drill-string top: (a) displacement, and (b) relative displacement between the platform and the drill-string top Fig. 28 (a) is the time-history response diagram of WOB, and it can be found that the WOB fluctuates sharply with the heave movement. Fig. 28 (b) - (d) shows the amplitude- frequency response curve of WOB. Since the slow drift motion is the same, there is little difference in the amplitudes of the corresponding frequency components, and the amplitude of the corresponding frequency component of the heave motion decreases as the period of the heave decreases. Figure 28. WOB time history curve and amplitude-frequency response curve 305 Fig. 29 shows the torsional time-history response curve of the bit, and Fig. 29 (a) - (c) shows the bit frictional torque, torsional displacement and rotational speed, respectively. It can be found that when the heave amplitude is 1.5m, the heave frequency has a significant influence on the stick slip vibration period of the bit. As the period of the heave decreases, the time of the inversion torque of the bit decreases, and the inversion phenomenon of the bit disappears. Only for the heave period of 11.2 s, the bit has a small amplitude inversion. Fig. 29 (d) - (f) shows the bit torsional amplitude- frequency response curve. It can be found that the frequency component with the largest amplitude corresponds to the heave frequency, indicating that the bit torsional direction resonates with the axial direction. The amplitude-frequency response curve with a heave period of 11.2 s has an additional high-energy frequency component compared to the two cases with smaller heave periods, which may affect the bit lifetime. Fig. 29 (g) - (I) shows the bit torsional phase diagram and Poincaré mapping. It can be found that the points of the Poincaré map of different heave periods are concentrated, indicating that the sites move quasi-periodically. With the increase of the heave period, the phase trajectory pattern expands into a ring, and the ring radius increases, which indicates that the relative angular displacement difference of the bit is larger, and the maximum rotation speed of the bit is larger, which will aggravate the damage of the bit. Figure 29. Bit torsional response: (a) friction torque, (b) displacement, (c) speed, (d)-(i) speed amplitude-frequency response, phase diagram and Poincaré map 306 Fig. 30 shows the torsional response at the midpoint of the drill-string (1250 m depth). It can be found that the torsional vibrations of the drill-string are also in resonance with the axial direction at different periods of the heave. At the same time, as shown in Fig. 30 (e) - (g), the phase trajectory expands with the increase of the heave period. Unlike the torsional oscillations of the bit, the Poincaré mapping points at 5.2 s are relatively scattered, indicating that the torsional oscillations have chaotic features. Figure 30. Torsion-response of drill-string at 1250m depth: (a) rotational speed, (b)-(g) amplitude-frequency response, phase diagram and Poincaré map 4. Conclusions and Recommendations In this paper, an ALT bidirectional coupling nonlinear drill-string element suitable for riserless drilling is established by energy method and finite element method, and a bidirectional coupling solution is proposed by using the Newmark-β method. The ALT coupling nonlinear vibration model of drilling string in deep-water riser-free drilling is established considering the platform heave and slow drift motion, wind current, ocean current and wave load, wall contact and drillbit rock interaction. The corresponding MATLAB calculation code is compiled. Field test results of Kyllingstad and Nessjøen (2009) and Teng (2017) and parameters of four example Wells in Legong, South China Sea are used to verify the model results. The influence of heave compensator on drill-string vibration is analyzed. The results show that heave compensator can effectively reduce bit pressure fluctuation, avoid skip drilling and bit reversal, and reduce drill-string lateral vibration in seawater section. The influence of platform heave and slow drift motion on drill-string vibration is analyzed. The results show that heave and slow drift motion should be considered simultaneously when simulating drill- string process. The influence of different heave amplitude and period, slow drift amplitude and period of the platform on drill-string vibration is analyzed. The results show that the drill-string moves periodically with the platform in the lateral direction, and the lateral fluctuation caused by the platform's heave motion is small due to the presence of heave compensator. The axial displacement of the top of the drill-string and the WOB are affected by both the slow drift and heave motion of the platform. The axial displacement of the top of the drill- string is mainly affected by the slow drift motion, and the WOB is mainly affected by the heave motion. When the amplitude, frequency of heave motion and amplitude of slow drift motion are large, the heave compensator should be added to compensate the cylinder stroke. The larger the amplitude of slow drift motion, the shorter the period, the stronger the nonlinear characteristics of drill-string and bit rotation. When the heave amplitude is large, it is easy to cause axial and torsion resonance, which is not conducive to the life of drill- string and bit. And when the heave amplitude and period are relatively large, it is easier to cause the drill-string torsion 307 deformation to increase, the torque to increase, and even lead to the bit reversal. So it is very important to control the heave motion of the platform. Compliance with Ethical Standards Conflict of interest The authors declare that they have no conflict of interest. Data availability The data used to support the findings of this study are included within the article. The processed data are available from the corresponding author upon request. Funding The project is supported by the National Natural Science Foundation of China (Grant No. 51875489) and Sichuan Science and Technology Program (Grant No. 2022YFQ0034). Appendix A1. TF The top of the drill-string is connected to the platform through the passive heave compensation device, which can be regarded as a gas-liquid spring on the top of the drill-string. The piston, piston rod, lower frame and hook of the compensator device are regarded as a concentrated mass at the top of the drill-string. The force acting on the heave compensation device is treated as a weightless spring with stiffness coefficient Tk and fluid viscous drag coefficient Tc (Woodacre et al., 2015; Niedzwecki and Thampi, 1988). Thus, the top tension at the upper end of the drill-string is expressed as: 0 0 0 0 2 0 0 1 T T r s z t s z t z c T T F F k w w c w w m w A P a k V c B w                             (25) where sw is the platform heave displacement, cA is the total cross-sectional area of compensator cylinders, 0P is the initial pressure, 0V is the total volume of compressed air in the pressure vessels and cylinders, a is the specific heat ratio of air , w  is the relative velocity between the cylinder and piston, the constants B and  are 2860 and 0.52, respectively, which are determined experimentally (Niedzwecki and Thampi, 1988). A2. hv , bv and ba The sea current is caused by tidal current and sea wind. The sea current velocity cwv caused by sea wind is proportional to the actual wind speed wv (Mao et al., 2016): =cw w wv k v (26) where wk is the dimensionless wind speed factor, the general value range is 0.024 ~ 0.05, and in this paper it is set to 0.03. According to Ekman's drifting theory, the current velocity can be expressed as: 1 7 1 7 h cw d h d M z H z v v v z M M H H z v v z M H                          (27) where M is Ekman depth which is generally determined from the range 200~300m. When the water depth exceeds M , the current velocity is considered negligible. According to Airy wave theory, the horizontal velocity and horizontal acceleration of water point velocity can be written as respectively (Wang and Liang, 2013):     2 e cos 2 e sin 2 b b k zb b b k zb b b h v k x t h a k x t               (28) where bh is the wave height, ( 2 / )bk   is the wave number in which  is the wavelength, and ( 2 / )T  is the wave circular frequency, T is the wave period. A3. )b （ and bW In the equation (22), the exponential decay friction coefficient of the drill can be determined by the following formula (Moharrami et al., 2021): ) ( ) b b c s c e             （ (29) where c and s are dynamic friction coefficient and static friction coefficient respectively,  is decay coefficients. According to different surfaces, the static friction coefficient can be greater than 1, and the dynamic friction coefficient is selected in the range 0~1. The weight on the drill-bit can be determined by the following formula: 0 e 0 e 0 0 e 0 0 ) ) 0 ) b b b b b b b b b b W k w w k w w W W k w w W         （ （ （ (30) where 0bW is the initial static weight on bit, bw is the longitudinal displacement of the bit, 0bw is the initial longitudinal displacement of the bit, and ek is the contact stiffness of the bit-rock. In this paper, according to the test data of the bit-rock contact stiffness presented in document (Liao, 2020), ek is taken as 20 kN/mm. 308 References [1] Liu, J., Liang, S., Chen, Y., & Zhang, M. (2023). Axial-lateral- torsional coupling nonlinear vibration of riserless drill-string considering both seawater and stratum parts. Applied Ocean Research, 140, 103742. [2] Baumart, A., 2000. Stick-slip and bit-bounce of deep-hole drillstrings. Journal of Energy Resources Technology 122(2), 78-82. [3] Christoforou, A.P., Yigit, A.S., 1997. Dynamic modelling of rotating drill-string with borehole interactions. Journal of Sound and Vibration 206(2), 243-260. [4] Cai, M.J., Mao, L.J., Xing, X.S., Zhang, H.Z., Li, J., 2022. Analysis on the nonlinear lateral vibration of drillstring in curved wells with beam finite element. Communications in Nonlinear Science & Numerical Simulation 104, 106065. [5] de Moraes, L.P.P., Savi, M.A., 2019. Drill-string vibration analysis considering an axial-torsional-lateral nonsmooth model. Journal of Sound and Vibration 438, 220-237. [6] Dykstra, M.W., 1996. Nonlinear drill string dynamics. The University of Tulsa. [7] Ghasemloonia, A., Geoff Rideout, D., Butt, S.D., 2015. A review of drillstring vibration modeling and suppression methods. Journal of Petroleum Science and Engineering 131, 150-164. [8] Guzek, A., Shufrin, I., Pasternak, E., Dyskin, A.V., 2015. Influence of drilling mud rheology on the reduction of vertical vibrations in deep rotary drilling. Journal of Petroleum Science and Engineering 135, 375-383. [9] Kyllingstad, A., Nessjøen P.J., 2009. A New Stick-Slip Prevention System. Society of Petroleum Engineers. [10] Kessai, I., Benammar, S., Doghmane, M.Z., 2022. Dynamic failure analysis and lifetime estimation of Tool-string in rotary drilling system under Torsional-Axial coupled vibrations. Engineering Failure Analysis 134, 106037. [11] Khoshroo, M., Eftekhari, M., 2022. Nonlinear dynamics of drill-string in horizontal well under axial load and supercritical regime of flowing fluid. International Journal of Non-Linear Mechanics 138, 103861. [12] Li, Z.F., Zhang C.Y., Song, G.M., 2017. Research advances and debates on tubular mechanics in oil and gas wells. Journal of Petroleum Science and Engineering 15, 194-212. [13] Liu, Y.S., Gao, D.L., 2017. A nonlinear dynamic model for characterizing downhole motions of drill-string in a deviated well. Journal of Natural Gas Science and Engineering 38, 466- 474. [14] Liu, J., Wang, J.X., Guo, X.Q., Dai, L.M., Zhang C., Zhu, H., 2022. Investigation on axial-lateral- torsion nonlinear coupling vibration model and stick-slip characteristics of drilling string in ultra-HPHT curved wells. Applied Mathematical Modelling 107, 182-206. [15] Li, W., Huang, H.L., Yu, F., Ni, H.J., Jiang, W., Zhang, X.B., 2020. Modeling and numerical study on drillstring lateral vibration for air drilling in highly-deviated wells. Journal of Petroleum Science and Engineering 195, 107913. [16] Liao, M.L, 2020. Nonlinear Dynamics of a Vibro-impact System for Indenter-rock Interaction. Journal of Mechanical Engineering 56(21), 121-130. [17] Moharrami M.J., Martins C.D., Shirih H., 2021. Nonlinear integrated dynamic analysis of drill-string under stick-slip vibration. Applied Ocean Research 108(1), 102521. [18] Morison, J.R., Johnson, J.W., Schaaf, S.A., 1950. The Force Exerted by Surface Waves on Piles. Journal of petroleum technology 2(5), 149-154. [19] Mao, L.J., Liu Q.Y., Zhou, S.W., Wang, G.R., Fu, Q., 2016. Deep water drilling riser mechanical behavior analysis considering actual riser string configuration. Journal of Natural Gas Science and Engineering 33, 240-254. [20] Qin, K., Di, Q., Zhou, X., He, Y., Wang, W., Chen, F., Zhang, H., 2022. Nonlinear dynamic characteristics of the drill-string for deep-water and ultra-deep water drilling. Journal of Petroleum Science and Engineering 209, 109905. [21] Ren, F.S., Wang B.J., Chen, S.L., 2018. Nonlinear modeling and qualitative analysis of coupled vibrations in a drill string. International Journal of Bifurcation and Chaos 28(10), 1850119. [22] Rajabali, F., Moradi, H., Vossoughi, G., 2020. Coupling analysis and control of axial and torsional vibrations in a horizontal drill-string. Journal of Petroleum Science and Engineering 195, 107534. [23] Robello, S., 2013. Modeling and analysis of drill-string vibration in riserless environment. Journal of Energy Resources Technology 135, 1-9. [24] Su, K.H., Wan, W., Liu, J.L., Qi C.W., 2013. Drilling Platform Motion Effects on the drill-string of Deepwater Riserless Operational Mode. Science Technology and Engineering 13(7), 1734-1739. [25] Tucker, W.R., Wang, C., 1999. An integrated model for drill- string dynamics. Journal of Sound and Vibration 224(1), 123- 165. [26] Teng, X.Q., Di, Q.F., Li, N., Chen, F., Zhou, B., Wang, M., 2017. Measurement and analysis of stick-slip characteristics of drill string in ultra-deep wells. Petroleum Drilling Techniques 45(2), 32-39. [27] Thakar, V., Nambiar, S., Shan, M., Sircar, A., 2018. A model on dual string drilling: on the road to deep waters. Modeling Earth Systems and Environment 4(2), 673-684. [28] Vromen, T., Detournay, E., Nijmeijer, H., Van de Wouw, N.. 2019. Dynamics of drilling systems with an antistall tool: effect on rate of penetration and mechanical specific energy. SPE Journal 24(05), 1982-1996. [29] Wada, R., Kaneko, T., Ozaki, M., Inoue, T., Senga, H., 2018. Longitudinal natural vibration of ultra-long drill-string during offshore drilling. Ocean Engineering 156, 1-13. [30] Wang, Y., Gao, D., Fang, J., 2015. Study on Lateral Nonlinear Dynamic Response of Deepwater Drilling Riser with Consideration of The Vessel Motions in Its Installation. CMC- Computers, Materials & Continua 48(1), 57-75. [31] Woodacre, J.K., Bauer, R.J., Irani, R.A., 2015. A review of vertical motion heave compensation systems, Ocean Engineering 104, 140-154. [32] Wang, S.Q., Liang, B.C., 2013. Wave Mechanics for Ocean Engineering. Ocean University of China Press. China. [33] Niedzwecki, J.M., Thampi, S.K., 1988. Heave compensated response of long multi-segment drill-string, Applied Ocean Research 10(4),181-190. [34] Xu, Q., Chen, G.M., Wang, G.D., Zhou, C.J., Yin, Z.M., 2011. Research and Application Prospect of Riserless Mud Recovery Drilling Technique. Drilling and Production Technology 34(01), 11-13. [35] Xing, Y.F., Liang, K., 2015. Nonlinear vibration analysis of longitudinal-transverse coupled beam. Journal of Beijing University of Aeronautics and Astronautics 41(8), 1359-1366. [36] Xia, K., Wang, Z., Chen, G., 2019. Assessment on Drilling Load for Riserless Drilling Technique. Drilling and Production Technology 43(5),16-19. 309 [37] Xie, K., Wang Y., Fug, T., 2019. Dynamic response of axially functionally graded beam with longitudinal-transverse coupling effect. Aerospace Science and Technology 85, 85-98. [38] Yin, Y.Q., 1987. Introduction to nonlinear finite element in solid mechanics. Peking University Press, Tsinghua University Press, Beijing, China. [39] Yang, J.H., Guo, X.X., Qiu, M.X., Sun, N.D., 2019. Disruptive innovations in the future oil and gas industry. World Petroleum Industry 26(6), 23-29. [40] Zhao, D., Hovda, S., Sangesland, S., 2016. The effect of stick slip vibration on the backward whirl of bottom hole assembly in drill-string, in: Proceedings of the ASME 35th International Conference on Ocean, Offshore and Arctic Engineering, V008T11A042.