基于FFT和润滑影响的RV减速器摆线针轮残余应力求解方法与流程

本发明属于摆线针轮优化设计领域，涉及基于fft和润滑影响的rv减速器摆线针轮残余应力求解方法。

背景技术：

作为rv减速器关键零部件，摆线针轮以其结构紧凑、传动比大、传动精度高和承载能力大等优点被广泛应用在机器人传动、测试设备及卫星接收装置中。摆线针轮行星传动中摆线轮与针齿通过润滑剂进行润滑，以减小啮合过程中的摩擦磨损并传递动力，其润滑性能的好坏直接影响到摆线针轮的传递精度、效率、疲劳寿命及动态特性。在摆线轮与针齿啮合过程中，它们的材料在外载荷作用下存在一个应变极限，在低速重载等恶劣工况下，摆线轮与针齿表面或许会产生残余应力。当残余应力较大时，摆线轮或针齿表面会产生凹坑或裂纹，进而引起摆线针轮振动加剧，最终引起摆线轮与针齿表面产生齿面磨损、点蚀和胶合等常见形式的失效。此外，低速重载状况下，由于和残余应力相对应的残余变形与最小膜厚度处于同一量级，因此摆线针轮残余变形(残余应力)会对其润滑性能产生较大影响，进一步摆线针轮润滑性能也会通过润滑膜压力对其残余应力产生影响。

传统的rv减速器摆线针轮残余应力大多采用成熟的商用有限元软件来计算，且忽略了摆线轮与针轮加工精度和流体润滑的影响，这样简化处理虽然为摆线针轮的弹塑性研究带来了很大的方便，但是它无法准确、真实地反映出实际工况下的残余应力，从而无法准确预测摆线针轮的疲劳寿命。此外，考虑润滑影响的摆线针轮残余应力的计算涉及了润滑膜压力、载荷、塑性应变增量和残余变形等4种迭代计算，计算量很大，使用传统的计算方法严重耗费了计算时间和资源，甚至无法求解。所以，传统的rv减速器摆线针轮残余应力计算方法更多的是具有统计与经验意义，其精度上不够准确、效率上不够高效。

技术实现要素：

有鉴于此，本发明的目的在于提供一种基于fft和润滑影响的rv减速器摆线针轮残余应力求解方法。

为达到上述目的，本发明提供如下技术方案：

基于fft和润滑影响的rv减速器摆线针轮残余应力求解方法，该方法包括以下步骤：

s1：将摆线针轮前半啮合过程均分为n个啮合子过程，即将摆线针轮啮合相位角ψ＝π分为n等份δψ＝π/n，并基于齿廓形成方法和齿轮啮合原理计算摆线针轮第n啮合子过程，即ψn＝n·δψ，的啮合力fn、卷吸速度uen和当量曲率半径rxn，n＝1,…,n，以模拟摆线针轮的加载和卸载过程；

s2：建立考虑表面残余变形v^r影响的摆线针轮线接触弹流润滑ehl模型，并基于jacobi线性迭代法计算润滑膜压力分布p；

s3：采用影响系数法icm和fft算法计算摆线针轮次表面弹性应力并与次表面残余应力共同构成摆线针轮次表面总应力σij，定义ij为xx、yy、zz、xy、xz和yz；

s4：基于von-mises屈服准则准确判断摆线针轮塑性发生区域；

s5：根据摆线针轮总应力σij，利用newton-raphson法和塑性流动法则求解摆线针轮塑性发生区域的塑性应变增量并更新塑性应变

s6：基于摆线针轮塑性应变利用icm及fft更新次表面残余应力并结合次表面弹性应力更新摆线针轮次表面总应力σij；

s7：判断摆线针轮塑性应变增量是否收敛，若不收敛，更新塑性应变增量并返回步骤s5；若收敛，则转入步骤s8；

s8：根据摆线针轮塑性应变利用icm和fft快速计算摆线针轮表面残余变形v^r；

s9：判断残余变形v^r是否收敛，若不收敛，更新摆线针轮表面形貌并返回步骤s2重新计算润滑膜压力分布；若收敛，则转入步骤s10；

s10：判断摆线针轮是否啮合完成(啮合步数n是否等于n)，即加载和卸载是否完成，若未完成，摆线针轮继续啮合，令n＝n+1并返回步骤s1；若啮合完成，则结束循环。

可选的，所述步骤s2中，摆线针轮ehl模型中考虑了摆线针轮塑性性能的影响，即在膜厚公式中增加残余变形项h＝hc+hg+δ+v^e+v^rl，其中，δa和δa分别表示摆线轮和针轮表面粗糙度，v^e和v^r分别为摆线针轮综合弹性变形和表面残余变形分布。

可选的，所述步骤s4中，基于von-mises屈服准则判断摆线针轮塑性区域时，材料强化准则采用混合强化法则，包括各向同性和运动强化法则。

本发明的有益效果在于：本方法可在考虑流体润滑和加工精度影响的前提下快速地计算rv减速器摆线针轮残余应力，可利用不同润滑剂流变模型和reynolds方程差分方案准确地获得不同流体润滑状态下rv减速器摆线针轮的润滑膜压力，可用于计算弹流润滑和混合润滑等状态下rv减速器摆线针轮的残余应力，本发明提出的方法可为工程实践中rv减速器摆线针轮的优化设计和使用提供理论指导，对延长摆线针轮使用寿命、提高机械系统可靠性是有益的。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作优选的详细描述，其中：

图1为rv减速器摆线针轮二维结构示意图；

图2为摆线针轮应力计算域网格划分示意图；

图3为摆线针轮残余应力快速求解流程示意图。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

其中，附图仅用于示例性说明，表示的仅是示意图，而非实物图，不能理解为对本发明的限制；为了更好地说明本发明的实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。

本发明实施例的附图中相同或相似的标号对应相同或相似的部件；在本发明的描述中，需要理解的是，若有术语“上”、“下”、“左”、“右”、“前”、“后”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此附图中描述位置关系的用语仅用于示例性说明，不能理解为对本发明的限制，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。

典型的rv减速器摆线针轮的二维图如图1所示。其中，摆线轮齿数为zc，针轮齿数为zp，针轮中心圆半径为rp，针齿半径为rp，摆线轮偏心距为a，摆线轮节圆半径为rc，ni为摆线针轮输入转速。通常，摆线轮齿廓的形成方法主要有两种：两圆外啮合形成法和两圆内啮合形成法，基于上述两种齿廓形成方法可得摆线轮实际齿廓方程，即

其中，ψ为摆线针轮啮合相位角，k1为摆线轮短幅系数，k1＝zpa/rp。根据微分中曲率半径求解公式，得到第n啮合子过程摆线针轮等量曲率半径rxn为

式中，θb为滚圆相对基圆的相对转角，且θb＝zcψn。

对摆线针轮进行受力分析，得到第n啮合子过程针齿作用在摆线轮上的啮合力fn为

fn＝fmaxsinθn(3)

式中，θn为第n啮合子过程中啮合点法线与线的夹角，fmax为摆线针轮啮合过程最大啮合力，且mv为摆线针轮输出力矩，na为rv减速器中摆线轮个数，ζ为考虑摆线轮加工和安装误差引入的力矩修正系数，本发明取ζ＝1.1。

对摆线针轮进行受力分析，得到第n啮合子过程的润滑剂卷吸速度uen为

rv减速器摆线针轮啮合过程中，对单一摆线轮齿而言，它承受着由小增大又由大变小的循环变化载荷，且摆线针轮当量曲率半径和润滑剂卷吸速度也处于周期性变化的过程中。由于塑性问题是非线性的，并且与加载路径有关，因而将摆线针轮啮合过程均分成n个啮合子过程(加载/卸载步)以模拟摆线针轮啮合过程，在每一啮合子过程中摆线针轮的啮合力、当量曲率半径和卷吸速度保持不变，摆线针轮第n啮合子过程的啮合力用fn表示，当量曲率半径用rxn表示，卷吸速度用uen表示，n＝1,…,n。

以摆线针轮第n啮合子过程啮合力fn、当量曲率半径rxn和卷吸速度uen作为输入量，一种基于fft且考虑润滑作用的rv减速器摆线针轮残余应力快速计算方法通过以下步骤实现：

s1：基于摆线针轮第n啮合子过程啮合力fn、当量曲率半径rxn和卷吸速度uen，并利用jacobi迭代法计算考虑摆线针轮齿面残余变形v^r影响的摆线针轮润滑计算域(如图1所示)中润滑膜压力分布p；

s2：基于摆线针轮润滑计算域中的液膜压力分布p，利用icm及fft计算摆线针轮应力计算域中次表面弹性应力分布即

式中，q^x为摆线针轮表面沿x向的剪切力分布，和分别表示摆线针轮齿面润滑膜压力和表面剪切力对次表面弹性应力的影响系数矩阵，可由下式计算。

式中，和分别表示作用于半空间表面原点上的单位法向力或单位x向剪应力引起的弹性应力分布；ij表示应力分量xx、yy、zz、xy、xz和yz；(2αδx,2βδy,2γδz)和(2ξδx,2ψδy)分别为次表面应力和表面润滑膜压力作用点，δx和δy分别表示应力计算域沿x和y向网格半宽，如图2所示。

进一步，叠加摆线针轮应力计算域中的弹性应力和残余应力计算总的应力分布σij，即

s3：基于摆线针轮应力计算域中总应力分布σij，根据不同的材料强化法则并利用von-mises屈服准则判断摆线针轮塑性发生区域；

s4：根据von-mises屈服函数定义式f＝σvm-g(λ)，得到下面的摆线针轮有效累计塑性应变增量δλ的公式，并计算第n迭代过程中的δλ⁽ⁿ⁾。

其中，μ为摆线针轮材料的切向模量，kp和γ为运动强化参数。

获得有效累计塑性应变增量δλ后，计算rv减速器摆线针轮塑性应变增量，即

进而，更新摆线针轮塑性应变分布

s5：基于摆线针轮塑性应变分布利用fft和icm计算摆线针轮应力计算域中残余应力分布即

式中，b和p为摆线针轮中相应的塑性应变对残余应力的影响系数。

s6：判断是否满足收敛条件若满足，则结束摆线针轮塑性应变增量迭代；若不满足，则更新下列相关参数并返回步骤s4。其中，rv减速器摆线针轮塑性应变增量收敛精度为1.0×10^-4。

s7：利用icm及fft计算塑性应变ε^p引起的摆线针轮残余变形v^r，即

s8：判断摆线针轮表面残余变形v^r是否满足条件若不满足，则更新摆线针轮表面形貌并返回步骤s1；若满足，则转入步骤s9。其中，摆线针轮表面残余变形收敛精度取5.0×10^-4。

s9：判断摆线针轮啮合(加载/卸载)是否结束，若未结束，则进入下一啮合点(n＝n+1)并返回步骤s1；若加载已结束，则输出计算结果并终止循环。

上述过程可用图3所示的rv减速器摆线针轮残余应力快速求解流程图表示。

在步骤s1中，润滑膜压力分布的求解过程如下：

s101：基于第n啮合子过程中的摆线针轮啮合力fn、卷吸速度uen和当量曲率半径rxn，利用hertz接触原理计算摆线针轮的线接触参数，利用dowson-higginson(d-h)膜厚公式初始化考虑摆线针轮残余变形v^r的刚体中心膜厚hc，并利用green函数法计算润滑计算域中摆线针轮综合弹性变形系数d^e；

s102：根据摆线针轮不同啮合子过程对润滑计算域压力分布p进行初始化：第1啮合子过程的压力分布采用hertz接触理论进行初始化，后续啮合子过程的压力分布采用前一啮合子过程的计算结果进行初始化；

s103：利用icm和fft加速计算润滑计算域中摆线针轮综合弹性变形v^e，即

v^e＝d^e:p(12a)

式中，d^e摆线针轮综合弹性变形影响系数矩阵。

进一步，计算考虑摆线针轮表面粗糙度δ、残余变形v^rl和弹性变形v^e的润滑膜厚h，即

h(x,y)＝hc+hg(x,y)+δ(x,y)+v^e(x,y)+v^rl(x,y)(12b)

式中，δ表示摆线轮和针齿综合表面粗糙度；v^rl为润滑计算域中的残余变形分布，由应力计算域中的残余变形v^r利用二元三点插值法插值获得；hg为摆线针轮几何间隙。

s104：利用roelands粘压公式计算rv减速器摆线针轮润滑计算域中的粘度分布η，利用dowson-higginson密压公式计算摆线针轮润滑计算域中的密度分布ρ；

s105：利用gauss消去法求解reynolds方程计算rv减速器摆线针轮润滑膜压力分布p；

需要注意的是，本发明中所涉及的润滑模型适用于牛顿流体(润滑油)和非牛顿流体(润滑脂和ree-eyring流体)。假设摆线针轮中润滑剂密度和压力不沿膜厚方向变化并忽略润滑剂彻体力和惯性力，上述不同类型润滑剂的reynolds方程可统一表达为相同的形式，即

其中，εx、εy和ε是与润滑剂流变参数有关的积分系数。需要注意的是，上述reynolds方程在摆线针轮不同润滑区域中具有不同的形式：在弹流润滑区域，εx＞0且εy＞0；在接触区域，εx＝0且εy＝0。摆线针轮润滑计算域的压力边界条件为

其中，xi和xo分别表示rv减速器摆线针轮润滑计算域沿x向的入口与出口坐标，yin和yout分别表示摆线针轮润滑计算域沿y向的入口和出口坐标。

s106：判断摆线针轮润滑膜压力迭代是否收敛，相应的收敛准则为

式中，和分别为上一次和本次迭代计算得到的点(s,t)处的脂膜压力，errp为压力收敛精度，本发明取errp＝5×10^-5。若收敛，则结束摆线针轮润滑膜压力迭代；否则对摆线针轮润滑膜压力p进行修正，并返回步骤s103，压力修正公式为

式中，fitp为压力松弛因子，且fitp＝0.1；itp为压力迭代次数。

s107：利用数值积分法计算摆线针轮润滑膜承载力f，并判断润滑膜承载力与摆线针轮啮合力fn是否平衡，相应的收敛准则为

式中，errf为摆线针轮载荷收敛精度，且errf＝1×10^-3。若满足上式，则结束摆线针轮载荷迭代循环，并向应力计算域输出润滑膜压力分布；否则根据下式对摆线针轮刚体中心膜厚hc进行修正，并返回步骤s103。

式中，fitf为摆线针轮载荷松弛因子，且fitf＝5.0×10^-4；itp为摆线针轮载荷迭代次数。

在步骤s3中，判断摆线针轮塑形发生区域的具体过程如下：

s301：根据不同的摆线针轮材料强化法则，计算摆线针轮材料von-mises等效应力σvm。

(1)对于摆线针轮材料各向同性强化法则和无强化模型(理想弹塑性模型)：

(2)对于摆线针轮材料运动强化法则：

其中，为摆线针轮材料的偏应力；xij为摆线针轮背应力分布，与塑性变形的历史有关。

s302：根据不同摆线针轮材料强化法则，计算摆线针轮材料局部屈服强度函数g(λ)，即

(1)对于摆线针轮材料无强化模型(理想弹塑性模型)或运动强化法则：

g(λ)＝σy(18a)

(2)对于摆线针轮材料各向同性强化法则：

式中，σy为摆线针轮材料初始屈服强度，b、c和e为摆线针轮材料指数强化参数，et为摆线针轮材料弹塑性切向模量，e为摆线针轮材料弹性模量，λ为摆线针轮材料的有效累积塑性应变，

s303：对于更加一般的摆线针轮混合强化法则(包括各向同性强化法则和运动强化法则)，计算摆线针轮材料的von-mises屈服函数f＝σvm-g(λ)，并判断屈服函数f是否满足f＞0，若满足则摆线针轮在此处发生塑性变形，否则摆线针轮在此处没有发生塑性变形。

在摆线针轮卸载/加载过程中，有效累积塑性应变增量dλ和von-mises屈服函数f满足kuhn-tucker条件，即

f≤0,dλ≥0,f·dλ＝0(19)

此外，摆线针轮材料塑性应变的变化服从塑性流动准则，即

需要注意的是，利用各向同性强化法则，用于表征材料屈服强度的π平面大小改变，而位置和形状保持不变；对于运动强化法则，π平面大小和形状保持不变，位置改变。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。