本发明属于齿轮传动技术领域，具体涉及一种螺旋锥齿轮的含误差齿面接触分析方法。
背景技术
螺旋锥齿轮因其重叠系数大、承载能力强、传动比高、传动平稳、噪声小等优点广泛应用于汽车、航空、矿山等机械传动领域。螺旋锥齿轮分为两种，一种是弧齿锥齿轮，其大轮轴线和小轮轴线相交；另一种是准双曲面螺旋锥齿轮，其大轮轴线和小轮轴线有一定偏置距。目前生产螺旋锥齿轮世界级公司主要是美国的格里森公司和瑞士的奥利康公司。
齿面接触分析方法(toothcontactanalysis，tca)作为模拟齿轮副齿面接触过程并判定齿面接触性能的关键技术之一，其数值求解的稳定性尤为重要，它与齿面接触点尤其是初始点接触点的确定尤为相关。在传统的齿面接触分析求解过程中，初值问题往往定义为估值问题，即都是人为估算设定，就有很大的偶然性和不确定性；另外，目前求解建立的非线性齿面接触分析方程的数值算法，忽略了制造装配误差，而在实际齿轮传动过程中，制造装配误差是不得不考虑的关键因素，它能直接影响齿面的接触性能如齿面接触区域、齿面传动误差等。
技术实现要素：
本发明的目的在于提供一种螺旋锥齿轮的含误差齿面接触分析方法，提出了精确的初值解决方案，将制造装配误差作为约束条件，以完成螺旋锥齿轮齿面接触点的精确求解，通过etca方程组完成对齿面接触性能的评价。
本发明提供的这种螺旋锥齿轮的含误差齿面接触分析方法，包括以下步骤：
(1)求解初始接触状态时的大轮齿面初值，设定含装配误差的约束条件，求解小轮齿面初值，将大轮齿面初值、小轮齿面初值作为齿面接触分析求解的精确初值，求解整个齿面的初始接触点；
(2)获取坐标转化矩阵，得到etca方程组，建立所述etca方程组的评价项：a)齿面接触印痕；b)传动误差；c)齿面误差，完成对螺旋锥齿轮齿面接触性能的评价。
在一个具体实施方式中，所述步骤(1)中，具体步骤为：
1.1)求解大轮齿面初值(φ2,θ2)；
以大轮齿面中点m作为参考点，齿面中点m到大轮轴线的距离为r2，齿面中点m到轴线的投影点与齿面交叉顶点o2的距离为l2，根据齿面的参数化离散方法，建立的基本坐标系{i,j,k}，可表示为：
i＝[100],j＝[010],k＝[001](1)
而pi是关于齿面中点的锥角δm的函数，有：
从齿面中点m到大端面的距离为f[t]b，b为齿面宽，ha2和hf2可根据以下公式求解：
ha2＝hae2-f[t]btanθa2(3)
hf2＝hfe2-f[t]btanθa2(4)
由于齿面点设定为齿面中点m，则齿面中点m到节锥的距离为：
齿面中点m到根锥的距离为：
齿面中点m到大轮轴线的距离为：
r2＝(re-fb)sinδ2-△h1cosδ2＝(re-f[t]b)sinδ2-1/2(hf2-ha2)cosδ2(7)
系数f[t]是判定条件，当取0.25，0.5，0.75时，分别对应的是齿轮大端、中端和小端的齿面参考点；
齿面中点m到轴线的投影点与齿面交叉顶点o2的距离为：
以大轮齿面中点m坐标为已知条件，来求解大轮齿面初值(φ2,θ2)，则存在以下关系：
显然，r2和l2都可以表示成关于φ2和θ2的参数化函数，齿面参数化建模中r2是具有非线性的，通过非线性迭代算法求解对应的(φ2，θ2)；
1.2)求解小轮齿面初值(φ1,θ1)；
在小轮齿面初值求解过程中，考虑制造装配因素作为齿面接触时约束条件，在安装误差的基本定义中，e为大小齿轮的轴线偏差，p为小轮沿轴向的偏移距离，g为大轮沿轴线的偏移距离，α为两轴线的轴交角；
当大小齿轮合理装配在一起时，设大轮绕轴线p2旋转，小轮绕轴线p1旋转，则它们的装配误差可以表示为：
而两个轴交点op和og的距离为：
在齿轮初始接触时刻，传动误差几乎为零，大轮和小轮的传动比是恒定的且等于z2/z1，由于传动误差定义为输出轮的实际转动位置和理论值之差
△φg，输出大轮的旋转角度差值φg可表示为输入小轮角度的φp近似抛物线函数，由此可得，齿面的装配直接影响齿面误差，则在初始接触时刻，小轮齿面点满足以下条件：
同理，利用非线性迭代算法求解该方程，得到小轮齿面初值(φ1,θ1)；
然后将大轮齿面初值、小轮齿面初值作为齿面接触分析方程求解的精确初值x(0)＝(φ1(0),φ2(0),θ1(0),θ2(0))。
在一个具体实施方式中，所述步骤(1)中，利用凸模型自适应信赖域算法，求解大轮齿面初值、小轮齿面初值和齿面初始接触点p0*(φ1,θ1,φ2,θ2)，求解这三个参数的过程相似。
在一个具体实施方式中，利用凸模型自适应信赖域算法求解小轮齿面初值，具体步骤为：
为了求解方程(14)，可将其转化成一个无约束非线性优化问题；
利用凸模型，信赖域的子问题可以表示为：
式中，fk＝f(xk)和gk＝g(xk)表示梯度，bk∈rn表示水平向量，且
为了选择凸模型的参数bk和bk，采用：
如果ρk≥0，有：
否则，βk＝1。在置换信赖域过程中，对于给定的足够大的正数b*>0，有
很显然，bk需要满足于：
另外，已经证明：当ρk>0,则dtk-1yk-1>0；当bk-1>0时，bk>0；在目前的搜索中，给定一个信赖域：
式中，rτ(rk)(τ∈(0,1))是一个关于信赖域半径增加或减小的松弛r-函数，有：
式中β，γ1，γ2和ω是相应的常数，关于r-函数的详细定义可以参考文献，另外，rk表示函数f和fk的增益比为：
在一个具体实施方式中，凸模型自适应信赖域算法的基本程序如下：
输入x0∈rn,△0>0,b0∈rn×1,0<β<1,0<γ1<1-β,γ2>0,m>1+γ2,ε>0,0<μ<1；
初始化:k＝0,△0＝1,x＝x0,b0＝0,b0＝i(单元矩阵)；
step1.评价fk＝f(xk)和
step2.ifk≥1,计算ρk通过求解方程(17)；ifρk≥0,求解βk根据方程(18)；
otherwise,βk＝1.利用方程(19)-(21)分别求解bk,yk-1,和bk.
step3.设定b＝bk,g＝gk,及b＝bk.获得近似解d*,通过求解方程(17).
step4.假定dk＝d*,根据方程(25)计算rk.
step5.ifrk<μ,令xk+1＝xk；otherwise,令xk+1＝xk+dk.
step6.修正参数bk和bk,求得bk+1和bk+1；通过方程(22)与(23)计算△k+1.
step7.令k＝k+1,返回初始化.
在一个具体实施方式中，所述步骤(2)中，具体步骤为：
2.1)获取坐标转化矩阵，得到etca方程组；
通过齿面建模进行齿面离散化，大小齿面点pi(φi,θi)(i＝1,小轮；i＝2,大轮)可在各自的坐标系oi(xi,yι,zι)(i＝1,2)求解出来，其点矢和法矢分别表示为ri(i＝1,2)和ni(i＝1,2)，然后两个齿面需要各自旋转一定角度达到齿面共轭点接触状态，完成齿面的接触传动；在此初始接触点p*(φ1,θ1,φ1,θ2)位置，满足齿面啮合原理，有：
其中，该旋转过程有旋转矩阵mi-f表示，etca接触条件主要包括：i)旋转后两个齿面点拥有共同的点矢，即(rf)1＝(rf)2；ii)旋转后两个齿面点的法矢共线，即n1＝χ0n2，χ0为常数，由此可得旋转矩阵mi-f(i＝1,2)；
首先，考虑大小齿轮是在各自坐标系下建立的齿面模型，则它们必须旋转一定的角度到达统一齿面啮合坐标系。则从齿轮原坐标系oi(xi,yi,zi)到啮合坐标系of(xf,yf,zf)的转换矩阵mf-i∈r4(i＝1,2)可表示为：
它的子矩阵mx-2为：
旋转角度为：
而且：
到达同一个啮合坐标系of(xf,yf,zf)后，由于制造装配误差的存在，大小齿面可能没有发生点接触，则需要进一步旋转变换，约束齿面接触运动过程，且满足tca接触条件，其间采用以下旋转方式：即旋转小轮齿面，而大轮齿面加入装配误差。在此旋转小轮齿面过程中，旋转矩阵可表示为：
(mg-x)1＝(mg-a)1·(ma-b)1·(mb-c)1·mc-x(32)
在实现法矢共线的旋转变换过程中，包含mg-a,ma-b,mb-c和mc-x等运动链，有：
式中，旋转角度
在实现点矢重合的转换过程中，其旋转矩阵(mt-f)1与矩阵mg-x的过程一致，包含的子运动链为mt-g,mg-a,ma-b,mb-c和mc-f，齐次变换过程为：
除了旋转角度
经过整个坐标变换后，小轮的齿面点可表示为：
nf:n1＝mx-1(mg-x)1·n1(38)
rf:(rf)1＝mf-1(mt-f)1·r1(39)
值得说明的是，在制造装配误差的考虑中，一般把误差值加在其中一个齿面中即可，其考虑装配误差的坐标变换矩阵为：
其用于计算法矢共线过程的子矩阵mm-a∈r3
把制造装配误差加入大轮齿面，通过所有旋转变换后，可得：
rf:(rf)2＝mm-amf-2·r2(42)
nf:n2＝mm-amx-2·ni(43)
综上所述，满足发现共线的etca条件，可得以下两个标量方程：
(fetca)1:(n1)1+(n2)1＝0(44)
(fetca)2:(n1)2+(n2)2＝0(45)
此处，反映的是ni＝[(ni)1,(ni)2,(ni)3]的其中两个分量的函数关系；
同理，满足点矢重合的etca条件时，有以下标量方程成立：
(fetca)3:[(rf)1]1-[(rf)2]1＝0(46)
(fetca)4:[(rf)1]2-[(rf)2]2＝0(47)
(fetca)5:[(rf)1]3-[(rf)2]3＝0(48)
此时，关于(rf)i:＝{[(rf)i]1,(rf)i]2,(rf)i]3}t(i＝1,2)的三个分量要同时考虑；
假定在齿轮传动中输入小轮的角速度为ω1，大轮的角速度为ω2，通过旋转变换后大轮的的角速度变为ωf2，有以下关系成立：
大轮与小轮的相对角速度为：
ω1-2＝ω1-ωf2(50)
考虑制造装配误差度的影响，则经过坐标转换后，有：
至此，大轮的线速度v2和小轮的线速度v1为：
则大轮与小轮的相对线速度v1-2为：
v1-2＝v1-v2(54)
在此初始接触状态，满足啮合原理，可得：
(fetca)6:(n2)1×(v1-2)1+(n2)2×(v1-2)2+(n2)3×(v1-2)3＝0(55)
综合以上，六个标量方程含有6个未知变量
2.2)建立etca方程组的评价项：a)齿面接触印痕；b)传动误差；c)齿面误差，完成对螺旋锥齿轮齿面接触性能的评价。
本发明的有益技术效果：
(1)本文所述螺旋锥齿轮齿面接触分析方法，提出了精确的etca初值解决方案，并在求解过程中考虑制造装配误差作为基本约束条件，匹配精确的etca方程组建立过程，给定精确的非线性求解算法，完成了整个齿面接触点尤其是初始接触点的鲁棒性精确求解。
(2)本发明所述螺旋锥齿轮齿面接触分析方法，采用凸模型自适应信赖域算法，具有稳定性、强收敛性和自适应性的特点，能有效解决雅克比矩阵的病态问题并只需较少迭代就可获得理想数值结果；同时，区别于传统二阶近似曲面模型，该凸模型在处理强非二阶行为、显著曲率变化等问题方面有更好的优势。
附图说明
图1为利用凸模型自适应信赖域法(cstar)求解大小轮齿面初值的三次计算的迭代质量评价。
图2为为三次迭代计算的信赖域半径置换情况。
图3为求解初始接触点p0*的基本流程。
图4为cstar、nr和ph法得到的初始点的结果比较。
图5为etca方程组的评价结果：(a.1)齿面接触点；(a.2)齿面印痕；(b)传动误差；(c)齿面误差。
图6为制造装配误差对etca评价结果的影响。
具体实施方式
下面将对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围，下面结合附图和具体实施例对本发明进一步说明。
本发明实施例提供一种螺旋锥齿轮的含误差齿面接触分析方法，包括以下步骤：
步骤一、求解初始接触状态时的大轮齿面初值，设定含装配误差的约束条件，求解小轮齿面初值，将大轮齿面初值、小轮齿面初值作为齿面接触分析求解的精确初值，求解整个齿面的初始接触点p0*(φ1,θ1,φ2,θ2)；
表1为端面铣削弧齿锥齿轮的齿面设计基本参数，表2为sgm调整卡的加工参数，即大轮采用展成法，小轮采用变性法加工，在齿轮传动过程中，大轮凹面和小轮凸面发生接触，这两个齿面的参数将被用来进行etca的计算；
表1弧齿锥齿轮齿面设计基本参数
表2弧齿锥齿轮sgm调整卡参数
1.1)求解大轮齿面初值(φ2,θ2)
以大轮齿面中点m作为参考点，齿面中点m到大轮轴线的距离为r2，齿面中点m到轴线的投影点与齿面交叉顶点o2的距离为l2，根据齿面的参数化离散方法，建立的基本坐标系{i,j,k}，可表示为：
i＝[100],j＝[010],k＝[001](1)
而pi是关于齿面中点的锥角δm的函数，有：
从齿面中点m到大端面的距离为f[t]b，b为齿面宽，ha2和hf2可根据以下公式求解：
ha2＝hae2-f[t]btanθa2(3)
hf2＝hfe2-f[t]btanθa2(4)
由于齿面点设定为齿面中点m，则齿面中点m到节锥的距离为：
齿面中点m到根锥的距离为：
齿面中点m到大轮轴线的距离为：
r2＝(re-fb)sinδ2-△h1cosδ2＝(re-f[t]b)sinδ2-1/2(hf2-ha2)cosδ2(7)
系数f[t]是判定条件，当取0.25，0.5，0.75时，分别对应的是齿轮大端、中端和小端的齿面参考点；
齿面中点m到轴线的投影点与齿面交叉顶点o2的距离为：
以大轮齿面中点m坐标为已知条件，来求解大轮齿面初值(φ2,θ2)，则存在以下关系：
显然，r2和l2都可以表示成关于φ2和θ2的参数化函数，齿面参数化建模中r2是具有非线性的，通过相应的迭代方法求解对应的(φ2,θ2)。
1.2)求解小轮齿面初值(φ1,θ1)
在小轮齿面初值求解过程中，考虑制造装配因素作为齿面接触时约束条件，在安装误差的基本定义中，e为大小齿轮的轴线偏差，p为小轮沿轴向的偏移距离，g为大轮沿轴线的偏移距离，α为两轴线的轴交角；
当大小齿轮合理装配在一起时，设大轮绕轴线p2旋转，小轮绕轴线p1旋转，则它们的装配误差可以表示为：
而两个轴交点op和og的距离为：
在齿轮初始接触时刻，传动误差几乎为零，大轮和小轮的传动比是恒定的且等于z2/z1，由于传动误差定义为输出轮的实际转动位置和理论值之差△φg，输出大轮的旋转角度差值φg可表示为输入小轮角度的φp近似抛物线函数。由此可得，齿面的装配直接影响齿面误差，则在初始接触时刻，小轮齿面点满足以下条件：
同理，利用非线性迭代算法求解该方程，得到小轮齿面初值(φ1,θ1)；
表3为利用凸模型自适应信赖域法(cstar)求解得到的大轮齿面初值、小轮齿面初值，它们被用作齿面接触分析的精确初值x(0)＝(φ1(0),φ2(0),θ1(0),θ2(0))；
表3齿面接触分析的精确初值
图1为利用凸模型自适应信赖域法(cstar)求解大小轮齿面初值的三次计算的迭代质量评价，在求解大轮齿面初值、小轮齿面初值的过程中，cstar法计算量小，收敛快，关于目标函数f(x)、迭代步长和一阶优化都能显示出较好的迭代质量。
在求解初始接触点p0*的迭代求解过程中，由于计算的复杂性，需要更大的计算量，但依然能保证较快的收敛速度。图2为三次迭代计算的信赖域半径置换情况，表示了在这次计算中的信赖域半径，在其初始化中所用到的基本系数为γ1＝γ2＝0.15，μ＝0.25，β＝0.1和ε＝10-14。
图3为求解初始接触点p0*的基本流程，为了进一步验证算法的正确性，将本发明与newton-raphson(nr)和powell混合(ph)算法进行对比，并与文献中的bfgs(broyden-fletcher-goldfarb-shanno)算法进行比较，bfgs算法求解的初始接触点为(0.993766,0.1796032,0.86776021,-0.09447587)作为参考值，图4为cstar、nr和ph法得到的初始点的结果比较，通过数值比较可得：本发明cstar算法更具优势，求得的初始点p0*(0.98991007,0.0177254,0.96769247,-0.0940923)及其φ2，θ2，φ1和θ1的计算精度分别是99.612％,98.692％,99.993％及99.594％。
步骤二、获取坐标转化矩阵，得到etca方程组，建立etca方程组的评价项：i)齿面接触印痕；ii)传动误差；iii)齿面误差，完成对螺旋锥齿轮齿面接触性能的评价；
通过齿面建模进行齿面离散化，大小齿面点pi(φi,θi)(i＝1,小轮；i＝2,大轮)可在各自的坐标系oi(xi,yι,zι)(i＝1,2)求解出来，其点矢和法矢分别表示为ri(i＝1,2)和ni(i＝1,2)，然后两个齿面需要各自旋转一定角度达到齿面共轭点接触状态，完成齿面的接触传动；在此初始接触点p*(φ1,θ1,φ1,θ2)位置，满足齿面啮合原理，有：
其中，该旋转过程有旋转矩阵mi-f表示，etca接触条件主要包括：i)旋转后两个齿面点拥有共同的点矢，即(rf)1＝(rf)2；ii)旋转后两个齿面点的法矢共线，即n1＝χ0n2，χ0为常数，由此可得旋转矩阵mi-f(i＝1,2)；
首先，考虑大小齿轮是在各自坐标系下建立的齿面模型，则它们必须旋转一定的角度到达统一齿面啮合坐标系。则从齿轮原坐标系oi(xi,yi,zi)到啮合坐标系of(xf,yf,zf)的转换矩阵mf-i∈r4(i＝1,2)可表示为：
它的子矩阵mx-2为：
旋转角度为：
而且：
到达同一个啮合坐标系of(xf,yf,zf)后，由于制造装配误差的存在，大小齿面可能没有发生点接触，则需要进一步旋转变换，约束齿面接触运动过程，且满足tca接触条件，其间采用以下旋转方式：即旋转小轮齿面，而大轮齿面加入装配误差。在此旋转小轮齿面过程中，旋转矩阵可表示为：
(mg-x)1＝(mg-a)1·(ma-b)1·(mb-c)1·mc-x(32)
在实现法矢共线的旋转变换过程中，包含mg-a,ma-b,mb-c和mc-x等运动链，有：
式中，旋转角度
在实现点矢重合的转换过程中，其旋转矩阵(mt-f)1与矩阵mg-x的过程一致，包含的子运动链为mt-g,mg-a,ma-b,mb-c和mc-f，齐次变换过程为：
除了旋转角度
经过整个坐标变换后，小轮的齿面点可表示为：
nf:n1＝mx-1(mg-x)1·n1(38)
rf:(rf)1＝mf-1(mt-f)1·r1(39)
值得说明的是，在制造装配误差的考虑中，一般把误差值加在其中一个齿面中即可，其考虑装配误差的坐标变换矩阵为：
其用于计算法矢共线过程的子矩阵mm-a∈r3
把制造装配误差加入大轮齿面，通过所有旋转变换后，可得：
rf:(rf)2＝mm-amf-2·r2(42)
nf:n2＝mm-amx-2·ni(43)
综上所述，满足发现共线的etca条件，可得以下两个标量方程：
(fetca)1:(n1)1+(n2)1＝0(44)
(fetca)2:(n1)2+(n2)2＝0(45)
此处，反映的是ni＝[(ni)1,(ni)2,(ni)3]的其中两个分量的函数关系；
同理，满足点矢重合的etca条件时，有以下标量方程成立：
(fetca)3:[(rf)1]1-[(rf)2]1＝0(46)
(fetca)4:[(rf)1]2-[(rf)2]2＝0(47)
(fetca)5:[(rf)1]3-[(rf)2]3＝0(48)
此时，关于(rf)i:＝{[(rf)i]1,(rf)i]2,(rf)i]3}t(i＝1,2)的三个分量要同时考虑；
假定在齿轮传动中输入小轮的角速度为ω1，大轮的角速度为ω2，通过旋转变换后大轮的的角速度变为ωf2，有以下关系成立：
大轮与小轮的相对角速度为：
ω1-2＝ω1-ωf2(50)
考虑制造装配误差度的影响，则经过坐标转换后，有：
至此，大轮的线速度v2和小轮的线速度v1为：
则大轮与小轮的相对线速度v1-2为：
v1-2＝v1-v2(54)
在此初始接触状态，满足啮合原理，可得：
(fetca)6:(n2)1×(v1-2)1+(n2)2×(v1-2)2+(n2)3×(v1-2)3＝0(55)
综合以上，六个标量方程含有6个未知变量
齿面初始接触点与各个齿面接触性能评价项的求解密切相关，确定精确齿面初始点的目的是为了进行etca结果评价，如图5所示，在小轮凸面，确定的齿面接触点和初始接触点p0*基本位于齿面中间区域，齿面接触印痕良好，主要分布在齿面中间部分，且没有出现过度的边缘接触，能够反映出一定的齿面承载能力；传动误差曲线呈对称分布，其最大幅值小于13arcsec，这能实现较低噪声和振动的齿轮传动；齿面误差分布合理，较小的误差集中在齿面中间区域，较大的主要分布在齿轮小端，平均绝对值为15.523μm，最大值为23.476μm，最小值为-24.138μm。
在本发明提出的齿面接触点的解析计算中，制造装配误差被论证为齿面初始接触点及齿面接触性能的一个主要影响因素，本发明重点给出制造装配误差对etca结果的影响，图6为制造装配误差对etca评价结果的影响，选取了两个评价项即齿面误差和传动误差，其中，g＝0作为约束条件且α一般为0，设定了5种关于装配误差项p和e的不同工况，得到的数据结果表明：制造装配误差对齿面误差和传动误差影响较大，尤其在工况2和工况3的取值区间内，变化最明显。因此，该方法可为实际制造过程中，考虑高性能齿轮传动和齿轮产品制造的制造装配误差的调整与优化提供思路。
以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例。对于本技术领域的技术人员来说，在不脱离本发明技术构思前提下所得到的改进和变换也应视为本发明的保护范围。