孙可平, 杨东超, 常旭, 朱衡, 鲁沛昕, 陈恳
清华大学 机械工程系, 北京 100084

收稿日期：2019-04-09
作者简介：孙可平(1979-), 男, 博士研究生
通信作者：杨东超, 副研究员, E-mail:ydc@tsinghua.edu.cn

摘要：轮式水平井爬行器驱动轮与套管管壁接触过程中会在管壁上留下塑性变形压痕，该压痕形貌对爬行器的运动性能影响较大，而轮齿与管壁材料的斜压接触过程非常复杂，目前还没有相应的理论模型来分析该过程。该文详细分析了正压力和转矩同时作用下，轮齿斜压管壁时不同载荷条件下压痕形成的规律，将载荷分为6个不同的作用区域。利用滑移线理论，分析计算了每个区域轮齿不同接触状态下的压力分布，建立了轮齿载荷与压痕形貌之间的力学计算模型。在此基础上，分析和计算了各作用区域间的临界载荷，并介绍了理论模型的正向和逆向应用方法。最后对4个稳定作用区域的理论模型从压深、滑移量、压痕形貌和接触面上压力分布等4个方面进行了仿真验证，并进行了压痕验证试验。仿真和试验结果均表明：对轮齿载荷作用的分区及所建立的理论计算模型是合理的。为今后爬行器载荷控制策略的优化以及驱动轮多齿牵引力的建模与优化奠定了基础。
关键词：爬行器压痕形貌塑性接触模型滑移线理论压力分布
Analysis of oblique indentations caused by the tractor driving wheel and the casing pipe wall
SUN Keping, YANG Dongchao, CHANG Xu, ZHU Heng, LU Peixin, CHEN Ken
Department of Mechanical Engineering, Tsinghua University, Beijing 100084, China


Abstract: Contact between the driving wheel of a horizontal well crawler and the casing pipe wall causes plastic deformation indentations in the pipe wall. The indentation morphology then greatly influences the crawler motion with the oblique indentation contact between the wheel teeth and the pipe wall being a very complex process. A theoretical model was developed to analyze the indentation formation for various loading conditions with positive pressure and torque on the gear tooth to the oblique indentation. The loading effects can be divided into 6 different action areas. Slip line theory was used to analyze the pressure distribution for various contact states of the gear tooth in each region in a mechanical calculation model relating the load on the gear tooth and the indentation morphology. The critical loads between the action zones are identified with forward and backward application methods of the theoretical model. Finally, theoretical models of the four stabilization zones are characterized based on their indentation depth, amount of slip, indentation morphology and pressure distribution on the contact surface, and the indentation verification test was also carried out. The simulations and tests show that the load division between the gear teeth is reasonable and the theoretical model is accurate. This lays a foundation for optimization of the crawler load control strategy and the modeling and optimization of the multi-tooth traction force on the drive wheel.
Key words: crawlerindentation morphologyplastic contact modelslip line theorypressure distribution
随着水平井的广泛应用^[1-2]，测井时需要使用爬行器来向井下输送仪器或设备^[3-10]。轮式爬行器一般由电机或液压装置通过曲柄滑块机构驱动支撑臂张开实现驱动轮与管壁接触并产生法向接触力，同时由潜油电机为驱动轮提供转矩，实现爬行器在水平井内沿套管轴向的前进运动。由于爬行器工作在井下油、水、泥混合外加高温高压的恶劣环境下，牵引力大小直接影响其负载能力、井下爬行距离、爬行速度、越障能力、井下脱困能力等。轮式爬行器主要依靠驱动轮与管壁的接触产生牵引力，由于正压力和牵引力均较大，驱动轮会在管壁上留下塑性变形压痕；而爬行器在井下工作时常会因控制不当出现驱动轮打滑或者电机堵转现象。测试实践发现，爬行器在套管内不同的载荷和运动状态下所产生的塑性变形压痕在形貌上有较大差别，故研究驱动轮轮齿与管壁之间的接触对于优化控制策略、提高牵引效率和爬行器的综合性能具有重要的意义。目前针对该类型接触的研究工作较少，基本上都是利用摩擦因数的测试经验值来处理正压力与牵引力的关系^[11-15]。文[16-17]主要利用有限元进行分析，但是没有建立完整的理论模型，无法进一步对滚动接触进行分析和对牵引力进行计算。文[18-19]认为：当材料的塑性应变非常大，以至于弹性应变相对微乎其微时，就可以看作符合理想刚塑性假设。在受载状态下，发生塑性流动区域的应力状态可以利用滑移线场来表示。而滑移线理论很早就被用来求解刚性压头压入问题，Grunzweig等^[20]在文[18]无摩擦模型的基础上考虑了不同摩擦因数的影响，建立了平均压力与楔形压头顶角的关系，并与实验结果进行了对比。Rice^[21]对平板压头的压入问题进行了计算，并考虑了材料各向异性的影响。BAI等^[22]考虑了滑移场中可能出现的应力不连续的现象，重新计算了钝角楔形压头的压入问题。文[23-24]利用滑移线理论得到了球形压头的平均压力方程，并分析了平面应变下理想刚塑性圆柱的压痕问题。但是以上模型均只考虑了法向力垂直正压的情况，没有考虑法向力和切向力同时作用的情况，且忽略了楔形压头圆角对应力和塑性变形的影响。常旭等^[25]利用滑移线理论分析了爬行器驱动轮正压套管管壁时，法向力和切向力同时作用下管壁压痕与载荷之间的关系。由于爬行器驱动轮运动过程中正压只是其中的一个特殊状态，绝大多数时刻均处于倾斜压入状态，因此有必要研究更普遍的法向力和切向力同时作用下的倾斜压入问题。
本文利用滑移线理论分析和推导正压力和转矩同时作用下，轮齿斜压管壁时压痕形貌与载荷之间的关系，并利用有限元仿真和试验对该理论模型进行验证，以期为爬行器驱动轮多齿牵引力的分析与建模以及爬行器控制策略的优化提供依据。
1 载荷作用分区1.1 载荷分析模型爬行器驱动轮在水平井套管内滚动时，接触轮齿的中轴线与竖直方向的夹角在不断变化，绝大多数时刻轮齿与管壁处于斜压状态，只有一个瞬时轮齿处于垂直正压状态。为提高所研究问题的普适性，本文选取轮齿以倾斜角度θ压入管壁时的状态进行分析。爬行器驱动轮在正压力F和转矩M的作用下，其力学模型如图 1所示，将其作用力分解为切向力F_T和法向力F_N，根据力平衡原理可得：

图 1 驱动轮力学模型

图选项

$\left\{ \begin{array}{l}{F_{\rm{N}}} = F,\\{F_{\rm{T}}} = \frac{M}{{R\cos \theta }} - F \cdot \tan \theta .\end{array} \right.$

(1)

其中R表示驱动轮半径。而轮齿在压入管壁的过程中，轮齿尖部既有垂直于套管轴线方向的位移，即压入深度d, 同时还有沿套管轴线方向的水平位移Δx，本文将研究建立轮齿所受的切向力F_T、法向力F_N与压深d、滑移量Δx之间的函数关系。
1.2 作用区域分界由式(1)可以看出，驱动轮所受到的F和M均可等效为F_T和F_N。轮齿受力后将引起管壁接触区域的塑性变形。轮齿的变形可忽略不计，但其运动轨迹将受到F_T和F_N的影响，进而决定了管壁接触区域的塑性变形。根据F_T和F_N之间的关系可以将斜压时轮齿的载荷分为6个作用区域，如图 2所示。其中，γ₁和γ₂区域为轮齿左右两侧的管壁材料均有塑性变形，γ₁区域水平滑移量Δx为正(图 1中向左)，γ₂区域水平滑移量Δx为负(图 1中向右)；β区域为轮齿只有左侧的管壁材料有塑性变形，轮齿沿右侧面向左下方滑动形成压痕；ω区域为轮齿只有右侧的管壁材料有塑性变形，轮齿沿左侧面向右下方滑动形成压痕；α和φ区域为2个不稳定区域，α区域由于正向(图 1中向左)的切向力F_T过大，轮齿将沿左侧面滑出，φ区域由于负向的切向力F_T过大，轮齿将沿右侧面滑出，各区域塑性变形仿真效果如图 3所示。

图 2 FT作用区域分界图

图选项

图 3 各区域塑性变形仿真效果图

图选项

2 斜压力学模型驱动轮轮齿与管壁的弹塑性变形问题涉及到接触非线性和材料非线性2个难点，目前尚无有效的理论模型可以精确刻画其变化规律。为了简化分析，进行如下假设：
假设1??由于驱动轮的材料为20CrMoTi，屈服强度(885 MPa)和硬度(217 HB)均大幅高于试验用6061铝板，故假设轮齿在接触过程中不发生变形。
假设2??由于轮齿厚度为10 mm，压痕深度为0.5 mm左右，二者相差较大，故假设在齿厚方向压力分布不发生变化，故转化为二维平面应变问题。
假设3??经过前期仿真表明，达到实际压深时弹性变形部分占比较小(2%~5%)，故假设管壁材料为理想刚塑性，采用Mises屈服准则下的Prandtl-Reuss本构关系和关联塑性流动法则进行分析。
假设4??由于摩擦力影响因素较多(各点接触状态、轮齿表面粗糙度、接触面温度等)，此处暂不考虑接触部分摩擦力的影响。
由此便可利用滑移线理论来分析处理该问题。
2.1 γ₁、γ₂区域力学建模当轮齿载荷位于γ₁或γ₂区域时，轮齿左右两侧的管壁材料均有塑性变形，建立压痕附近的滑移线场如图 4所示，其中p为接触部分的压力分布。根据Hencky方程^[19]得到：

图 4 轮齿两侧均有塑性变形时的力学分析模型

图选项

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

(2)

其中：Δ表示两点间的变化量；σ表示静水压力；φ表示经过该点的α₀线的切向与x轴的夹角，逆时针为正；k是材料的屈服切应力。则轮齿接触部分压力分布可表示为：

$p\left( x \right) = \left\{ \begin{array}{l}{p_1}(x) = 2k\left( {1 + \left( {{\theta _0} + \theta } \right) - {\lambda _1}} \right),\\\;\;\;x \in \left[ { - {x_{{B_1}}}, - r\cos \left( {{\theta _0} + \theta } \right)} \right];\\{p_2}(x) = 2k\left( {1 + \frac{{\rm{ \mathsf{ π} }}}{2} + \arcsin \frac{x}{r} - {\lambda _1}} \right),\\\;\;\;x \in \left[ { - r\cos \left( {{\theta _0} + \theta } \right),0} \right];\\{p_3}(x) = 2k\left( {1 + \frac{{\rm{ \mathsf{ π} }}}{2} - \arcsin \frac{x}{r} - {\lambda _2}} \right),\\\;\;\;\;x \in \left[ {0,r\cos \left( {{\theta _0} - \theta } \right)} \right];\\{p_4}(x) = 2k\left( {1 + \left( {{\theta _0} - \theta } \right) - {\lambda _2}} \right),\\\;\;\;\;x \in \left[ {r\cos \left( {{\theta _0} - \theta } \right),{x_{{B_2}}}} \right].\end{array} \right.$

(3)

其中：x表示接触点的横坐标，x_B1表示左侧隆起最高点(图 4中B₁点)的横坐标绝对值，x_B2表示右侧隆起最高点(图 4中B₂点)的横坐标绝对值，λ₁表示∠B₁C₁D₁的角度值(rad)，λ₂表示∠B₂C₂D₂的角度值(rad)，θ₀表示驱动轮齿顶角的一半，θ表示压入角度，r表示圆角半径。
通过将压力分布函数沿轮齿接触部分分别按式(4)—(5)进行积分可得到轮齿法向力F_N和切向力F_T，其中，b表示驱动轮齿厚。

${F_{\rm{N}}} = {f_1}\left( {{x_{{B_1}}},{x_{{B_2}}},{\lambda _1},{\lambda _2}} \right) =\\ \left\{ \begin{array}{l}\left\{ \begin{array}{l}\int_{ - {x_{{B_1}}}}^0 {{p_2}\left( x \right) \cdot b \cdot {\rm{d}}x} + \int_0^{{x_{{B_2}}}} {{p_3}\left( x \right) \cdot b \cdot {\rm{d}}x} ,\\\;\;\;\;0 < {x_{{B_1}}} \le r\cos \left( {{\theta _0} + \theta } \right),0 < {x_{{B_2}}} \le r\cos \left( {{\theta _0} - \theta } \right);\\\int_{ - {x_{{B_1}}}}^0 {{p_2}\left( x \right) \cdot b \cdot {\rm{d}}x} + \int_0^{r\cos \left( {{\theta _0} + \theta } \right)} {{p_3}\left( x \right) \cdot b \cdot {\rm{d}}x} + \int_{r\cos \left( {{\theta _0} + \theta } \right)}^{{x_{{B_2}}}} {{p_4}\left( x \right) \cdot b \cdot {\rm{d}}x} ,\\\;\;\;\;0 < {x_{{B_1}}} \le r\cos \left( {{\theta _0} + \theta } \right),{x_{{B_2}}} > r\cos \left( {{\theta _0} - \theta } \right);\end{array} \right.\\\left\{ \begin{array}{l}\int_{ - {x_{{B_1}}}}^{ - r\cos \left( {{\theta _0} + \theta } \right)} {{p_1}} \left( x \right) \cdot b \cdot {\rm{d}}x + \int_{ - r\cos \left( {{\theta _0} + \theta } \right)}^0 {{p_2}} (x) \cdot b \cdot {\rm{d}}x + \int_0^{{x_{{B_2}}}} {{p_3}} (x) \cdot b \cdot {\rm{d}}x,\\\;\;\;\;{x_{{B_1}}} > r\cos \left( {{\theta _0} + \theta } \right),0 < {x_{{B_2}}} \le r\cos \left( {{\theta _0} - \theta } \right);\\\int_{ - {x_{{B_1}}}}^{ - r\cos \left( {{\theta _0} + \theta } \right)} {{p_1}} (x) \cdot b \cdot {\rm{d}}x + \int_{ - r\cos \left( {{\theta _0} + \theta } \right)}^0 {{p_2}} (x) \cdot b \cdot {\rm{d}}x \cdots \\ \cdots + \int_0^{r\cos \left( {{\theta _0} - \theta } \right)} {{p_3}} (x) \cdot b \cdot {\rm{d}}x + \int_{\left. {r\cos {\theta _0} - \theta } \right)}^{{x_{{B_2}}}} {{p_4}} (x) \cdot b \cdot {\rm{d}}x,\\\;\;\;\;{x_{{B_1}}} > r\cos \left( {{\theta _0} + \theta } \right),{x_{{B_2}}} > r\cos \left( {{\theta _0} - \theta } \right).\end{array} \right.\end{array} \right.$

(4)

${F_{\rm{T}}} = {f_2}\left( {{x_{{B_1}}},{x_{{B_2}}},{\lambda _1},{\lambda _2}} \right) =\\ \left\{ \begin{array}{l}\left\{ \begin{array}{l}\int_{ - {x_{{B_1}}}}^0 {{p_2}} (x) \cdot b \cdot \frac{{ - x}}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x - \int_0^{{x_{{B_2}}}} {{p_3}} (x) \cdot b \cdot \frac{x}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x,\\\;\;\;\;0 < {x_{{B_1}}} \le r\cos \left( {{\theta _0} + \theta } \right),0 < {x_{{B_2}}} \le r\cos \left( {{\theta _0} - \theta } \right);\\\int_{ - {x_{{B_1}}}}^0 {{p_2}} (x) \cdot b \cdot \frac{{ - x}}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x - \int_0^{r\cos \left( {{\theta _0} - \theta } \right)} {{p_3}} (x) \cdot b \cdot \frac{x}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x - \\\;\;\;\;\int_{r\cos \left( {{\theta _0} - \theta } \right)}^{{x_{{B_2}}}} {{p_4}} (x) \cdot b \cdot \frac{{\cos \left( {{\theta _0} - \theta } \right)}}{{\sin \left( {{\theta _0} - \theta } \right)}} \cdot {\rm{d}}x,\\\;\;\;\;0 < {x_{{B_1}}} \le r\cos \left( {{\theta _0} + \theta } \right),{x_{{B_2}}} > r\cos \left( {{\theta _0} - \theta } \right);\end{array} \right.\\\left\{ \begin{array}{l}\int_{ - {x_{{B_1}}}}^{ - r\cos \left( {{\theta _n} + \theta } \right)} {{p_1}} (x) \cdot b \cdot \frac{{\cos \left( {{\theta _0} + \theta } \right)}}{{\sin \left( {{\theta _0} + \theta } \right)}} \cdot {\rm{d}}x + \int_{ - r\cos \left( {{\theta _0} + \theta } \right)}^0 {{p_2}} (x) \cdot b \cdot \frac{{ - x}}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x - \\\;\;\;\;\;\int_0^{{x_{{B_2}}}} {{p_3}} (x) \cdot b \cdot \frac{x}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x,\\{x_{{B_1}}} > r\cos \left( {{\theta _0} + \theta } \right),0 < {x_{{B_2}}} \le r\cos \left( {{\theta _0} - \theta } \right);\\\int_{ - {x_{{B_1}}}}^{ - r\cos \left( {{\theta _0} + \theta } \right)} {{p_1}} (x) \cdot b \cdot \frac{{\cos \left( {{\theta _0} + \theta } \right)}}{{\sin \left( {{\theta _0} + \theta } \right)}} \cdot {\rm{d}}x + \int_{ - r\cos \left( {{\theta _0} + \theta } \right)}^0 {{p_2}} (x) \cdot b \cdot \frac{{ - x}}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x \cdots \\\;\;\;\;\; \cdots - \int_0^{r\cos \left( {{\theta _0} + \theta } \right)} {{p_3}} (x) \cdot b \cdot \frac{x}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x - \int_{r\cos \left( {{\theta _0} + \theta } \right)}^{{x_{{B_2}}}} {{p_4}} (x) \cdot b \cdot \frac{{\cos \left( {{\theta _0} - \theta } \right)}}{{\sin \left( {{\theta _0} - \theta } \right)}} \cdot {\rm{d}}x,\\\;\;\;\;\;{x_{{B_1}}} > r\cos \left( {{\theta _0} + \theta } \right),{x_{{B_2}}} > r\cos \left( {{\theta _0} - \theta } \right).\end{array} \right.\end{array} \right.$

(5)

2.2 α区域力学建模当轮齿载荷位于α区域时，只有轮齿左侧的管壁材料有塑性变形，而右侧面材料不受力，建立压痕附近的滑移线场如图 5所示，其接触部分压力分布可表示为：

图 5 轮齿左侧有塑变时的力学分析模型

图选项

$p\left( x \right) = \left\{ \begin{array}{l}{p_1}\left( x \right) = 2k\left( {1 + \left( {{\theta _0} + \theta } \right) - \lambda } \right),\;\;\;\;\;\;\;\;\;\;\;\;\;\;x \in \left[ { - {x_{{B_1}}}, - r\cos \left( {{\theta _0} - \theta } \right)} \right];\\{p_2}\left( x \right) = 2k\left( {1 + \frac{{\rm{ \mathsf{ π} }}}{2} + \arcsin \frac{x}{r} - {\lambda _1}} \right),\;\;\;\;\;\;x \in \left[ { - r\cos \left( {{\theta _0} - \theta } \right),0} \right];\\{p_3}\left( x \right) = 2k\left( {1 + \frac{{\rm{ \mathsf{ π} }}}{2} + \arcsin \frac{x}{r} - {\lambda _1}} \right),\;\;\;\;\;\;x \in \left[ {0,r\cos \left( {{\theta _0} - \theta } \right)} \right].\end{array} \right.$

(6)

通过将压力分布函数沿轮齿接触部分分别按式(7)—(8)进行积分可得到轮齿法向力F_N和切向力F_T。

$\begin{array}{*{20}{c}}{{F_{\rm{N}}} = {f_{{\rm{1L}}}}\left( {{x_{{B_1}}},{\lambda _1}} \right) = }\\{\left\{ \begin{array}{l}\int_{ - {x_{{B_1}}}}^0 {{p_2}} (x) \cdot b \cdot {\rm{d}}x + \int_0^{r\cos \left( {{\theta _0} - \theta } \right)} {{p_3}} (x) \cdot b \cdot {\rm{d}}x,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0 < {x_{{B_1}}} \le r\cos \left( {{\theta _0} + \theta } \right);\\\int_{ - {x_{{B_1}}}}^{ - r\cos \left( {{\theta _0} + \theta } \right)} {{p_1}} (x) \cdot b \cdot {\rm{d}}x + \int_{ - r\cos \left( {{\theta _0} + \theta } \right)}^0 {{p_2}} (x) \cdot b \cdot {\rm{d}}x + \int_0^{r\cos \left( {{\theta _0} + \theta } \right)} {{p_3}} (x) \cdot b \cdot {\rm{d}}x,\;\;\;\;{x_{{B_1}}} > r\cos \left( {{\theta _0} + \theta } \right).\end{array} \right.}\end{array}$

(7)

$\begin{array}{*{20}{c}}{{F_{\rm{T}}} = {f_{2{\rm{L}}}}\left( {{x_{{B_1}}},{\lambda _1}} \right) = }\\{\left\{ \begin{array}{l}\int_{ - {x_{{B_1}}}}^0 {{p_2}} (x) \cdot b \cdot \frac{{ - x}}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x - \int_0^{r\cos \left( {{\theta _0} + \theta } \right)} {{p_3}} (x) \cdot b \cdot \frac{x}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0 < {x_{{B_1}}} \le r\cos \left( {{\theta _0} + \theta } \right);\\\int_{ - {x_{{B_1}}}}^{ - r\cos \left( {{\theta _0} + \theta } \right)} {{p_1}} (x) \cdot b \cdot \frac{{\cos \left( {{\theta _0} + \theta } \right)}}{{\sin \left( {{\theta _0} + \theta } \right)}} \cdot {\rm{d}}x + \int_{ - r\cos \left( {{\theta _0} + \theta } \right)}^0 {{p_2}} (x) \cdot b \cdot \frac{x}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x - \\\int_0^{r\cos \left( {{\theta _0} + \theta } \right)} {{p_3}} (x) \cdot b \cdot \frac{x}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{x_{{B_1}}} > r\cos \left( {{\theta _0} + \theta } \right).\end{array} \right.}\end{array}$

(8)

2.3 φ区域力学建模当轮齿载荷位于φ区域时，只有轮齿右侧的管壁材料有塑性变形，而左侧面材料不受力，建立压痕附近的滑移线场如图 6所示，其接触部分压力分布可表示为：

图 6 轮齿右侧有塑变时的力学分析模型

图选项

$p\left( x \right) = \left\{ \begin{array}{l}{p_2}(x) = 2k\left( {1 + \frac{{\rm{ \mathsf{ π} }}}{2} + \arcsin \frac{x}{r}} \right),\;\;\;\;\;\;\;\;\;\;\;x \in \left[ { - r\cos \left( {{\theta _0} + \theta } \right),0} \right];\\{p_3}(x) = 2k\left( {1 + \frac{{\rm{ \mathsf{ π} }}}{2} - \arcsin \frac{x}{r} - {\lambda _2}} \right),\;\;\;\;\;x \in \left[ {0,r\cos \left( {{\theta _0} - \theta } \right)} \right];\\{p_4}(x) = 2k\left( {1 + \left( {{\theta _0} - \theta } \right) - {\lambda _2}} \right),\;\;\;\;\;\;\;\;\;\;\;\;x \in \left[ {r\cos \left( {{\theta _0} - \theta } \right),{x_{{B_2}}}} \right].\end{array} \right.$

(9)

通过将压力分布函数沿轮齿接触部分分别按式(10)—(11)进行积分可得到轮齿法向力F_N和切向力F_T。

$\begin{array}{*{20}{c}}{{F_{\rm{N}}} = {f_{1{\rm{R}}}}\left( {{x_{{B_2}}},{\lambda _2}} \right) = }\\{\left\{ \begin{array}{l}\int_{ - r\cos \left( {{\theta _0} + \theta } \right)}^0 {{p_2}} (x) \cdot b \cdot {\rm{d}}x + \int_0^{{x_{{B_2}}}} {{p_3}} (x) \cdot b \cdot {\rm{d}}x,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0 < {x_{{B_2}}} \le r\cos \left( {{\theta _0} - \theta } \right);\\\int_{ - r\cos \left( {{\theta _1} + \theta } \right)}^0 {{p_2}} (x) \cdot b \cdot {\rm{d}}x + \int_0^{r\cos \left( {{\theta _0} - \theta } \right)} {{p_3}} (x) \cdot b \cdot {\rm{d}}x + \int_{r\cos \left( {{\theta _0} - \theta } \right)}^{{x_{{B_s}}}} {{p_4}} (x) \cdot b \cdot {\rm{d}}x,\;\;\;\;\;\;{x_{{B_2}}} > r\cos \left( {{\theta _0} - \theta } \right).\end{array} \right.}\end{array}$

(10)

$\begin{array}{*{20}{c}}{{F_{\rm{T}}} = {f_{2{\rm{R}}}}\left( {{x_{{B_2}}},{\lambda _2}} \right) = }\\{\left\{ \begin{array}{l}\int_{ - r\cos \left( {{\theta _0} + \theta } \right)}^0 {{p_2}} (x) \cdot b \cdot \frac{{ - x}}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x - \int_0^{{x_{{B_2}}}} {{p_3}} (x) \cdot b \cdot \frac{x}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x,\;\;\;\;\;\;\;\;\;\;\;0 < {x_{{B_2}}} \le r\cos \left( {{\theta _0} - \theta } \right);\\\int_{ - r\cos \left( {{\theta _0} + \theta } \right)}^0 {{p_2}} (x) \cdot b \cdot \frac{{ - x}}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x - \int_0^{r\cos \left( {{\theta _0} - \theta } \right)} {{p_3}} (x) \cdot b \cdot \frac{x}{{\sqrt {{r^2} - {x^2}} }} \cdot {\rm{d}}x - \\\int_{r\cos \left( {{\theta _0} - \theta } \right)}^{{x_{{B_2}}}} {{p_4}} (x) \cdot b \cdot \frac{{\cos \left( {{\theta _0} - \theta } \right)}}{{\sin \left( {{\theta _0} - \theta } \right)}} \cdot {\rm{d}}x,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{x_{{B_2}}} > r\cos \left( {{\theta _0} - \theta } \right).\end{array} \right.}\end{array}$

(11)

2.4 β区域力学建模当轮齿载荷位于β区域时，只有轮齿左侧的管壁材料有塑性变形，右侧材料虽没有塑性变形，但是承受压力，建立压痕附近的滑移线场如图 5所示，其中F_R为轮齿右侧面直线段部分所承受的总压力，方向垂直于右侧接触面指离轮齿。接触部分压力分布同样可表示为式(6)。
由于此时轮齿右侧面材料也承受了一部分载荷，故将压力分布函数沿轮齿接触部分积分得到的法向力和切向力不能单独与载荷形成平衡，而是与F_R一起与轮齿载荷形成平衡。记式(7)的积分结果为F_NLD，记式(8)的积分结果为F_TLD，则在材料塑性变形过程中可建立如下平衡方程组：

$\left\{ {\begin{array}{*{20}{l}}{{F_{{\rm{NLD}}}} + {F_{\rm{R}}} \cdot \sin \left( {{\theta _0} - \theta } \right) = {F_{\rm{N}}},}\\{{F_{{\rm{TLD}}}} - {F_{\rm{R}}} \cdot \cos \left( {{\theta _0} - \theta } \right) = {F_{\rm{T}}}.}\end{array}} \right.$

(12)

由此可推导得到β区域的力学模型如下：

${F_{{\rm{NLD}}}} + \tan \left( {{\theta _0} - \theta } \right) \cdot {F_{{\rm{TLD}}}} = {F_{\rm{N}}} + {F_{\rm{T}}} \cdot \tan \left( {{\theta _0} - \theta } \right).$

(13)

2.5 ω区域力学建模当轮齿载荷位于ω区域时，只有轮齿右侧的管壁材料有塑性变形，左侧材料虽没有塑性变形，但是承受压力，建立压痕附近的滑移线场如图 6所示，其中F_L为轮齿左侧面直线段部分所承受的总压力，方向垂直于左侧接触面指离轮齿。接触部分压力分布同样可表示为式(9)。由于此时轮齿左侧面材料也承受了一部分载荷，故将压力分布函数沿轮齿接触部分积分得到的法向力和切向力同样不能单独与载荷形成平衡，而是与F_L一起与轮齿载荷形成平衡。记式(10)的积分结果为F_NRD，记式(11)的积分结果为F_TRD，则在材料塑性变形过程中可建立如下平衡方程组：

$\left\{ {\begin{array}{*{20}{l}}{{F_{{\rm{NRD}}}} + {F_{\rm{L}}} \cdot \sin \left( {{\theta _0} + \theta } \right) = {F_{\rm{N}}},}\\{{F_{{\rm{TRD}}}} + {F_{\rm{L}}} \cdot \cos \left( {{\theta _0} + \theta } \right) = {F_{\rm{T}}}.}\end{array}} \right.$

(14)

由此可推导得到ω区域的力学模型如下：

${F_{{\rm{NRD}}}} - \tan \left( {{\theta _0} + \theta } \right) \cdot {F_{{\rm{TRD}}}} = {F_{\rm{N}}} - {F_{\rm{T}}} \cdot \tan \left( {{\theta _0} + \theta } \right).$

(15)

综上所述，建立了轮齿载荷F_N和F_T与压痕隆起位置参量x_B1、x_B2和压痕隆起高度参量λ₁、λ₂之间的函数关系。
3 压痕形貌和压深与滑移量的关系3.1 隆起位置和高度与压深和滑移量的关系为了进一步计算，希望找出压痕形貌参量x_B1、x_B2、λ₁、λ₂与轮齿压入深度d、滑移量Δx之间的函数关系。如图 7所示，记压深d=0时轮齿的轮廓线方程为f(x)，取D₁E₃=a₁，C₁E₃=c₁·a₁，D₂E₃=a₂，C₂E₃=c₂·a₂根据材料塑性变形过程中体积不变的定理有：S_ΔB1C1D1=S_D1FE1E=S₁和S_ΔB2C2D2=S_D2GE1E=S₂。经推导得到压痕形貌参量x_B1、x_B2、λ₁、λ₂与压深d、滑移量Δx的函数关系如下：

图 7 斜压时压痕形貌图

图选项

$\left\{ {\begin{array}{*{20}{l}}{{S_1} = \left( {{a_1} + \frac{{\Delta x}}{2}} \right)d - \int_{ - {a_1}}^0 f (x) \cdot {\rm{d}}x,}\\{{x_{{B_1}}} = \frac{{2{S_1} \cdot \tan \left( {{\theta _0} + \theta } \right)}}{{\left( {{c_1} - 1} \right){a_1}}} + {a_1},}\\{{\lambda _1} = \frac{{2{S_1}}}{{{{\left( {{c_1} - 1} \right)}^2}a_1^2 - 2{S_1} \cdot \tan \left( {{\theta _0} + \theta } \right)}},}\end{array}} \right.$

(16)

$\left\{ {\begin{array}{*{20}{l}}{{S_2} = \left( {{a_2} - \frac{{\Delta x}}{2}} \right)d - \int_0^{{a_2}} f (x) \cdot {\rm{d}}x,}\\{{x_{{B_2}}} = \frac{{2{S_2} \cdot \tan \left( {{\theta _0} - \theta } \right)}}{{\left( {{c_2} - 1} \right){a_2}}} + {a_2},}\\{{\lambda _2} = \frac{{2{S_2}}}{{{{\left( {{c_2} - 1} \right)}^2}a_2^2 - 2{S_2} \cdot \tan \left( {{\theta _0} - \theta } \right)}}.}\end{array}} \right.$

(17)

由图 7可解得a₁、a₂与压深d的函数关系表示如下：

${a_1} = \left\{ \begin{array}{l}\sqrt {(2r - d)d} ,\\\;\;\;0 \le d < r - r\sin \left( {{\theta _0} + \theta } \right);\\(d - r)\tan \left( {{\theta _0} + \theta } \right) + \frac{r}{{\cos \left( {{\theta _0} + \theta } \right)}},\\\;\;\;d \ge r - r\sin \left( {{\theta _0} + \theta } \right).\end{array} \right.$

(18)

${a_2} = \left\{ \begin{array}{l}\sqrt {\left( {2r - d} \right)d} ,\\\;\;\;0 \le d < r - r\sin \left( {{\theta _0} - \theta } \right);\\\left( {d - r} \right)\tan \left( {{\theta _0} - \theta } \right) + \frac{r}{{\cos \left( {{\theta _0} - \theta } \right)}},\\\;\;\;d \ge r - r\sin \left( {{\theta _0} - \theta } \right).\end{array} \right.$

(19)

如图 7所示，当压深d=0时，轮齿的轮廓线方程可表示如下：

$f\left( x \right) = \left\{ \begin{array}{l}{f_1}\left( x \right) = \frac{{ - x}}{{\tan \left( {{\theta _0} + \theta } \right)}} + r\left( {1 - \frac{1}{{\sin \left( {{\theta _0} + \theta } \right)}}} \right),\;\;\;\;x < - r\cos \left( {{\theta _0} + \theta } \right);\\{f_2}\left( x \right) = {f_3}\left( x \right) = r - \sqrt {{r^2} - {x^2}} ,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;x \in \left[ { - r\cos \left( {{\theta _0} + \theta } \right),r\left( {{\theta _0} - \theta } \right)} \right];\\{f_4}\left( x \right) = \frac{x}{{\tan \left( {{\theta _0} - \theta } \right)}} + r\left( {1 - \frac{1}{{\sin \left( {{\theta _0} - \theta } \right)}}} \right),\;\;\;x > \cos \left( {{\theta _0} - \theta } \right).\end{array} \right.$

(20)

因而式(16)和(17)中塑性变形区的积分部分面积可分别按下式计算：

$\int_{ - {a_1}}^0 f (x) \cdot {\rm{d}}x = \left\{ \begin{array}{l}\int_{ - {a_1}}^0 {{f_2}} (x) \cdot {\rm{d}}x,\quad \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0 \le {a_1} \le r\cos \left( {{\theta _0} + \theta } \right);\\\int_{ - {a_1}}^{ - r\cos \left( {{\theta _0} + \theta } \right)} {{f_1}(x) \cdot {\rm{d}}x} + \int_{ - r\cos \left( {{\theta _0} + \theta } \right)}^0 {{f_2}} (x) \cdot {\rm{d}}x,\;\;\;\;{a_1} > r\cos \left( {{\theta _0} + \theta } \right),\end{array} \right.$

(21)

$\int_0^{{a_2}} f (x) \cdot dx = \left\{ {\begin{array}{*{20}{l}}{\int_0^{{a_2}} {{f_3}} (x) \cdot {\rm{d}}x,}&{0 \le {a_2} \le r\cos \left( {{\theta _0} - \theta } \right);}\\{\mathop \smallint \nolimits_0^{r\cos \left( {{\theta _0} - \theta } \right)} {f_3}(x) \cdot {\rm{d}}x + \int_{r\cos \left( {{\theta _0} - \theta } \right)}^{{a_2}} {{f_4}} (x) \cdot {\rm{d}}x,}&{{a_2} > r\cos \left( {{\theta _0} - \theta } \right).}\end{array}} \right.$

(22)

当轮齿载荷位于β、α区域时，右侧材料没有塑性变形，此时右侧塑性变形区面积S₂=0，x_B2=a₂，λ₂=0。相应地，当轮齿载荷位于ω、φ区域时，左侧材料没有塑性变形，此时左侧塑性变形区面积S₁=0，x_B1=a₁，λ₁=0。
3.2 宽度系数的求解根据楔形压头滑移线解法^[19]解得宽度系数c₁、c₂表示如下：

$\left\{ {\begin{array}{*{20}{l}}{{c_1}\left( {{c_1} - 2} \right)\sin \left( {{\theta _0} + \theta } \right) = \sqrt {c_1^2 - {{\cos }^2}\left( {{\theta _0} + \theta } \right)} ,}\\{{c_2}\left( {{c_2} - 2} \right)\sin \left( {{\theta _0} - \theta } \right) = \sqrt {c_2^2 - {{\cos }^2}\left( {{\theta _0} - \theta } \right)} .}\end{array}} \right.$

(23)

综上，通过式(3)—(23)的推导过程建立了轮齿斜压管壁时所受的法向力F_N和切向力F_T与轮齿压深d和滑移量Δx之间的函数关系。
4 区域临界线求解及模型的应用4.1 区域临界线求解式(4)和(5)所示的γ₁、γ₂区域力学模型反映的是轮齿左右两侧材料同时有塑性变形时载荷与压痕形貌之间的关系。当轮齿切向力F_T增大导致轮齿开始沿右侧面向左下方运动时, 此时右侧材料开始无塑性变形但承受压力时，轮齿载荷作用开始从γ₁区域进入β区域，其临界线F_T0L计算如下：

$\left\{ \begin{array}{l}令\;\Delta {x_{0{\rm{L}}}} = {d_{0{\rm{L}}}} \cdot \tan \left( {{\theta _0} - \theta } \right),\\由\;{F_{\rm{N}}} = {f_1}\left( {{d_{0{\rm{L}}}},\Delta {x_{0{\rm{L}}}}} \right)\;解出\;{d_{0{\rm{L}}}},\\则\;{F_{{\rm{T0L}}}} = {f_2}\left( {{d_{0{\rm{L}}}},\Delta {x_{0{\rm{L}}}}} \right) = \\{f_2}\left( {{d_{0{\rm{L}}}},{d_{0{\rm{L}}}} \cdot \tan \left( {{\theta _0} - \theta } \right)} \right).\end{array} \right.$

(24)

当切向力F_T进一步增大导致轮齿只有左侧材料有塑性变形，轮齿右侧面虽未与管壁脱离接触，但开始不受力时，轮齿载荷作用开始从β区域进入α区域，其临界线F_T1L计算如下：

$\left\{ \begin{array}{l}令\;\Delta {x_{1{\rm{L}}}} = {d_{1{\rm{L}}}} \cdot \tan \left( {{\theta _0} - \theta } \right),\\由\;{F_{\rm{N}}} = {f_{{\rm{1L}}}}\left( {{d_{1{\rm{L}}}},\Delta {x_{1{\rm{L}}}}} \right)\;解出\;{d_{1{\rm{L}}}},\\则\;{F_{{\rm{T1L}}}} = {f_{{\rm{2L}}}}\left( {{d_{1{\rm{L}}}},\Delta {x_{1{\rm{L}}}}} \right) = \\{f_{{\rm{2L}}}}\left( {{d_{1{\rm{L}}}},{d_{1{\rm{L}}}} \cdot \tan \left( {{\theta _0} - \theta } \right)} \right).\end{array} \right.$

(25)

在γ₁和γ₂区域，轮齿左右两侧管壁材料均有塑性变形，所不同的是γ₁区域轮齿的Δx为正值(图 4中向左)，γ₂区域轮齿的Δx为负值(图 4中向右)，故γ₁和γ₂区域的临界线F_T0Z表征的状态为压痕两侧材料均有塑性变形且轮齿垂直压入管壁，水平滑移量Δx=0，F_T0Z计算如下：

$\left\{ \begin{array}{l}令\;\Delta {x_{0C}} = 0;\\由\;{F_{\rm{N}}} = {f_1}\left( {{d_{0{\rm{C}}}},0} \right)\;解出\;{d_{0C}};\\则\;{F_{{\rm{T0Z}}}} = {f_2}\left( {{d_{0C}},0} \right).\end{array} \right.$

(26)

当轮齿切向力F_T减小导致轮齿开始沿左侧面向右下方运动时, 此时左侧材料开始无塑性变形但承受压力时，轮齿载荷作用开始从γ₂区域进入ω区域，其临界线F_T0R计算如下：

$\left\{ {\begin{array}{*{20}{l}}{令\;\Delta {x_{0{\rm{R}}}} = - {d_{0{\rm{R}}}} \cdot \tan \left( {{\theta _0} + \theta } \right),}\\{由\;{F_{\rm{N}}} = {f_1}\left( {{d_{{\rm{0R}}}},\Delta {x_{0{\rm{R}}}}} \right)\;解出\;{d_{0{\rm{R}}}},}\\{则\;{F_{{\rm{T0R}}}} = {f_2}\left( {{d_{{\rm{0R}}}},\Delta {x_{0{\rm{R}}}}} \right) = }\\{{f_2}\left( {{d_{{\rm{0R}}}}, - {d_{0{\rm{R}}}} \cdot \tan \left( {{\theta _0} + \theta } \right)} \right).}\end{array}} \right.$

(27)

当切向力F_T进一步减小导致轮齿只有右侧材料有塑性变形，轮齿左侧面虽未与管壁脱离接触，但开始不受力时，轮齿载荷作用开始从ω区域进入φ区域，其临界线F_T1R计算如下：

$\left\{ {\begin{array}{*{20}{l}}{令\;\Delta {x_{1{\rm{R}}}} = - {d_{1{\rm{R}}}} \cdot \tan \left( {{\theta _0} + \theta } \right),}\\{由\;{F_{\rm{N}}} = {f_{1{\rm{R}}}}\left( {{d_{1{\rm{R}}}},\Delta {x_{1{\rm{R}}}}} \right)\;解出\;{d_{1{\rm{R}}}},}\\{则\;{F_{{\rm{T}}1{\rm{R}}}} = {f_{2{\rm{R}}}}\left( {{d_{1{\rm{R}}}},\Delta {x_{1{\rm{R}}}}} \right) = }\\{{f_{2{\rm{R}}}}\left( {{d_{1{\rm{R}}}}, - {d_{1{\rm{R}}}} \cdot \tan \left( {{\theta _0} + \theta } \right)} \right).}\end{array}} \right.$

(28)

4.2 力学模型的应用利用上述斜压力学模型可以完成正向和逆向问题的求解。当已知轮齿的法向载荷F_N和切向载荷F_T需要预测和计算压深d和滑移量Δx时，首先需要根据式(24)—(28)计算出各作用区域的临界值，据此判断当前载荷处在哪个区域，然后按照该区域的计算模型求解压深d和滑移量Δx，具体流程如图 8所示。

图 8 已知载荷求解压痕形貌流程图

图选项

当已知压深d和滑移量Δx需要追溯形成该压痕的载荷时，首先计算Δx与d之间的比例系数C，合理C值需满足：-tan(θ₀+θ)≤C≤tan(θ₀-θ), 不合理的压痕是由α区域或φ区域或其他非正常过程产生的。当C处于正常区间内值时，对应γ₁或γ₂作用区域，此时轮齿两侧面材料均有塑性变形，载荷有唯一解。当C为正常区间的边界值时，对应β或ω作用区域。从式(13)可知，当已知轮齿载荷F_N和F_T时，由于在β区域轮齿运动满足：Δx=d·tan(θ₀-θ), 式左侧为单一参数d的函数，可据此求解出压痕形貌的唯一解。但当已知压痕的形貌求解轮齿载荷时，式(13)右侧为轮齿载荷F_N和F_T满足左侧约束的一个解的集合，此时不能依据式(13)求解出β区域轮齿载荷的唯一解。同理，不能依据式(15)求解出ω区域轮齿载荷的唯一解。已知压痕形貌求解载荷的具体流程如图 9所示。

图 9 已知压痕形貌求解载荷流程图

图选项

5 仿真验证5.1 验证方法为了验证所建立的斜压力学模型，用ABAQUSP有限元分析软件对该问题进行了仿真计算。仿真采用显式算法(Explicit)；针对接触区域，利用自适应网格划分技术(ALE)以解决大变形问题导致的网格畸变；驱动轮轮齿采用解析刚体, 齿顶角2θ₀=100°，斜压角度θ=5°，圆角半径r=0.05 mm；轮齿厚度b=10 mm，管壁材料为铝合金，材料屈服强度σ_s=205.78 MPa(材料达到屈服后屈服强度保持不变)，Young's模量E=71 000 MPa，Possion比μ=0.33。网格划分用的是平面应变、减缩积分四边形单元(CPE4R)。边界条件如图 10所示，管壁下表面固定，上表面与轮齿接触，轮齿施加法向力F_N和切向力F_T。根据式(23)求解的压痕宽度系数C₁=3.201，C₂=3.383。据Mises屈服准则，取k=118.81 MPa。为对比分析各不同作用区域内轮齿载荷与压痕形貌的不同变化规律，各作用区域的法向力F_N统一按下式加载：

图 10 有限元仿真模型

图选项

${F_{\rm{N}}} = 50000t,\;\;\;\;0 \le t \le 0.1.$

(29)

其中：t为仿真时间，单位为s；法向力F_N和切向力F_T的单位为N。
5.2 γ₁区域验证结果为了验证γ₁区域的力学计算模型，按下式施加切向力：

${F_{\rm{T}}} = 6000t,\;\;\;\;0 \le t \le 0.1.$

(30)

当轮齿载荷处于γ₁区域时，轮齿左右两侧面材料均有塑性变形，轮齿将向左下方移动，滑移量为正值，仿真验证结果如图 11所示。

图 11 γ1区域仿真验证结果

图选项

5.3 γ₂区域验证结果为了验证γ₂区域的力学计算模型，按下式施加切向力：

${F_{\rm{T}}} = - 100t,\;\;\;\;0 \le t \le 0.1.$

(31)

当轮齿载荷处于γ₂区域时，轮齿左右两侧面材料均有塑性变形，轮齿将向右下方移动，滑移量为负值，仿真验证结果如图 12所示。

图 12 γ2区域仿真验证结果

图选项

对比γ₁、γ₂区域的载荷、压深和滑移量可知，当轮齿左右两侧面材料均有塑性变形时，法向力主要对压深产生影响，而切向力主要对水平滑移量产生影响。同时由于γ₁、γ₂区域所能施加的切向力相对较小，则轮齿所能输出的牵引力也相对较小，故当轮齿载荷处于γ₁、γ₂区域时其牵引效率较低。
5.4 β区域验证结果为了验证β区域的力学计算模型，按下式施加切向力：

${F_{\rm{T}}} = \left\{ {\begin{array}{*{20}{l}}{2000t,}&{0 \le t \le 0.005;}\\{10 + 18000(t - 0.005),}&{0.005 < t \le 0.01;}\\{100 + 40000(t - 0.01),}&{0.01 < t \le 0.015;}\\{300 + 16000(t - 0.015),}&{0.015 < t \le 0.02;}\\{380 + 19375(t - 0.02),}&{0.02 < t \le 0.1.}\end{array}} \right.$

(32)

当轮齿载荷处于β区域时，轮齿将沿右侧面向左下方移动，滑移量为正值，仿真验证结果如图 13所示。

图 13 β区域仿真验证结果

图选项

对比β区域和γ₁区域的载荷、压深和滑移量可知，当轮齿只有左侧材料有塑性变形时，除法向力主要对压深产生影响外，切向力也会令压深发生变化。由于β区域能施加较大的切向力，轮齿能输出较大的牵引力，故为提高潜油电机的利用效率，应控制合适的法向力F_N，在避免电机出现堵转的前提下，使轮齿所受的切向力F_T位于β区域，并尽可能靠近β区域与α区域的临界线F_T1L。
5.5 ω区域验证结果为了验证ω区域的力学计算模型，按下式施加切向力：

${F_{\rm{T}}} = \left\{ {\begin{array}{*{20}{l}}{8000t,}&{0 \le t \le 0.005;}\\{40 + 18000(t - 0.005),}&{0.005 < t \le 0.01;}\\{130 + 9666.67(t - 0.01),}&{0.01 < t \le 0.1.}\end{array}} \right.$

(33)

当载荷处于ω区域时，轮齿将沿左侧面向右下方移动，滑移量为负值，仿真验证结果如图 14所示。

图 14 ω区域仿真验证结果

图选项

对比ω区域和γ₂区域的载荷、压深和滑移量可知，当轮齿只有右侧材料有塑性变形时，法向力主要对压深产生影响，而切向力除了对水平滑移量产生影响外也会令压深发生变化。
6 试验验证6.1 试验条件及方法试验装置如图 15所示，其中丝杆螺母驱动移动平台沿竖直方向运动，用以施加法向力；拉紧螺栓驱动铝板沿水平方向运动，用以施加切向力，通过调整驱动轮的固定角度设置轮齿压入铝板的角度。

图 15 试验装置实物图

图选项

试验前采用WDW3020型电子万能材料试验机对试验用铝板材料进行拉伸试验，测得试验铝板塑性阶段的平均应力值为327.9 MPa。驱动轮齿顶角为80°，圆角半径为0.05 mm，轮齿厚度为10 mm，斜压角度为10°。选取爬行器井下工作常用负载状态对应的β区域进行验证。共进行2组压痕试验，其F_N分别为6 129 N和6 251 N，F_T分别为2 015 N和1 653 N，试验后铝板压痕如图 16所示。

图 16 铝板压痕试验结果

图选项

6.2 试验结果及分析试验结束后，采用NexView三维白光干涉形貌仪对铝板压痕进行测量，根据测量得到的压痕表面形貌的三维点云，在MATLAB中选取具有相同Y坐标的数据点，即可绘制该截面下的压痕轮廓。根据压痕轮廓修正后的轮齿压入角度为9.8°，据式(23)解得压痕宽度系数C₁=3.284，C₂=3.940，采用图 8所示流程求解已知载荷条件下的压痕形貌，压痕形貌试验测量结果和理论计算结果的对比如图 17所示。从压痕形状上看，铝板材料左侧隆起较高，右侧几乎没有隆起发生，这与前述对β区域的理论分析是一致的。从压痕具体形貌参数上看，理论模型较为准确地计算了压入深度，2组压入深度实际值分别为0.357 mm和0.362 mm，而理论模型计算值分别为0.401 mm和0.395 mm，误差分别为12.1%和9.1%。

图 17 压痕试验验证结果

图选项

但是从图 17中可以看出，左侧隆起高度的计算结果误差较大，分析原因主要有以下几个方面：
1) 模型假设材料是理想刚塑性，而铝板实际变形是存在一定的弹性变形的，在卸载之后会发生回弹现象，塑性变形量的减小将导致隆起高度降低；此外，通过材料拉伸试验发现，在应变逐渐增大的过程中材料会发生硬化现象，强度不断增大会导致相同力的作用下压入深度减小，塑性变形区域减小也将导致隆起高度降低。
2) 由于试验装置中导轨等并非无摩擦，摩擦力会导致力传感器测量的数值偏大，该部分误差也会导致理论计算的隆起高度相比实际值偏大。
3) 无论是有限元仿真还是试验结果均表明，两侧隆起并非是直线相交，因为交点处应力值会发生奇异，而是会产生一定圆弧过渡，从而导致实际隆起高度降低。
7 结论与展望本文根据F_T和F_N之间的不同比例关系，将轮齿倾斜压入管壁材料时的载荷条件分为6个不同作用区域，利用滑移线理论分析和计算了各作用区域内轮齿接触部分的压力分布函数，并建立了压痕形貌与载荷之间的计算模型。通过有限元仿真和压痕试验发现：当轮齿左右两侧面材料均有塑性变形时，法向力主要对压深产生影响，而切向力主要对水平滑移量产生影响；当轮齿只有单侧材料有塑性变形时，除法向力主要对压深产生影响外，切向力也会令压深发生变化。
斜压过程的力学模型为多齿牵引力建模提供了齿间塑性变形位移约束方程，为研究某一个转动周期内轮齿的受力与压痕轮廓的函数关系奠定了基础。但是由于对材料弹性变形和硬化的忽略以及试验装置摩擦等的影响，导致试验实际隆起高度与计算值存在较大误差，后续研究中，可以将模型中的材料参数视为关于Young's模量、屈服强度和硬化指数的变量，并利用有限元仿真结果进行拟合，对模型进行改进。

参考文献

[1]	HEDDLESTON D C. Horizontal well production logging deployment and measurement techniques for US land shale hydrocarbon plays[R]. Oklahoma City: Society of Petroleum Engineers, 2009.
[2]	HALLUNDBEK J. Well tractors for highly deviated and horizontal wells[R]. Oklahoma City: Society of Petroleum Engineers, 1994.
[3]	姜岳庆.水平井测井及仪器输送方法研究[D].大庆: 东北石油大学, 2013. JIANG Y Q. Study on horizontal well logging tools and transport method of tools[D]. Daqing: Northeast Petroleum University, 2013. (in Chinese)
[4]	刘清友, 李维国. Sondex水平井井下爬行器的研究与应用[J]. 石油钻采工艺, 2008, 30(5): 115-117. LIU Q Y, LI W G. Research and application of Sondex downhole tractor for Horizontal wells[J]. Oil Drilling & Production Technology, 2008, 30(5): 115-117. DOI:10.3969/j.issn.1000-7393.2008.05.027 (in Chinese)
[5]	李湘涛, 秦羽乔, 陈四平, 等. 水平井测井仪器输送技术及其应用[J]. 石油机械, 2014, 42(8): 98-102. LI X T, QIN Y Q, CHEN S P, et al. The horizontal well logging instrument conveying technology and its applications[J]. China Petroleum Machinery, 2014, 42(8): 98-102. DOI:10.3969/j.issn.1001-4578.2014.08.023 (in Chinese)
[6]	周劲辉, 张勇, 李翠. 水平井自扶正式电缆牵引器的设计[J]. 石油机械, 2015, 43(2): 79-82. ZHOU J H, ZHANG Y, LI C. Design of self-righting cable tractor in horizontal well[J]. China Petroleum Machinery, 2015, 43(2): 79-82. (in Chinese)
[7]	SHEIRETOV T. Wireline tractors and mechanical services tools: comparative study of technical solutions[C]//SPE/ICOTA Coiled Tubing and Well Intervention Conference and Exhibition, Houston, USA: SPE, 2016.
[8]	OMARI M A, PLESSING H. Innovation in coiled tubing tractor technology extends the accessibility of coiled tubing in horizontal wells, allowing better possibilities for well intervention[C]//SPE Middle East Oil and Gas Show and Conference. Manama, Bahrain: SPE, 2007.
[9]	HENDERSON B, HOPWOOD C, HAMILTON C, et al. Cost saving benefits of using a fully bi-directional tractor system[C]//SPE/CIM International Conference on Horizontal Well Technology. Calgary, Alberta, Canada: SPE, 2000.
[10]	AL-SHAWLY A S, AL-BUALI M H, AL-OMRAN M R, et al. Wireline well tractor technology experience in extended reach horizontal well[C]//Abu Dhabi International Petroleum Exhibition and Conference. Abu Dhabi, UAE: SPE, 2010.
[11]	刘振.水平井牵引器机构分析和建模仿真技术研究[D].哈尔滨: 哈尔滨工业大学, 2008. LIU Z. Research on mechanism analysis and modeling simulation of horizontal well tractor[D]. Harbin: Harbin Institute of Technology, 2008. (in Chinese)
[12]	曾华军.水平井牵引器驱动系统关键技术研究[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)
[13]	于淼.水平井牵引器拖动力数值分析[D].大庆: 东北石油大学, 2011. YU M. The mule downhole mumerical analysis of drag force on[D]. Daqing: Northeast Petroleum University, 2011. (in Chinese)
[14]	侯松伟.水平井牵引器运动学仿真与应用研究[D].大庆: 东北石油大学, 2017. HOU S W. Research on kinematics simulation and application for horizantal well tractor[D]. Daqing: Northeast Petroleum University, 2017. (in Chinese)
[15]	张勇.水平井电缆牵引器的机械设计研究[D].北京: 中国石油大学(北京), 2016. ZHANG Y. Research on mechanical design of cable tractor in horizontal wells[D]. Beijing: China University of Petroleum (Beijing), 2016. (in Chinese)
[16]	刘清友, 李雨佳, 任涛, 等. 水平井爬行器驱动轮力学分析[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)
[17]	秦浩.水平井爬行器驱动机构结构优化[D].成都: 西南石油大学, 2017. QIN H, Optimization of horizontal crawler's driving mechanism[D]. Chengdu: Southwest Petroleum University, 2017. (in Chinese)
[18]	HILL R. The mathematical theory of plasticity[M]. Oxford: Larendon Press, 1950.
[19]	JOHNSON K L. Contact mechanics[M]. Beijing: Higher Education Press, 1992.
[20]	GRUNZWEIG J, LONGMAN I M, PETCH N J. Calculations and measurements on wedge-indentation[J]. Journal of the Mechanics and Physics of Solids, 1954, 2(2): 81-86. DOI:10.1016/0022-5096(54)90002-2
[21]	RICE J R. Plane strain slip line theory for anisotropic rigid/plastic materials[J]. Journal of the Mechanics and Physics of Solids, 1973, 21(2): 63-74. DOI:10.1016/0022-5096(73)90030-6
[22]	BAI Y L, 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
[23]	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): 47. DOI:10.1007/s11249-015-0524-3
[24]	JACKSON R L. A solution of rigid perfectly plastic cylindrical indentation in plane strain and comparison to elastic-plastic finite element predictions with hardening[J]. Journal of Applied Mechanics, 2018, 85(2): 24501. DOI:10.1115/1.4038495
[25]	常旭, 杨东超, 孙可平, 等. 爬行器驱动轮正压过程分析[J]. 清华大学学报(自然科学版), 2019, 59(7): 537-543. CHANG X, YANG D C, SUN K P, et al. Analysis of the normal indentation process caused by tractor wheels[J]. Journal of Tsinghua University (Science and Technology), 2019, 59(7): 537-543. (in Chinese)