轮式爬行器^[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。

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

$ \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为接触部分的压力分布。

根据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中实线所示。

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

据此可将正压力和转矩同时作用的情况划分为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时，轮齿沿右侧边滑下，与轮齿接触的左侧达到塑性极限而右侧未达到，这与之前的分析一致，两侧隆起高度也呈现出明显的不同。

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

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

对于接触部分的压力分布，从图 8可以看出，压力分布的规律与理论模型相同：在轮齿斜边处呈现均匀分布，在圆角处应力增大。并且理论计算值也与仿真结果比较一致。由于有限元网格划分和迭代计算接触的问题，局部呈现出压力值的波动；F_T＞F_Ta时为了计算压痕形状忽略了圆角部分，导致压力分布的峰值位置略有差异。

5 试验结果

通过水平井爬行器在5.5寸套管井中的爬行试验^[27]，获得了切向力和法向力的测量数据如表 1所示。将测量数据的平均值代入理论模型，可以计算出理论压痕形状。

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

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

6 结论

本文针对爬行器驱动轮同时受正压力和转矩作用的特点，利用滑移线理论建立了压痕形状与作用力关系的计算模型，发现：

1) 法向力一定时，切向力处于不同区间会导致塑性变形区域发生变化，进而导致压痕形状的不同和各处压力分布的改变。

2) 对于轮齿类似的楔形压头，法向力主要影响压入深度，切向力主要影响水平偏移，但也会使压入深度发生一些变化。

3) 该模型可用于计算爬行器驱动轮压入深度和压痕隆起高度，但是无法刻画整个滚动过程中压痕最终的宽度。

本文建立的正压模型可为爬行器牵引力建模与分析、驱动轮形状的优化设计、管壁材料选择和损伤分析等提供理论依据；同时也可以用于多个力作用下的楔形压头压痕分析等。由于该模型仅考虑了驱动轮正压情况，导致与实际压痕形状存在偏差。下一步需要在此模型基础上考虑斜压接触以及摩擦、材料参数等的影响，进一步完善接触模型。