常旭, 杨东超, 孙可平, 朱衡, 杨淇耀
清华大学 机械工程系, 北京 100084

收稿日期：2018-11-28
作者简介：常旭(1994-), 男, 硕士研究生
通信作者：杨东超, 副研究员。E-mail:ydc@tsinghua.edu.cn

摘要：轮式水平井爬行器依靠驱动轮与管壁的接触产生牵引力。目前已有的相关研究均忽略了接触导致的管壁塑性变形压痕，简单地利用摩擦因数的经验值来处理力的关系，导致结果误差较大。该文针对压痕形成过程利用滑移线理论，在楔形压头压痕模型的基础上，建立了正压力和转矩同时作用下的正压（轮齿对称轴线垂直于管壁）压痕计算模型。有限元仿真反映了在不同受力条件下塑性变形区域的变化规律，也验证了主动力与压痕形状的定量关系。压痕测量试验证明了理论模型用于计算驱动轮压入深度和压痕隆起高度的正确性，为今后的驱动轮斜压建模、牵引力建模等研究奠定了基础。
关键词：滑移线理论水平井爬行器塑性接触模型压痕试验
Analysis of the normal indentation process caused by tractor wheels
CHANG Xu, YANG Dongchao, SUN Keping, ZHU Heng, YANG Qiyao
Department of Mechanical Engineering, Tsinghua University, Beijing 100084, China


Abstract: The traction force of wheeled horizontal-well tractors is mainly produced by the contact between the wheels and the pipe walls. Previous studies ignored the plastic deformation caused by the contact between the wheels and the walls and used empirical friction coefficients to model the forces, which has led to significant errors. This paper derives a wedge indentation model (the wedge symmetry axis is perpendicular to the pipe wall) based on the simultaneous action of the normal force and the torque to predict the indentation formation according to slip line theory. Finite element analysis results not only show the plastic deformation for various forces, but also verify the quantitative relationship between the forces and the indentation shape. Comparisons with measured indentation shapes shows that the theoretical model can predict the indentation depth and the pile-up height, which can be used to improve sideling indentation modeling of wheels for better traction modeling.
Key words: slip line theoryhorizontal-well tractorplastic contact modelindentation testing
轮式爬行器^[1-4]主要用于向油井水平段输送仪器，其牵引力依靠驱动轮的楔形轮齿与管壁之间的接触产生。目前针对爬行器与管壁接触的研究工作相对较少，相关文献简单地利用摩擦因数的经验值来处理正压力与牵引力的关系^[5-11]。然而，在实际工作中由于正压力和牵引力可达3 000 N以上，驱动轮会在管壁产生塑性变形压痕。这就导致接触状态与传统的干摩擦区别较大，牵引力的理论上界也存在差异。刘清友等^[10]和秦浩^[11]考虑到该情况并利用有限元法进行了分析，虽然提出了压痕区域面积计算方法，但是没有建立完整受力条件下的理论模型，无法进一步对滚动接触进行分析和对牵引力进行计算。因此，根据实际发生的塑性接触过程建立新的正压理论模型十分必要。
滑移线理论最早用于求解满足平面应变假设和理想刚塑性假设的塑性接触问题。Hill^[12]和Johnson^[13]介绍了该理论的原理，并用于求解平板刚性压头压入问题。Grunzweig等^[14]在Hill等提出的无摩擦模型的基础上考虑了不同摩擦因数的影响，建立了平均压力与楔形压头顶角的关系，并与试验结果进行了对比。Rice^[15]对平板压头的压入问题进行了计算，其中考虑了材料各向异性的影响。Bai等^[16]考虑了滑移场中可能出现的应力不连续的现象，重新计算了钝角楔形压头的压入问题。
有限元等数值方法可以有效地求解各种非线性问题，被广泛用于处理弹塑性接触过程。但是由于计算时间比较长、计算与分析效果依赖使用者的经验，故该方法常被用于理论模型的验证及优化。例如，Komvopoulos等^[17]利用有限元结果标定了球形压头平均压力分布和压深关系式中的各系数值。Taljat等、Hernot等利用有限元结果分析了球形压头^[18-19]、锥形压头^[20]压痕两侧隆起高度的变化规律。Jackson^[21]利用有限元结果与滑移线理论得到的结果进行对比，得到了球形压头的平均压力方程。Bomidi^[22]分析了滚动过程弹塑性接触状态。王战江等^[23-26]则利用影响系数的半解析方法来求解弹塑性接触的问题，其中使用了最速下降法和快速Fourier变换，计算速度比有限元方法有了较大提高。但以上模型均未考虑压头受到转矩作用的情况，难以直接用于分析爬行器的驱动轮压痕。
本文针对爬行器驱动轮同时受正压力和转矩作用的实际情况，利用滑移线理论推导了正压情况下压痕形状与接触反力之间的关系，并利用有限元模型和实际压痕试验对该模型进行验证，以期为驱动轮斜压建模、牵引力建模和爬行器结构的优化设计等提供理论依据和奠定基础。
1 驱动轮力学模型水平井爬行器驱动轮受力如图 1所示。图 1a表示原动件施加的对管壁的正压力F和驱动轮转矩M; 图 1b表示管壁对驱动轮的接触反力，可将其分解为切向力F_T和法向力F_N。

图 1 爬行器驱动轮受力示意图

图选项

仅考虑正压情况，根据力的平衡可以得到

$\left\{ {\begin{array}{*{20}{l}}{F = {F_{\rm{N}}}, }\\{M = {F_{\rm{T}}}R.}\end{array}} \right.$

(1)

其中R表示驱动轮半径。
2 轮齿仅受到正压力的压痕计算模型驱动轮齿与管壁的接触实际上是三维的弹塑性变形问题，涉及到接触非线性和材料非线性两个难点，目前尚无有效的理论模型可以精确刻画其中的规律。为了简化分析，作出如下几点假设：
1) 由于驱动轮的材料为20CrMoTi，屈服强度(885 MPa)和硬度(217 HB)大幅高于试验用6061铝板，故假设轮齿在接触过程中不发生变形；
2) 由于轮齿厚度为10 mm，压痕深度为0.5 mm左右，二者相差较大，故假设在齿厚方向压力分布不发生变化，故而转化为二维平面应变问题；
3) 经过初步仿真表明，达到实际压深时弹性变形部分占比较小(1%~2%)，故假设管壁材料为理想刚塑性；
4) 由于摩擦力影响因素(各点接触状态、轮齿表面粗糙度等)较多，此处暂不考虑摩擦力的影响；
由此便可利用滑移线理论来处理该问题。首先，针对仅受正压力作用的情况，根据式(1)可知，F_N＞0, F_T=0。建立压痕附近的滑移线场：α线和β线如图 2所示。其中p为接触部分的压力分布。

图 2 正压情况滑移线场

图选项

根据Hencky方程^[10]可以得到

$ \left\{ {\begin{array}{*{20}{l}}{\Delta \sigma - 2k\Delta \varphi = 0, {\rm{点在}}\alpha {\rm{线上;}}}\\{\Delta \sigma + 2k\Delta \varphi = 0, {\rm{点在}}\beta {\rm{线上}}{\rm{.}}}\end{array}} \right. $

(2)

其中: Δσ和Δφ分别表示两点间σ和φ的变化量，σ表示静水压力，φ表示经过该点α线的切向与x轴的夹角(逆时针为正); k是材料的屈服切应力。对于管壁AB段和BC段分别列出σ和p的关系为

$\left\{ {\begin{array}{*{20}{l}}{\sigma = - p + k, \;\;\;{\rm{点在}}AB{\rm{上;}}}\\{\sigma = - k, \;\;\;\;\;\;\;\;\;{\rm{点在}}BC{\rm{上}}{\rm{.}}}\end{array}} \right.$

(3)

联立式(2)和(3)，并根据压头形状可以计算出接触部分AB段的压力分布：

$p(x) = \left\{ {\begin{array}{*{20}{l}}{2k\left( {1 + \frac{\pi }{2} + \arcsin \frac{x}{r} - \gamma } \right), }\\{x \in \left[ { - r\cos {\theta _0}, 0} \right];}\\{2k\left( {1 + {\theta _0} - \gamma } \right), }\\{x \in \left[ { - {x_B}, - r\cos {\theta _0}} \right).}\end{array}} \right.$

(4)

其中：x是接触点的横坐标；γ为∠BCD的角度，单位为rad；θ₀为驱动轮齿齿顶角的一半，即图 2中BI与y轴夹角；r是圆角半径；x_B为隆起最高点(图 2中点B)的横坐标绝对值。通过将压力分布从-x_B积分到0便可以得到力的大小，积分结果如式(5)、(6)所示。

${F_{\rm{N}}} = \left\{ {\begin{array}{*{20}{l}}{2kl\left[ {\left( {1 + {\theta _0} - \gamma } \right){x_B} + r\left( {1 - \sin {\theta _0}} \right)} \right], }&{{x_B} > r\cos {\theta _0};}\\{2kl\left[ {\left( {1 + \frac{\pi }{2} - \gamma } \right){x_B} + r - {x_B}\arcsin \frac{{{x_B}}}{r} - \sqrt {{r^2} - x_B^2} } \right], }&{{x_B} \le r\cos {\theta _0};}\end{array}} \right.$

(5)

${F_{\rm{T}}} = \left\{ {\begin{array}{*{20}{l}}{2kl\left[ {\left( {1 + {\theta _0} - \gamma } \right)\frac{{{x_B}\cos {\theta _0} - r}}{{\sin {\theta _0}}} + r\left( {1 + \frac{\pi }{2} - \gamma } \right) - r\cos {\theta _0}} \right], }&{{x_B} > r\cos {\theta _0};}\\{2kl\left[ {\left( {1 + \frac{\pi }{2} - \gamma } \right)\left( {r - \sqrt {{r^2} - x_B^2} } \right) + \sqrt {{r^2} - x_B^2} \arcsin \frac{{{x_B}}}{r} - {x_B}} \right], }&{{x_B} \le r\cos {\theta _0}.}\end{array}} \right.$

(6)

其中l为驱动轮齿厚。由于仅受正压力作用，驱动轮齿两侧压力值对称分布，则两侧切向力大小相等、方向相反，切向力合力为0；两侧法向力大小相等、方向相同，合力为式(5)计算值的2倍。
为了进一步计算，接下来推导压痕形状参量γ、x_B与点A纵坐标绝对值即压痕深度d的关系。根据之前的假设，针对理想刚塑形材料平面应变问题，利用体积守恒可以得到

${{S}_{\vartriangle BCD}}={{S}_{AIDE}}\triangleq S.$

(7)

不妨令DE=a，CE=c^*a，其中c^*为常数。根据Johnson给出的楔形压头滑移线解法^[13]，可以推导出c^*与θ₀的关系：

${{c}^{*}}\left( {{c}^{*}}-2 \right)\sin {{\theta }_{0}}=\sqrt{{{c}^{*2}}-{{\cos }^{2}}{{\theta }_{0}}}.$

(8)

若d=0，轮齿的轮廓线方程为f(x)，则可以将式(7)化为

$\left\{ \begin{array}{*{35}{l}} {{x}_{B}}=\frac{2S\tan {{\theta }_{0}}}{\left( {{c}^{*}}-1 \right)a}+a, \\ \gamma =\arctan \frac{2S}{{{\left( {{c}^{*}}-1 \right)}^{2}}{{a}^{2}}-2S\tan {{\theta }_{0}}}, \\ S=ad-\int_{0}^{a}{f}(x)\text{d}x. \\\end{array} \right.$

(9)

其中a与d有如下关系：

$ a=\left\{ \begin{array}{*{35}{l}} \sqrt{(2r-d)d}, & \ \ \ 0\le d<r-r\sin {{\theta }_{0}}; \\ (d-r)\tan {{\theta }_{0}}+\frac{r}{\cos {{\theta }_{0}}}, & \quad d\ge r-r\sin {{\theta }_{0}}. \\\end{array} \right. $

(10)

由此便建立了仅在正压力作用时，压痕形状(包括深度、隆起高度、隆起角度等)与力的大小的关系。
3 轮齿同时受到正压力和转矩作用的压痕计算模型当正压力和转矩同时作用于驱动轮时，根据式(1)可知，F_N＞0, F_T＞0，故轮齿不仅会形成y方向的压深d，还会形成x方向的偏移Δx。这样在计算左侧隆起时，式(9)中S变为图 3中阴影部分的面积，进而导致两侧隆起高度的不一致。新的压痕形状示意图如图 3中实线所示。

图 3 压痕形状示意图(Δx≤dtanθ0)

图选项

当保持正压力F不变，转矩M逐渐增大时，切向力F_T和偏移Δx会逐渐增大。但是在F_T增大到某个临界值F_Ta之后，右侧压力值会减小到塑性极限以下，滑移线理论将不再适用。此时的压痕形状将如图 4所示，轮齿将沿着右侧边进行滑移。当右侧边的压力值减小到0，对应切向力为F_Tb时，如果进一步增大转矩，则轮齿右侧边与管壁将不再接触，轮齿整体向左滑移即出现打滑现象。

图 4 压痕形状示意图(Δx=dtanθ0)

图选项

据此可将正压力和转矩同时作用的情况划分为3种情况：
1) 当0≤F_T≤F_Ta时，两侧边均接触且达到塑性极限，Δx≤dtanθ₀；
2) 当F_Ta＜F_T≤F_Tb时，两侧边均接触但左侧达到塑性极限，右侧未达到，Δx=dtanθ₀；
3) 当F_T＞F_Tb时，左侧接触且达到塑性极限，右侧不接触，Δx＞dtanθ₀。
由于第3种情况属于非正常工作状态，故此处不进一步讨论，以下将分别针对前2种情况分析压痕形状与作用力大小的关系。
当0≤F_T≤F_Ta, Δx≤dtanθ₀时，由于管壁左右两侧边均处于塑性极限状态，压力分布形式可用式(4)表示，但是区别在于两侧隆起的形状不同导致γ和x_B的不同，进而积分得到的力也存在差异。在轮齿形状和管壁材料已知的情况下，不妨将式(5)和(6)分别表示为F_N(γ, x_B)和F_T(γ, x_B)，则接触反力为

$\left\{ \begin{array}{*{35}{l}} {{F}_{\text{N}}}={{F}_{\text{N}1}}\left( {{\gamma }_{1}}, {{x}_{B1}} \right)+{{F}_{\text{N}2}}\left( {{\gamma }_{2}}, {{x}_{B2}} \right), \\ {{F}_{\text{T}}}={{F}_{\text{T}1}}\left( {{\gamma }_{1}}, {{x}_{B1}} \right)-{{F}_{\text{T}2}}\left( {{\gamma }_{2}}, {{x}_{B2}} \right). \\\end{array} \right.$

(11)

利用式(9)可以计算γ和x_B的数值，但是此时由于Δx不为0，因此计算面积的方程有所区别，新的方程如式(12)，其中i=1代表左侧，i=2代表右侧。D′E′长度a与d的关系仍为式(10)所示。

$\left\{ \begin{array}{*{35}{l}} {{x}_{Bi}}=\frac{2{{S}_{i}}\tan {{\theta }_{0}}}{\left( {{c}^{*}}-1 \right)a}+a, \quad i=1, 2; \\ {{\gamma }_{i}}=\arctan \frac{2{{S}_{i}}}{{{\left( {{c}^{*}}-1 \right)}^{2}}{{a}^{2}}-2{{S}_{i}}\tan {{\theta }_{0}}}, \\ {{S}_{i}}=\left[ a-{{(-1)}^{i}}\cdot \Delta x \right]d-\int_{{{(-1)}^{i}}\cdot \Delta x}^{a}{f}(x)\text{d}x. \\\end{array} \right.$

(12)

与仅受正压力情况不同，式(8)、(10)和(12)建立了(Δx, d)与(γ, x_B)的映射关系。由此也可以计算出临界状态Δx=dtanθ₀，且两侧边均处于塑性极限下的F_N和F_T，此时F_T即为F_Ta。
当F_Ta＜F_T≤F_Tb，Δx=dtanθ₀时，两侧边均接触但左侧达到塑性极限，右侧未达到，因此左侧边的压力分布仍可以利用滑移线理论进行求解。为了简化计算，假设两侧边的作用反力均沿该边法线方向(忽略圆角)，如图 4中F₁和F₂所示，则根据式(11)可得到

${{F}_{\text{N}1}}\left( {{\gamma }_{1}}, {{x}_{B1}} \right)=\frac{1}{2}\left( {{F}_{\text{N}}}+{{F}_{\text{T}}}\tan {{\theta }_{0}} \right).$

(13)

因为此时有Δx=dtanθ₀，所以由式(12)可得

$\left\{ \begin{array}{*{35}{l}} {{x}_{B1}}=\frac{2{{S}_{1}}\tan {{\theta }_{0}}}{\left( {{c}^{*}}-1 \right)a}+a, \\ {{\gamma }_{1}}=\arctan \frac{2{{S}_{1}}}{{{\left( {{c}^{*}}-1 \right)}^{2}}{{a}^{2}}-2{{S}_{1}}\tan {{\theta }_{0}}}, \\ {{S}_{1}}=\left( a+d\tan {{\theta }_{0}} \right)d-\int_{-d\tan {{\theta }_{0}}}^{a}{f}(x)\text{d}x. \\\end{array} \right.$

(14)

式(14)得到的γ₁和x_B1均为d的函数。联立式(5)、(8)、(10)、(13)和(14)，便可建立F_N1与d的关系。
由此便在仅受正压力的压痕模型基础上，推导出了正压力和转矩同时作用下的压痕形状与接触反力的关系。下面将通过仿真及试验对该模型进行分析和验证。
4 仿真结果为了验证提出的压痕理论模型，利用ABAQUS有限元分析软件对该问题进行了仿真计算。仿真采用显式(explicit)算法；针对接触区域利用了自适应网格划分技术(ALE)，用于解决大变形问题导致的网格畸变；轮齿采用解析刚体，齿顶角80°，圆角半径0.1 mm；管壁材料为铝合金，屈服强度为205.78 MPa，Young模量为71 000 MPa，达到屈服后屈服强度不变；网格划分采用平面应变、减缩积分四边形单元(CPE4R)。
为验证0≤F_T≤F_Ta和F_Ta＜F_T≤F_Tb这2种不同的情况，分别按式(15)和式(16)施加载荷仿真。

$\left\{ \begin{array}{l}{F_{\rm{N}}} = 8\;000t\;{\rm{且}}{F_{\rm{T}}} = 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;t \le 0.05;\\{F_{\rm{N}}} = 36\;000t - 1\;400{\rm{且}}{F_{\rm{T}}} = 2\;000t - 100, \;\;\;\;\;\;\;0.05 < t \le 0.15;\end{array} \right.$

(15)

$\left\{ \begin{array}{l}{F_{\rm{N}}} = 8\;000t\;{\rm{且}}{F_{\rm{T}}} = 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;t \le 0.05;\\{F_{\rm{N}}} = 36\;000t - 1\;400{\rm{且}}{F_{\rm{T}}} = 12\;000t - 600, \;\;\;\;\;\;\;0.05 < t \le 0.15.\end{array} \right.$

(16)

其中t为仿真时间，单位为s。当t=0.15 s时，式(15)对应的F_N=4 000 N，F_T=1 200 N；式(16)对应的F_N=4 000 N，F_T=1 200 N。此时，计算结果Mises应力云图如图 5所示，深色区域为达到屈服的区域，虚线表示加载前楔形压头的位置。可以看出，当切向力小于F_Ta时，与轮齿接触的两侧边均达到塑性极限，左右隆起高度有一定差别；当切向力大于F_Ta时，轮齿沿右侧边滑下，与轮齿接触的左侧达到塑性极限而右侧未达到，这与之前的分析一致，两侧隆起高度也呈现出明显的不同。

图 5 有限元仿真结果应力云图

图选项

从有限元仿真计算结果中获取d、Δx随t的变化数据，再与理论模型计算得到的d和Δx进行对比，得到的结果如图 6所示。可以看出，两种情况下理论模型均与仿真值能够较好地吻合，局部的波动是受到了网格划分和加载的影响。进一步分析可知，当法向力一定时，切向力是决定水平偏移量的关键，但也会对压入深度造成一定影响，这与式(13)相一致。

图 6 仿真模型与理论模型对比

图选项

提取t=150 ms时压痕形状和表面压力分布的仿真结果，与理论计算的结果进行对比，结果如图 7、8所示。针对图 7中的压痕形状，可以看出理论模型基本刻画了两侧隆起高度不一致的情况，并且较为精确地指出了接触区域的范围；但是由于理论模型中假设两侧隆起为直线，所以计算值与仿真结果之间存在一些偏差。

图 7 压痕形状对比

图选项

图 8 压力分布对比

图选项

对于接触部分的压力分布，从图 8可以看出，压力分布的规律与理论模型相同：在轮齿斜边处呈现均匀分布，在圆角处应力增大。并且理论计算值也与仿真结果比较一致。由于有限元网格划分和迭代计算接触的问题，局部呈现出压力值的波动；F_T＞F_Ta时为了计算压痕形状忽略了圆角部分，导致压力分布的峰值位置略有差异。
5 试验结果通过水平井爬行器在5.5寸套管井中的爬行试验^[27]，获得了切向力和法向力的测量数据如表 1所示。将测量数据的平均值代入理论模型，可以计算出理论压痕形状。
表 1 爬行试验测量数据

类别	1	2	3	4	平均值
法向力/N	2 925.3	2 350.1	3 220.1	3 677.4	3 043.2
切向力/N	731.3	666.4	652.9	633.3	664.7

表选项






爬行试验也获得了驱动轮在实际工作过程中留下的压痕痕迹，如图 9a所示。为了获取压痕形状的具体参数，试验利用了Nexview三维白光干涉形貌仪对其进行测量，如图 9b所示。该仪器利用非接触式白光扫描干涉原理可以精确测量物体表面的形貌，垂直分辨率可达0.1 nm。

图 9 压痕形貌测量试验

图选项

试验测得铝板表面5个不同压痕的形状截面，并将理论值与实际值进行了对比，结果如图 10所示。可以看出，理论模型较为准确地反映了压入深度和隆起高度，也解释了在转矩作用下两侧隆起高度出现的明显不一致现象。实际压痕宽度比理论值大的原因主要是驱动轮在实际滚动过程中既有正压状态也有斜压状态。若正压力和转矩均不变，则正压时正压力在法向投影最大，压入深度也达到最大；斜压时正压力分力与转矩效果类似，还用于产生水平偏移，导致实际压痕宽度变大。因此对于完整刻画实际压痕形状仍需要建立更为精确的斜压模型。

图 10 压痕形状对比

图选项

6 结论本文针对爬行器驱动轮同时受正压力和转矩作用的特点，利用滑移线理论建立了压痕形状与作用力关系的计算模型，发现：
1) 法向力一定时，切向力处于不同区间会导致塑性变形区域发生变化，进而导致压痕形状的不同和各处压力分布的改变。
2) 对于轮齿类似的楔形压头，法向力主要影响压入深度，切向力主要影响水平偏移，但也会使压入深度发生一些变化。
3) 该模型可用于计算爬行器驱动轮压入深度和压痕隆起高度，但是无法刻画整个滚动过程中压痕最终的宽度。
本文建立的正压模型可为爬行器牵引力建模与分析、驱动轮形状的优化设计、管壁材料选择和损伤分析等提供理论依据；同时也可以用于多个力作用下的楔形压头压痕分析等。由于该模型仅考虑了驱动轮正压情况，导致与实际压痕形状存在偏差。下一步需要在此模型基础上考虑斜压接触以及摩擦、材料参数等的影响，进一步完善接触模型。

参考文献

[1]	GIEM G, SHEIRETOV T, COUBLE Y. Wireline tractor for through-tubing intervention in wells with barefoot openhole completions[C]//SPE/ICoTA Coiled Tubing and Well Intervention Conference and Exhibition. Richardson, TX: Society of Petroleum Engineers, 2018, SPE-189923-MS.
[2]	FOUCHER P A, POORTEN R V. Wireline tractor advanced restriction navigation[C]//SPE/ICoTA Coiled Tubing and Well Intervention Conference and Exhibition. Richardson, TX: Society of Petroleum Engineers, 2018, SPE-189926-MS.
[3]	LI Y J, LIU Q Y, CHEN Y H, et al. Design and analysis of an active helical drive downhole tractor[J]. Chinese Journal of Mechanical Engineering, 2017, 30(2): 1-10.
[4]	TOKI T, BENYGZER M, WAHEDI K A, et al. Optimization of the operation time required for data gathering requirement in extended-reach drilling well by the use of open hole tractor in the United Arab Emirates[C]//Abu Dhabi International Petroleum Exhibition & Conference. Richardson, TX: Society of Petroleum Engineers, 2017, SPE-188680-MS.
[5]	刘振.水平井牵引器机构分析和建模仿真技术研究[D].哈尔滨: 哈尔滨工业大学, 2008. LIU Z. Research on mechanism analysis and modeling simulation of horizontal well tractor[D]. Harbin: Harbin Institute of Technology, 2008. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10213-2009229678.htm
[6]	曾华军.水平井牵引器驱动系统关键技术研究[D].哈尔滨: 哈尔滨工业大学, 2010. ZENG H J. Research on key technology of drive system of horizontal well tractor[D]. Harbin: Harbin Institute of Technology, 2010. (in Chinese) http://cdmd.cnki.com.cn/article/cdmd-10213-1011261813.htm
[7]	于淼.水平井牵引器拖动力数值分析[D].大庆: 东北石油大学, 2011. YU M. The mule downhole mumerical analysis of drag force on[D]. Daqing: Northeast Petroleum University, 2011. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10220-1011210254.htm
[8]	侯松伟.水平井牵引器运动学仿真与应用研究[D].东北石油大学, 2017. HOU S W. Research on kinematics simulation and application for horizontal well tractor[D]. Daqing: Northeast Petroleum University, 2017. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10220-1017096828.htm
[9]	张勇.水平井电缆牵引器的机械设计研究[D].北京: 中国石油大学, 2016. ZHANG Y. Research on mechanical design of cable tractor in horizontal wells[D]. Beijing: China University of Petroleum, 2016. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-11414-1018700630.htm
[10]	刘清友, 李雨佳, 任涛, 等. 水平井爬行器驱动轮力学分析[J]. 钻采工艺, 2014, 37(1): 68-71. LIU Q Y, LI Y J, REN T, et al. Mechanical analysis of driving wheels of horizontal wells tractor[J]. Drilling & Production Technology, 2014, 37(1): 68-71. DOI:10.3969/J.ISSN.1006-768X.2014.01.20 (in Chinese)
[11]	秦浩.水平井爬行器驱动机构结构优化[D].成都: 西南石油大学, 2017. CQIN H, Optimization of horizontal crawler's driving mechanism[D]. Chengdu: Southwest Petroleum University, 2017. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10615-1017108470.htm
[12]	希尔. 塑性数学理论[M]. 北京: 世界图书出版公司, 2013. HILL R. The mathematical theory of plasticity[M]. Beijing: World Book Inc, 2013. (in Chinese)
[13]	JOHNSON K L.接触力学[M].徐秉业, 罗学富, 刘信声, 宋国华, 孙学伟, 译.北京: 高等教育出版社, 1992. JOHNSON K L. Contact mechanics[M]. XU B Y, LUO X F, LIU X S, SONG G H, SUN X W, Trans. Beijing: Higher Education Press, 1992. (in Chinese)
[14]	GRUNZWEIG J, LONGMAN I M, PETCH N J. Calculations and measurements on wedge-indentation[J]. Journal of the Mechanics & Physics of Solids, 1954, 2(2): 81-86.
[15]	RICE J R. Plane strain slip line theory for anisotropic rigid/plastic materials[J]. Journal of the Mechanics & Physics of Solids, 1973, 21(2): 63-74.
[16]	BAI Y, DODD B. A slip-line field solution for plane-strain indentation by an obtuse-angle wedge[J]. International Journal of Mechanical Sciences, 1982, 24(2): 119-123. DOI:10.1016/0020-7403(82)90042-X
[17]	KOMVOPOULOS K, KOGUT L. Analysis of the spherical indentation cycle for elastic-perfectly plastic solids[J]. Journal of Materials Research, 2004, 19(12): 3641-3653. DOI:10.1557/JMR.2004.0468
[18]	TALJAT B, PHARR G M. Development of pile-up during spherical indentation of elastic-plastic solids[J]. International Journal of Solids & Structures, 2004, 41(14): 3891-3904.
[19]	HERNOT X, BARTIER O, BEKOUCHE Y, et al. Influence of penetration depth and mechanical properties on contact radius determination for spherical indentation[J]. International Journal of Solids & Structures, 2015, 43(14): 4136-4153.
[20]	HERNOT X, BARTIER O, MAUVOISIN G, et al. A universal formulation for indentation whatever the indenter geometry[J]. Mechanics of Materials, 2015, 81: 101-109. DOI:10.1016/j.mechmat.2014.11.006
[21]	JACKSON R L, GHAEDNIA H, POPE S. A solution of rigid-perfectly plastic deep spherical indentation based on slip-line theory[J]. Tribology Letters, 2015, 58(3): 1-7.
[22]	BOMIDI, JOHN A R SADEGHI F. Three-dimensional finite element elastic-plastic model for subsurface initiated spalling in rolling contacts[J]. Journal of Tribology, 2013, 136(1): 011402. DOI:10.1115/1.4025841
[23]	WANG Z J, JIN X Q, KEER L M, et al. Numerical methods for contact between two joined quarter spaces and a rigid sphere[J]. International Journal of Solids & Structures, 2012, 49(18): 2515-2527.
[24]	WANG Z J, JIN X Q, LIU S B, et al. A new fast method for solving contact plasticity and its application in analyzing elasto-plastic partial slip[J]. Mechanics of Materials, 2013, 60(1): 18-35.
[25]	ZHOU Q, JIN X Q, WANG Z J, et al. An efficient approximate numerical method for modeling contact of materials with distributed inhomogeneities[J]. International Journal of Solids & Structures, 2014, 51(19-20): 3410-3421.
[26]	WANG Z J, YU H, WANG Q. Analytical solutions for elastic fields caused by eigenstrains in two joined and perfectly bonded half-spaces and related problems[J]. International Journal of Plasticity, 2016, 76: 1-28. DOI:10.1016/j.ijplas.2015.07.006
[27]	孙可平, 杨东超, 常旭, 等. 水平井电缆牵引器滚轮与管壁压力的测量[J]. 石油学报, 2018, 39(6): 712-718. SUN K P, YANG D C, CHANG X, et al. Measurement of the normal contact force between the cable tractor's wheels and the pipe wall of horizontal well[J]. Acta Petrolei Sinica, 2018, 39(6): 712-718. (in Chinese)