管内可控铰接柔性结构双层接触非线性力学分析

外管内柔性结构是一个伴随机构与结构之间相互转化的三体两层梁梁接触非线性问题，柔性结构主要包括内杆(柔性杆)和弯外壳(导向管)，内杆为多段可控万向铰接结构，弯外壳为单侧开槽结构，管截面出现偏心，且多处单侧开槽的细长弯外壳与铰接柔性杆和井壁构成三体两层接触，如何合理描述柔性结构受力状态，揭示其柔度与载荷传递规律，对超短半径水平井钻井具有重要意义。

对于柔性铰接杆，人们开展了相关的力学分析。苏义脑[1]、赵俊平等[2]将铰链接头视为平面铰链，建立了铰接式柔性钻具的力学模型，通过静力学分析，推导出任一铰链处的受力状况，确定了钻压与大钩载荷之间的关系。帅健等[3-4]采用凝聚自由度方法建立了铰接式钻具组合的二维有限元模型，并对其进行了单因素影响分析。张绍槐等[5]、狄勤丰等[6]采用有限元迭代法研究了柔性短节位置对底部钻具动力学特性的影响，建立了钻柱动力学系统平衡方程；高加强[7]和LIU等[8]考虑柔性钻杆的结构特点和受力状态，建立了不同转角的柔性钻杆三维有限元模型，得到柔性钻杆所能承受的极限拉力和扭矩。

对于梁梁接触问题，WRIGGERS等[9]首先提出了用于研究圆截面梁接触动力学问题的无摩擦接触模型。此后，他们将该模型推广到含库仑摩擦的圆截面梁和矩形截面梁的接触问题中[10-11]，通过反演算法来计算摩擦力。WEYLER等[12]分析了罚函数法和拉格朗日法的优缺点，指出两种方法在处理梁梁接触问题时表现相似。LITEKA提出了梁梁摩擦接触的埃尔米特多项式平滑法[13]和多点接触法[14]。KHUDE等[15]结合绝对节点坐标方法与离散元法，提出了大变形细梁的多点接触模型，利用小球算法进行接触检测并利用并行计算提高效率。REPUPILLI[16]通过使两接触梁单元形成的四面体的体积等于零而提出了计算梁接触问题的鲁棒方法。NETO等[17]利用小球算法提出了梁的自接触模型。MEIER等[18-19]基于基尔霍夫梁理论开发了一种全新的全角梁接触公式，该公式在大接触角的范围内应用点接触公式，在小接触角范围内应用线接触公式。刘延强等[20]采用罚函数法对梁梁接触边界进行处理。刘巨保等[21-22]通过在节点位置放置间隙元，使梁梁之间形成一个连续体，来描述梁梁之间随机接触状态。庞东晓等[23]利用元胞自动机的基本理论，提出了一种基于元胞自动机的梁梁多向接触非线性有限元方程的求解方法。董世民等[24]应用可移动的双向弹簧元描述梁梁接触状态，建立了梁在梁内弯曲变形规律与接触压力的混合有限元仿真模型。

综上所述，人们对铰接连接的变截面柔性细长杆与外管的接触问题研究较少，且现有研究中，铰接处没有考虑相对转角限制，接触梁均是同心截面梁。鉴于此，本文对可控铰接柔性杆与偏心截面管和外管的双层接触非线性问题开展研究。

1 管内柔性结构双层接触模型的建立

假设外管的轴线为任意圆弧与直线的组合，外管内铰接柔性结构如图1所示。

将柔性杆、导向管和外管看作三维有限应变梁单元，单节柔性杆之间通过万向铰连接，导向管为多处单侧开槽结构，在外载作用下柔性杆伴随有机构与结构之间的相互转化。同时柔性杆与导向管、导向管与外管均产生随机接触，各自产生的接触力又作为外载影响着接触体对的接触行为。柔性杆、导向管及外管构成三体两层接触，两两接触行为相互耦合，单独的空间梁单元模型不能解决铰接柔性杆与导向管、导向管与外管边界的接触非线性问题。本文采用空间梁单元、万向铰单元、刚性梁单元和接触单元多种单元组合的有限元模型解决可控万向铰接柔性杆、开槽不完整导向管与外管接触问题。

1.1 可控万向铰接柔性杆的分析模型

采用有限元法，用两节点12个自由度的空间梁单元对柔性杆进行离散，整体平衡方程可表示为

式中，F为节点力列阵；K为总体刚度矩阵；δ为柔性杆节点位移列阵。

基于可控万向铰可在2个方向上发生相对转动的特点，通过2个节点对其进行定义，节点I和J重合，每个节点具有3个平动位移以及3个转动位移，万向铰单元坐标表示如图2所示，2个节点需定义局部直角坐标系来描述运动约束，y轴方向为传动轴方向。可控万向铰是在万向铰的基础上，在2个方向上只能相对转动固定的角度，当达到转动限制度数时，连接方式就变为刚性连接，其运动由机构问题转化为结构问题。

任意给定时刻万向铰单元的运动约束表示如下：

式中：uI、uJ为节点I和节点J的线位移向量；

均为沿各局部坐标系相应方向的单位向量。

节点I和节点J的局部坐标系相对位置的第一和第三卡登角[25]分别为

两局部坐标系相对转角的改变量Δφ和Δψ可表示为

式中，φ0、ψ0分别为节点I和节点J的局部坐标系相对位置的初始第一和第三卡登角。

当Δφ和Δψ达到转动限制度数时，万向铰单元的运动约束表示为

式中，θI、θJ分别为节点I和节点J的角位移向量。

由上述分析可知，万向铰单元由于节点位置是重合的，因此没有形函数，只需通过给定运动约束来实现其功能，将万向铰单元与梁单元的平衡方程相结合，经过坐标转换和拼装过程就能得到万向铰连接的梁单元的总体平衡方程：

式中,F1为铰接柔性杆节点力列阵；Kn为整体坐标系下万向铰连接的梁单元的总体刚度矩阵；δ1为铰接柔性杆节点位移列阵。

1.2 导向管的分析模型

导向管被离散为空间梁单元，基于导向管多处开槽的结构特点，分别定义完整圆截面1-1及偏心截面2-2，建立导向管单元的模型如图3所示。以梁单元ei和ei+1为例，ei单元的截面为1-1，ei+1单元的截面为2-2，两截面的内外径相等，ra为截面内径，rd为截面外径，θ为偏心截面的圆心角，单元ei由节点i和节点i+1组成，单元ei+1由节点i+2和节点i+3组成，在2个梁单元截面交界处，由于形心位置不同，导致交界处节点不共用，也就是图中节点i+1和i+2位置不同，要引入连接单元ej以满足交界处位移和载荷协调，显然，连接单元ej由节点i+1和i+2组成，且该单元有12个自由度，因此，连接单元采用刚性梁单元来实现拉压、弯曲、扭转等力素的传递[26]。

节点i+1和节点i+2的位移可分别表示为

任意给定时刻连接单元ej满足的关系式为

导向管由多个相同结构的子结构连接而成，各子结构之间同样通过连接单元ej+1建立连接。

根据上述分析，将刚性梁单元与梁单元的平衡方程相结合，经过坐标转换和拼装过程就能得到具有偏心截面的变截面梁总体平衡方程：

式中,F2为导向管的节点力列阵；K2为梁单元总体刚度矩阵；Kd为连接单元的总体刚度矩阵；δ2为导向管节点位移列阵。

1.3 三体两层接触分析模型

为分析柔性杆、导向管和外管的三体两层接触问题，创建接触单元和目标单元，分别依附在柔性杆的外表面和导向管的内表面以及导向管的外表面和外管的内表面，接触单元和目标单元与所依附的梁单元满足变形协调条件，当接触体对未发生接触时，它不影响柔性杆和导向管的运动；当接触时，其刚度为某一给定值，阻抗接触体之间的相互侵入。

当发生接触时，接触面有法向接触压力，若接触体发生相对滑动，则有轴向摩擦阻力存在，本文采用经典的库仑摩擦模型，表示如下：

式中,P为法向接触压力；μ为摩擦因数。

圆截面梁是否发生接触，可通过引入间隙函数h来定义两物体的接触状态，当间隙函数小于零时可判定接触存在。梁在梁内接触模型如图4所示。

设柔性杆与导向管之间的间隙函数为h1，导向管与外管之间的间隙函数为h2，间隙函数可由几何关系确定为

式中,ra、rd分别为导向管的内外半径；rw为外管的内径；rc为柔性杆的外径；ζ1、ζ2分别为柔性杆中心线与导向管中心线的最小距离，以及导向管中心线与外管中心线之间的最小距离。

接触体对在接触边界上满足变形协调，基于接触表面沿法线方向的接触间隙，采用增广拉格朗日算法计算法向接触力。

通过以上分析，综合可控铰接连接的机构问题及接触非线性问题，结合梁单元、万向铰接单元和梁梁接触单元，可得到实现分析导向管内可控万向铰接柔性杆接触非线性的整体平衡方程：

式中，Ft1为柔性杆与导向管之间的接触力列阵，是节点位移的函数；Kt1为柔性杆与导向管之间的接触刚度矩阵。

同理，结合梁单元、连接单元和梁梁接触单元，可得到实现分析带有双层接触的偏心截面导向管的整体平衡方程：

式中，Ft2为导向管与外管之间的接触力列阵；Kt2为导向管与外管之间的接触刚度矩阵。

式(14)与式(15)通过柔性杆与导向管之间的接触力Ft1耦合，由于第一次求得的解通常不能满足接触判定条件，需要反复迭代计算，直到得出全部满足接触判别条件的解为止，因此，选用修正的牛顿辛普森法进行求解。

2 模型验证

为了验证可控万向铰接模型和双层接触分析模型以及求解算法的正确性，建立图5所示的平面模型，该模型基于小变形假设。内梁中间有一个可控铰连接，可控铰的转动限制度数为β，内梁直径为d，内梁与中间梁的环空间隙为Δ1，中间梁与外梁的环空间隙为Δ2，中间梁的厚度为t，内梁与中间梁的弹性模量均为E，惯性矩分别为 I1和I2，杆长为2l，内梁左端全固定，右端为可动铰支座，中间梁左端铰支，外梁全固定，内梁铰接点作用集中载荷F。

先以内梁为研究对象，设铰接处刚好达到转动限制度数时，临界集中力为Fr，铰接点在Fr作用下产生的挠度为δr，要使内梁铰接处达到转动限制度数后再与中间梁发生接触，应使Δ1≥|δr|，那么当F取值不同时，3根梁的接触状态不同。设内梁C点刚好与中间梁发生接触时的临界力F=F1，此时，C点挠度为-Δ1；中间梁刚好与外梁接触时的临界力F=F2，此时，内梁铰接点C处的挠度为-(Δ1+Δ2)，中间梁对应点C处的挠度为-Δ2。令内梁与中间梁的接触力为FC，中间梁与外梁的接触力为F′C，铰接点的弯矩为MC，挠度为δC，那么根据3根梁的接触状态，将其分为4种情况，如表1所示，其中虚线为梁在集中载荷作用下的弯曲变形。根据叠加原理、平衡方程及变形协调条件，可求得F在不同取值范围内铰接点的挠度和弯矩以及梁梁之间的接触力。需要指出的是，当F>F2时，中间梁与外梁的接触实际上为线接触，为了计算方便，令F略大于F2，假设中间梁与外梁为点接触，接触边界仍处理为集中接触力，将所求参量的表达式汇总如表1所示。

令E=210 MPa，d=20 mm，l =1 m，β=1 °，t=3 mm，Δ1=8 mm，Δ2=3 mm，可求得Fr=34.48 N，F1=57.49 N，F2=390.96 N，为了得到3根梁的4种接触状态，F分别取20 N、45 N、200 N和400 N，将数值计算方法得到的结果与表1所列的理论解进行对比，误差分析如表2所示。由表2可知，数值解与理论解最大误差仅为3.8%，验证了可控铰接模型及双层接触模型的可行性。

3 工程应用

采用管内柔性结构的双层接触模型和计算方法对超短半径水平井造斜段柔性钻具进行接触非线性分析。

造斜钻具总成由柔性钻杆、导向管、钻具三连接和钻头组成，其中钻具三连接由钻具上三连接和钻具下三连接构成，钻具上三连接将常规钻杆、柔性杆和导向管连接在一起，钻具下三连接将钻头、柔性钻杆和导向管连接在一起。在造斜钻进时，常规钻杆将轴向力和扭矩传递给钻具上三连接，钻具上三连接将轴向力传递给导向管，将扭矩传递给柔性钻杆，从而实现载荷的分类传递；钻具下三连接将导向管轴向力和柔性钻杆扭矩传递给钻头，从而实现不同载荷的集中传递。

钻具上三连接是柔性钻具的关键部件，为内外嵌套式结构，内轴和外套通过只传递轴向载荷的轴承连接，以达到造斜钻进时内轴只传递扭矩、外套只传递轴向力的目的。

钻具上三连接端部有与常规钻杆相匹配的锥螺纹接头，在钻具上三连接关键连接部位采用O形圈进行密封。钻具上三连接与导向管之间采用螺纹连接，与柔性钻杆之间采用键连接。上三连接结构如图6所示。

根据造斜段柔性钻具的钻进过程，对曲率半径ρ=3.2 m、造斜钻进井斜角γ=60°的柔性钻具进行力学分析，柔性钻具结构如图7所示。

3.1 柔性钻具力学模型的建立

选取造斜钻具总成和井壁为研究对象，柔性钻杆、导向管、三连接和钻头的弹性模量为210 GPa，泊松比为0.3，密度为7 850 kg/m3，屈服极限为835 MPa。导向管的结构尺寸如图8所示。

柔性钻杆可控铰接连接部件结构如图9a所示，铰接处限制度数为3°，由于研究的重点在于细长柔性钻具与井壁的接触状态，并非工具的强度，故为了提高计算效率，考虑柔性钻杆的内部结构，将其简化为相应截面尺寸的圆柱结构，简化后的结构如图9b所示。单节柔性钻杆各截面S1～ S5对应的长度分别为33 mm、28 mm、38 mm、32 mm和27 mm，对应的外径分别为54 mm、84 mm、84 mm、91 mm和69 mm，内径均为30 mm，铰接位置位于截面S2和截面S3之间。上下三连接内径为30 mm，外径为110 mm，长度为ls=80 mm，钻头内径为30 mm，外径为118 mm，长度为lz=76 mm，井壁内径为124 mm，外径为140 mm，柔性钻具总长度l=ρπ/2+2ls+lz。柔性钻杆与导向管的摩擦因数取0.1，导向管、三连接、钻头与井壁之间的摩擦因数取0.3。

柔性钻杆、导向管和井壁均采用梁单元进行离散，导向管圆截面梁与偏心截面梁之间采用刚性梁连接，柔性钻杆铰接连接采用万向铰单元进行离散，由于单节柔性钻杆有5个不同的截面，因此，将单节柔性钻杆划分为5个梁单元和1个万向铰单元。单节导向管有2个不同的截面，且截面偏心，将单节导向管划分为2个梁单元和2个刚性梁单元，井壁划分为梁单元且单元长度与导向管相对应。

3.1.3 造斜段柔性钻杆万向铰初始卡登角的确定

根据井眼曲率半径和柔性钻杆结构参数可确定柔性钻杆通过弯曲井眼时万向铰处的弯曲角度(初始第一和第三卡登角)，而万向铰处的弯曲角度由局部坐标系定义。设井眼轴线所在的平面为XY平面，那么各万向铰处初始第一卡登角(绕X轴方向)均为零，柔性钻杆万向铰处局部坐标系的确定如图10所示。

设第n个铰由节点In和节点Jn组成，且两节点重合，图中

和

分别表示节点In和节点Jn 的局部坐标系，ey方向为传动轴方向，ex与ey在XY平面内相互垂直，ez垂直于XY平面，ξIn和ξJn分别表示节点In和节点Jn的ey方向与整体坐标系Y轴方向的夹角，也就是节点In和节点Jn的初始第三卡登角(绕Z轴方向)，θ为铰接处与曲率中心的连线与水平方向的夹角。

单节柔性钻杆的长度为a，且满足

单节柔性钻杆长度相对于曲率半径值很小，可认为单节柔性钻杆对应的弧长等于杆长a。以造斜段起点为基准，第一个铰与造斜段起点连线的弧长应为l1+l2，以单节柔性钻杆长度a为弦对应的中心角为φ，弦与半径的夹角为η。下面根据几何关系进行分析。

当n=1时，在△OAB中：

根据式(17)和式(18)可确定各个万向铰的初始第三卡登角，也就确定了万向铰处的局部坐标系。

3.1.4 载荷及边界条件

为真实地描述上下三连接载荷传递功能，根据上下三连接的结构尺寸，以上三连接为例，将上三连接简化为2个不同截面、相同长度的梁，分别与柔性钻杆和导向管顶端相连，与柔性钻杆顶端相连的三连接外径等于与导向管顶端相连的三连接的内径。下三连接的简化方法与上三连接相同，钻头与导向管下三连接相连。根据柔性钻具的受力特点，将边界条件设为：耦合柔性钻杆和导向管的上三连接的三个方向线自由度，钻头处的弯曲角自由度不受限制，横向自由度为有初始间隙的接触边界[27]，扭转和轴向自由度为固定边界；柔性钻杆下三连接的边界条件与钻头处的边界条件一致，井壁全约束，柔性钻杆与导向管之间为接触摩擦边界，导向管与井壁、三连接与井壁以及钻头与井壁之间为接触摩擦边界。在柔性钻杆上三连接处施加T′=2 kN·m的扭矩，在导向管上三连接处施加F′=40 kN的轴向力，柔性钻杆、导向管、三连接和钻头承受自重q。柔性钻具力学模型如图11所示。

3.2 柔性钻具力学分析结果

通过对柔性钻具的力学分析，得到柔性钻杆在载荷作用下发生变形。为了研究可控万向铰对柔性钻杆变形及载荷传递规律的影响，将柔性钻杆变形放大30倍，得到柔性钻杆在不同坐标平面下的变形形态如图12所示，柔性钻杆与导向管之间的接触压力如图13所示，柔性钻杆万向铰处第三卡登角和第一卡登角如图14所示，柔性钻杆的扭矩和轴力分布规律如图15和图16所示。

由图12～图16可知，柔性钻杆顶端向X轴正方向倾斜，这是由于导向管在轴向力的作用下，割缝闭合，使其产生侧向弯曲，而柔性钻杆与导向管通过上三连接将顶端节点线位移自由度耦合，三个方向的线位移相等，因此柔性钻杆顶端会随着导向管产生侧向弯曲。柔性钻杆在XY平面内竖直段的“折线形”变形相对于造斜段较明显，这是由于万向铰竖直段初始第三卡登角为0，造斜段初始第三卡登角不为0，在3°的转动限制下，竖直段柔性钻杆的相对转角大于造斜段的相对转角，因此变形越明显。而柔性钻杆在YZ平面内竖直段和造斜段的“折线形”变形都很明显，这是由于竖直段和造斜段初始第一卡登角均为0，相对转角均较大。柔性钻杆与导向管之间的接触力呈随机分布状态，只有柔性钻杆最大截面处与导向管发生接触。竖直段柔性钻杆的扭矩在2 kN·m附近波动，且变化不大。造斜段柔性钻杆的扭矩整体呈逐渐减小的趋势，这是由于竖直段的万向铰卡登角正负值相当，也就是柔性钻杆扭矩增大的值与减小的值基本相等，而造斜段第三卡登角均为正值。柔性钻杆的轴力出现了中和点，这是由于柔性钻杆与导向管缩短的距离相等，当导向管割缝完全闭合时，柔性钻杆还有没达到转动限制度数的万向铰，柔性钻杆同时又承受自重作用，因此导向管顶端对柔性钻杆起支撑作用，导致柔性钻杆上部分受拉，下部分受压。

将导向管变形放大30倍，得到导向管变形形态如图17所示，导向管等效应力云图见图18，导向管与井壁之间的接触力如图19所示，导向管轴力分布规律如图20所示。

由图17～图20可知，导向管顶端上三连接向X轴正方向倾斜，这是由于导向管在轴向力的作用下，割缝闭合产生侧向弯曲。导向管偏心截面位置的等效应力大于完整截面位置，最大等效应力位于导向管最下端偏心截面位置处，其值为725 MPa，满足强度要求。导向管上三连接与井壁发生接触，上三连接附近的导向管不与井壁发生接触，随着竖直段长度的增大，导向管与井壁再次发生接触，接触力先增大后减小又趋于平稳，然后又逐渐增大，由于钻头的支撑作用，与钻头相连的下三连接及其附近的导向管不与井壁发生接触；竖直段导向管轴力逐渐增大，造斜段导向管顶端附近轴力先平稳，随着长度的增大，轴力逐渐减小，后又有增大的趋势。

4 结论

(1)基于柔性结构的内杆为可控万向铰接的结构特点构造了万向铰单元并给出运动约束；根据柔性结构弯外壳多处单侧开槽的结构特点，定义了开槽不完整的偏心梁截面，将其通过刚性梁单元与完整截面梁建立连接；综合可控铰接结构的机构问题及接触非线性问题，将梁单元、万向铰接单元、刚性梁单元和梁梁接触单元进行组合，建立了可控万向铰接柔性结构与外管双层接触非线性有限元模型，并选取修正的牛顿辛普森法进行数值求解。

(2)建立了具有一个可控铰接的内梁、中间梁和外梁双层接触平面模型，根据叠加原理、平衡方程和变形协调条件，推导了内梁在可控铰接处施加不同横向载荷下铰接处的挠度和弯矩以及各梁之间的接触力，将理论解与数值解进行对比，最大误差仅为3.8%，验证了外管内可控铰接柔性结构的接触模型和数值计算方法的可行性。

(3)铰接柔性钻杆的变形相对于初始构型呈现“折线形”的变形，铰接处多数达到转动限制度数；柔性钻杆与导向管之间的接触力呈随机分布状态，且只有柔性钻杆最大截面处与导向管接触；随着长度的增大，竖直段柔性钻杆的扭矩变化不大，造斜段柔性钻杆的扭矩逐渐减小。

(4)竖直段导向管与井壁之间的接触力小于造斜段，随着长度的增大，竖直段导向管轴力逐渐增大，造斜段导向管轴力逐渐减小。

[1] 苏义脑，赵俊平. 短半径水平井铰链式井下动力钻具组合力学模型的建立与分析[J]. 石油钻采工艺， 1995，17(4)：25-33.

SU Yinao, ZHAO Junping. Establishment and Analysis of Mechanical Model of Articulated Downhole Motor Assembly for Short Radius Horizontal Well[J]. Oil Drilling & Production Technology, 1995, 17(4): 25-33.

[2] 赵俊平, 王晓华. 短半径水平井铰接式钻柱串整体力学分析[J]. 东北石油大学学报, 1998, 22(3): 26-27.

ZHAO Junping, WANG Xiaohua. Global Force Analysis of Hinged Drill String in Short Radius Horizontal Well [J]. Journal of Northeast Petroleum University, 1998, 22(3): 26-27.

[3] 帅健，刘春. 铰接式钻具组合的动力学分析[J]. 石油矿场机械，2001，30(5)：15-18.

SHUAI Jian, LIU Chun. Dynamic Analysis of Hinged Drill String Assembly[J]. Oil Field Machinery, 2001, 30(5): 15-18.

[4] 帅健，刘春. 新型铰接式钻具组合特性分析[J]. 中国海上油气，2001，13(5)：35-38.

SHUAI Jian, LIU Chun. Static Analysis of New Type of Hinged Bottom Hole Assembly [J]. China Offshore Oil & Gas, 2001, 13(5): 35-38.

[5] 张绍槐，狄勤丰. 用旋转导向钻井系统钻大位移井[J].石油学报，2000，21(1)：76-80.

ZHANG Shaohuai, DI Qinfeng. Drilling Extended Reach Well with Rotary Steering Drilling System[J]. Acta Petrolei Sinica, 2000, 21(1): 76-80.

[6] 狄勤丰，王明杰，胡以宝，等. 柔性短节位置对带旋转导向工具底部钻具组合动力学特性的影响[J]. 中国石油大学学报(自然科学版)，2012，36(5)：84-88.

DI Qinfeng, WANG Mingjie, HU Yibao, et al. Effect of Flex Sub’s Position on Bottom Hole Assembly with Rotary Steering Tool[J]. Journal of China University of Petroleum (Edition of Natural Science), 2012, 36(5): 84-88.

[7] 高加强. 超短半径水平井柔性钻具关键部件结构设计及力学分析[D]. 大庆：东北石油大学，2012.

GAO Jiaqiang. Structural Design and Mechanical Analysis of Key Components of Flexible Drilling Tools for Ultra-short Radius Horizontal Wells[D]. Daqing: Northeast Petroleum University, 2012.

[8] LIU H, HUANG S, SUN Q, et al. Application of Ultra Short Radius Lateral Drilling Technique in Top Thick Reservoir Exploitation after Long Term Water Flooding[C]//Abu Dhabi International Petroleum Exhibition & Conference. Abu Dhabi, 2016: 1-11.

[9] WRIGGERS P, ZAVARISE G. On Contact between Three-dimensional Beams Undergoing Large Deflections[J]. Communications in Numerical Methods in Engineering, 1997, 13(6): 429-438.

[10] ZAVARISE G, WRIGGERS P. Contact with Friction between Beams in 3-D Space[J]. International Journal for Numerical Methods in Engineering, 2000, 49(8): 977-1006.

[11] LITEWKA P, WRIGGERS P. Contact between 3D Beams with Rectangular Cross-sections[J]. International Journal for Numerical Methods in Engineering, 2002, 53(9): 2019-2041.

[12] WEYLER R, OLIVER J, SAIN T, et al. On the Contact Domain Method: a Comparison of Penalty and Lagrange Multiplier Implementations[J]. Computer Methods in Applied Mechanics & Engineering, 2012, 205(1): 68-82.

[13] LITEWKA P. Hermite Polynomial Smoothing in Beam-to-beam Frictional Contact[J]. Computational Mechanics, 2007, 40(5): 815-826.

[14] LITEWKA P. Enhanced Multiple-point Beam-to-beam Frictionless Contact Finite Element[J]. Computational Mechanics, 2013, 52(6): 1365-1380.

[15] KHUDE N, MELANZ D, STANCIULESCU I, et al. A Parallel GPU Implementation of the Absolute Nodal Coordinate Formulation with a Frictional/contact Model for the Simulation of Large Flexible Body Systems[C]//ASME Conference on Multibody Systems and Nonlinear Dynamics. Baltimore, 2011: 1-11.

[16] REPUPILLI M. A Robust Method for Beam-to-beam Contact Problems Based on a Novel Tunneling Constraint[J]. Biological Control, 2012, 27(3): 300-309.

[17] NETO A G, PIMENTA P M, WRIGGERS P. Self-contact Modeling on Beams Experiencing Loop Formation[J]. Computational Mechanics, 2015, 55(1): 193-208.

[18] MEIER C, POPP A, WALL W A. A Finite Element Approach for the Line-to-line Contact Interaction of Thin Beams with Arbitrary Orientation[J]. Computer Methods in Applied Mechanics & Engineering, 2016, 308: 377-413.

[19] MEIER C, WALL W A, POPP A. A Unified Approach for Beam-to-beam Contact[J]. Computer Methods in Applied Mechanics & Engineering, 2017, 315: 972-1010.

[20] 刘延强，吕英民. 环空钻柱结构三维非线性分析[J]. 应用数学和力学，1994，15(3)：259-272.

LIU Yanqiang, LYU Yingmin. Three-dimensional Nonlinear Analysis of Drill String Structure in Annulus[J]. Applied Mathematics and Mechanics, 1994, 15(3): 259-272.

[21] 刘巨保，岳欠杯，李治淼，等. 井筒内杆管柱双层接触有限元分析及应用[J]. 力学与实践，2009，31(2)：35-39.

LIU Jubao, YUE Qianbei, LI Zhimiao, et al. The Analysis and Application of Double Layer Contact Finite Element about Rod and Tubing String in Well Bore[J]. Mechanics in Engineering, 2009, 31(2): 35-39.

[22] 刘巨保，栾绍信，张学鸿. 水平井压裂管柱受力变形分析的间隙元法[J]. 石油学报，1994(1)： 135-140.

LIU Jubao, LUAN Shaoxin, ZHANG Xuehong. Mechanics and Deformation Analysis for the Fracturing String in a Horizontal Well Using Gap Element Method[J]. Acta Petrolei Sinica, 1994, 15(1): 135-140.

[23] 庞东晓，刘清友，孟庆华，等. 三维弯曲井眼钻柱接触非线性问题求解方法[J]. 石油学报，2009，30(1)：121-124.

PANG Dongxiao, LIU Qingyou, MENG Qinghua, et al. Solving Method for Nonlinear Contact Problem of Drill Strings in 3D Curved Borehole[J]. Acta Petrolei Sinica, 2009, 30(1): 121-124.

[24] 董世民，张万胜，张红，等. 定向井有杆抽油系统杆管分布接触压力的研究[J]. 工程力学，2011，28(10)：179-184.

DONG Shimin, ZHANG Wansheng, ZHANG Hong, et al. Research on the Distribution Contact Pressure between Sucker Rod and Tubing String of Rod Pumping System in Directional Wells[J].Engineering Mechanics, 2011, 28(10):179-184.

[25] SHAH S V, SAHA S K, DUTT J K. Euler-Angle-Joints(EAJs)[J]. Dynamics of Tree-Type Robotic Systems, 2013, 62: 27-55.

[26] LAULUSA A, BAUCHUAU O A, CHOI J-Y, et al. Evaluation of some Shear Deformable Shell Elements[J]. International Journal of Solids and Structures, 2006, 43(17): 5033-5054．

[27] 刘巨保，岳欠杯. 石油钻采管柱力学[M]. 大庆：石油工业出版社，2011：18-19.

LIU Jubao, YUE Qianbei. Oil Drilling and Production String Mechanics [M]. Daqing: Petroleum Industry Press, 2011: 18-19.