技术领域本发明涉及一种用于由弧齿锥齿轮轴上功率谱向齿面应力谱精确转换的方法，尤其适用于复杂工况下的齿轮传动系统的齿面应力谱分布规律的计算分析。
背景技术：
在弧齿锥齿轮传动系统的寿命设计与预测分析工作中，需要以真实工况弧齿锥齿轮的齿面应力谱分布为依据，目前较为流行的弧齿锥齿轮传动系统齿面应力谱获取方法通常为实验测量法。实验测量法为在试验台上的弧齿锥齿轮的轮齿根部粘贴应变片，通过对运行中的轮齿根部应变信号进行采集，并经处理得到弧齿锥齿轮的齿根弯曲应力谱分布，但实验测量法所需设备昂贵，人力、物力投入大，试验周期长。针对弧齿锥齿轮的齿面接触应力谱的测量方面，还没有实际可行的实验测量方法。与此同时，在理论计算方面，目前尚缺乏一种基于实际工况下的精确、简单方便的实际工况下的弧齿齿轮传动系统的齿根弯曲和齿面接触应力谱分布计算分析方法，在弧齿锥齿轮传动系统的设计阶段就可以计算得到弧齿锥齿轮的齿面应力谱分布规律，进而可以为传动系统的寿命预测和设计工作提供重要的依据。
技术实现要素：
本发明要解决的技术问题为：针对给定复杂工况下的弧齿锥齿轮载荷谱分布分析问题，综合考虑弧齿锥齿轮在真实工况下的边界条件，通过建立弧齿锥齿轮的动力学模型的建立，实现了弧齿锥齿轮动载系数的求解，利用弧齿锥齿轮LTCA加载接触ABAQUS有限元模型，计算出真实工况下的弧齿锥齿轮齿面接触应力和齿根弯曲应力，结合工况时间比，利用混合正态分布参数估计法对齿轮随机应力谱分布规律进行估计，为齿轮传动系统的寿命设计和预测提供了一种简便的由弧齿锥齿轮传动系统轴上功率谱向齿面应力谱精确转换的方法，有效降低了载荷谱测试实验所带来的难度和成本。本发明采用的技术方案是：一种由弧齿锥齿轮传动系统轴上功率谱向齿面应力谱精确转换的方法，该方法步骤如下：步骤(1)、根据轴上功率谱计算得到各个工况下的弧齿锥齿轮工作扭矩；步骤(2)、计算各个工况下的弧齿锥齿轮动载系数；步骤(3)、建立弧齿锥齿轮加载接触ABAQUS有限元分析模型，添加载荷，偏移量和动载系数边界条件，进行有限元模型计算；步骤(4)、从各个工况的有限元计算结果中提取齿面接触应力值和齿根弯曲应力值；步骤(5)、根据功率谱中的各个工况时间比，将齿面接触应力值和齿根弯曲应力值分别扩展得到齿面整个工作过程的随机应力谱；步骤(6)、对齿面随机应力谱分布利用混合正态分布法进行估计，得到齿面应力谱分布规律。进一步的，所述步骤(1)中根据功率、转速和扭矩的换算公式，将功率谱各个工况下的功率、转速计算得到相应的扭矩T：T＝9550P/n式中P为功率(kw)；n为转速(r/min)；T为扭矩(N*m)。进一步的，所述步骤(2)中的弧齿锥齿轮动载荷系数式中Fn为齿轮间的动态啮合力，T为齿轮工作额定扭矩，r为齿轮啮合点半径。进一步的，所述步骤(2)中弧齿锥齿轮动载系数的计算中，齿轮间的动态啮合力Fn可以通过建立弧齿锥齿轮的动力学模型并求解计算得到，动力学模型中主动轮1、被动轮2都具有沿着x、y、z方向的移动自由度和绕各自轴向的转动自由度θ1x、θ2y，齿轮动力学方程为：m1X··1+cx1X·1+kx1X1=Fxm1Y··1+cy1Y·1+ky1Y1=Fym1Z··1+cz1Z·1+kz1Z1=FzJ1θ1x··=T1-Fzr1m2X··2+cx2X·2+kx2X2=-Fxm2Y··2+cy2Y·1+ky1Y1=-Fym2Z··2+cz2Z·2+kz2Z2=-FzJ2θ2y··=-T2+Fzr2]]>式中m1、m2为主、被动轮的质量，J1、J2为两齿轮的转动惯量，cij、kij(i＝x,y,z；j＝1,2)为主、被动轮沿x、y、z方向的阻尼系数和刚度系数，T1、T2为输入力矩和负载力矩，Fx、Fy、Fz为齿轮间的动态啮合力Fn在x、y、z方向的分力：Fn=kh(t)λn+chλ·nFx=-Fn(sinαncosδ2+cosαnsinβ2sinδ2)Fy=Fn(sinαnsinδ2-cosαnsinβ2cosδ2)Fz=Fncosαncosβ2]]>式中kh(t)为弧齿锥齿轮的时变啮合刚度，λn为齿轮啮合点间由振动和误差产生的位移，ch为啮合阻尼，αn为法面压力角，β2为被动轮的中点螺旋角，δ2为被动轮的节锥角。进一步的，所述步骤(2)中的弧齿锥齿轮动载系数的计算中，齿轮啮合点间由振动和误差产生的位移为：λn＝(x1-x2)cosδ1sinαn-(y1-y2)cosδ1cosαnsinβ1-[(z1+r1θ2y)-(z2+r2θ1x)]cosαncosβ1-en(t)式中en(t)为齿轮副的法向传动静态误差，xi、yi、zi(i＝1为主动轮，i＝2为被动轮)为齿轮沿x、y、z方向的移动位移，θ1x为主动轮绕轴x转动位移，θ2y为被动轮绕轴y的转动位移，δ1为主动轮的节锥角，αn为法面压力角，β1为主动轮的中点螺旋角。进一步的，所述步骤(2)中弧齿锥齿轮动载系数的计算中，弧齿锥齿轮动力学模型的方程组为非线性微分方程组，编制五阶变步长自适应龙格库塔求解程序，设定动态响应求解时间t，进行计算得到λn和进而可以得到Fn，通过kv和Fn的关系可以得到齿轮的动载系数。进一步的，所述步骤(3)中的弧齿锥齿轮加载接触ABAQUS有限元分析模型为五齿分析模型，网格单元为六面体单元C3D8R，载荷边界条件为各个工况下的扭矩，幅值曲线设置为动载系数的周期变化曲线，分析步设为稳态分析步。进一步的，所述步骤(4)中从有限元分析结果中提取动载系数分别处于最大值、最小值、平均值的齿面点的一个啮合周期内的接触应力和齿根弯曲应力变化值，分别进行叠加得到该工况下的齿面接触应力值和齿根弯曲应力值。进一步的，所述步骤(5)中各工况时间比η是各工况时间占总时间的百分比，将步骤(4)中提取的各工况下的齿面接触应力值和齿根弯曲应力值扩展η×10000次，将扩展后的齿面接触应力值和齿根弯曲应力值分别进行叠加，得到弧齿锥齿轮整个工作过程的齿面接触随机应力谱和齿根弯曲随机应力谱。进一步的，所述步骤(6)中对齿面随机应力谱的分布拟合采用混合正态分布参数估计法，需要设定的系统参数有：l为需要估计的正态分布的个数，r_w为单峰值极小值位置求解半径，espi为精度容差系数。进一步的，所述步骤(6)中利用混合正态分布参数估计法计算得出齿面随机应力谱概率分布公式为式中：w_li为各正态分布的权重系数，μi为各正态分布的期望，σi为各正态分布的标准差，l为需要估计的正态分布的个数。本发明的原理：利用各工况的下的转速n、矩T结合齿轮的参数通过动力学模型求解各工况下的动载系数kv，通过建立弧齿锥齿轮LTCA加载接触ABAQUS有限元模型，控制边界条件，精确地计算出各工况下的齿面接触应力和齿根弯曲应力，根据各工况的时间百分比，将齿面接触应力和齿根弯曲应力进行扩展叠加得到相应的随机应力谱，再通过混合正态分布参数估计法对齿面随机应力谱分布进行拟合。本发明与现有技术相比的有益效果是：首先，应用本发明不需要经过专门的齿轮载荷谱试验直接通过该理论分析计算就可以精确地获取齿面应力谱分布规律，为齿轮传动系统的寿命设计与预测工作提供依据，有效地减少了试验所带来的人力和物力成本；其次，目前尚未有较成熟的由功率谱向齿轮应力谱精确转换方法，本发明可操作性强、精确度高，只需经过相应的分析计算，即可得到弧齿锥齿轮在真实工况下的齿面应力谱分布规律。附图说明图1为弧齿锥齿轮动力学模型。图2为工况1的弧齿锥齿轮动载系数响应图。图3为弧齿锥齿轮LTCA有限元加载接触ABAQUS有限分析模型。图4为工况1的齿面应力数据图，其中a为齿面接触应力数据，b为齿根弯曲应力数据。图5为齿面应力谱的混合正态分布参数估计图，其中a为齿面接触应力谱的概率分布估计图，b为齿根弯曲应力谱的概率分布估计图。图6为本发明的方法流程图。具体实施方式下面结合附图以及具体实施例进一步说明本发明。本发明一种由弧齿锥齿轮传动系统轴上功率谱向齿面应力谱精确转换的方法，其方法流程如下：步骤(1)、根据轴上功率谱计算得到各个工况下的弧齿锥齿轮工作扭矩；步骤(2)、计算各个工况下的弧齿锥齿轮动载系数；步骤(3)、建立弧齿锥齿轮加载接触ABAQUS有限元分析模型，添加载荷，偏移量和动载系数边界条件，进行有限元模型计算；步骤(4)、从各个工况的有限元计算结果中提取齿面接触应力值和齿根弯曲应力值；步骤(5)、根据功率谱中的各个工况时间比，将接触应力和齿根弯曲应力变化值分别扩展得到齿面整个工作过程的随机应力谱；步骤(6)、对齿面随机应力谱分布利用混合正态分布法进行估计，得到齿面应力谱分布规律。所述步骤(1)中根据功率、转速和扭矩的换算公式，将功率谱各个工况下的功率、转速计算得到相应的扭矩T：T＝9550P/n式中P为功率(kw)；n为转速(r/min)；T为扭矩(N*m)。所述步骤(2)中的弧齿锥齿轮动载荷系数式中Fn为齿轮间的动态啮合力，T为齿轮工作额定扭矩，r为齿轮啮合点半径。所述步骤(2)中弧齿锥齿轮动载系数的计算中，齿轮间的动态啮合力Fn可以通过建立弧齿锥齿轮的动力学模型并求解计算得到，动力学模型中主动轮1、被动轮2都具有沿着x、y、z方向的移动自由度和绕各自轴向的转动自由度θ1x、θ2y，齿轮动力学模型的动力学方程为：m1X··1+cx1X·1+kx1X1=Fxm1Y··1+cy1Y·1+ky1Y1=Fym1Z··1+cz1Z·1+kz1Z1=FzJ1θ1x··=T1-Fzr1m2X··2+cx2X·2+kx2X2=-Fxm2Y··2+cy2Y·1+ky1Y1=-Fym2Z··2+cz2Z·2+kz2Z2=-FzJ2θ2y··=-T2+Fzr2]]>式中m1、m2为主、被动轮的质量，J1、J2为两齿轮的转动惯量，cij、kij(i＝x,y,z；j＝1,2)为主、被动轮沿x、y、z方向的阻尼系数和刚度系数，T1、T2为输入力矩和负载力矩，Fx、Fy、Fz为齿轮间的动态啮合力Fn在x、y、z方向的分力：Fn=kh(t)λn+chλ·nFx=-Fn(sinαncosδ2+cosαnsinβ2sinδ2)Fy=Fn(sinαnsinδ2-cosαnsinβ2cosδ2)Fz=Fncosαncosβ2]]>式中kh(t)为弧齿锥齿轮的时变啮合刚度，λn为齿轮啮合点间由振动和误差产生的位移，ch为啮合阻尼，αn为法面压力角，β2为被动轮的中点螺旋角，δ2为被动轮的节锥角。所述步骤(2)中的弧齿锥齿轮动载系数的计算中，齿轮啮合点间由振动和误差产生的位移为：λn＝(x1-x2)cosδ1sinαn-(y1-y2)cosδ1cosαnsinβ1-[(z1+r1θ2y)-(z2+r2θ1x)]cosαncosβ1-en(t)式中en(t)为齿轮副的法向传动静态误差，xi、yi、zi(i＝1为主动轮，i＝2为被动轮)为齿轮沿x、y、z方向的移动位移，θ1x为主动轮绕轴x转动位移，θ2y为被动轮绕轴y的转动位移，δ1为主动轮的节锥角，αn为法面压力角，β1为主动轮的中点螺旋角。所述步骤(2)中弧齿锥齿轮动载系数的计算中，弧齿锥齿轮动力学模型的方程组为非线性微分方程组，编制五阶变步长自适应龙格库塔求解程序，设定动态响应求解时间t，进行计算得到λn和进而可以得到Fn，通过kv和Fn的关系，可以得到齿轮的动载系数。所述步骤(3)中的弧齿锥齿轮加载接触ABAQUS有限元分析模型为五齿分析模型，网格单元为六面体单元C3D8R，载荷边界条件为各个工况下的扭矩，幅值曲线设置为动载系数的周期变化曲线，分析步设为稳态分析步。所述步骤(4)中从有限元分析结果中提取动载系数分别处于最大值、最小值、平均值的齿面点的一个啮合周期内的接触应力和齿根弯曲应力变化值，分别进行叠加得到齿面接触应力和齿根弯曲应力3个啮合周期应力变化值曲线。所述步骤(5)中各工况时间比η是各工况时间占总时间的百分比，将步骤(4)中提取的各工况下的接触应力和齿根弯曲应力值扩展η×10000次，将扩展后的齿面接触应力值和齿根弯曲应力值分别进行叠加，得到弧齿锥齿轮整个工作过程中的齿面接触随机应力谱和齿根弯曲随机应力谱。所述步骤(6)中对齿面随机应力谱的分布拟合采用混合正态分布参数估计法，需要设定的系统参数有：l为需要估计的正态分布的个数，r_w为单峰值极小值位置求解半径，espi为精度容差系数。所述步骤(6)中利用混合正态分布参数估计法计算得出齿面随机应力谱概率分布公式为式中：w_li为各正态分布的权重系数，μi为各正态分布的期望，σi为各正态分布的标准差，l为需要估计的正态分布的个数。具体的，本发明的流程图如图6所示。下面以表1所示的某弧齿锥齿轮传动系统为例，具体说明本发明方法，但本发明的保护范围不限于下述实例：表1某弧齿锥齿轮传动系统齿轮参数步骤(1)：根据轴上功率谱计算得到各个工况下的弧齿锥齿轮工作扭矩；某弧齿锥齿轮传动系统的功率谱具有12个工况，根据功率、转速和扭矩的换算公式，利用各个工况下的功率、转速计算得到相应的扭矩：T＝9550P/n式中P为功率(kw)；n为转速(r/min)；T为扭矩(N*m)。该弧齿锥齿轮传动系统的功率谱和扭矩如表2所示。表2某弧齿锥齿轮传动系统功率谱和扭矩步骤(2)：计算各个工况下的弧齿锥齿轮动载系数；以工况1为例，弧齿锥齿轮动载荷系数式中Fn为齿轮间的动态啮合力，T为齿轮工作额定扭矩，其中工况1的输入力矩T1＝622.6KN，负载力矩T2＝1842KN；r为齿轮啮合点半径，主动轮半径r1＝26.2mm，被动轮半径r2＝185.3mm。齿轮间的动态啮合力Fn可以通过建立弧齿锥齿轮的动力学模型并求解计算得到，动力学模型如图1所示，主动轮1、被动轮2都具有沿着x、y、z方向的移动自由度和绕各自轴向的转动自由度θ1x、θ2y，齿轮动力学模型的动力学方程为：m1X··1+cx1X·1+kx1X1=Fxm1Y··1+cy1Y·1+ky1Y1=Fym1Z··1+cz1Z·1+kz1Z1=FzJ1θ1x··=T1-Fzr1m2X··2+cx2X·2+kx2X2=-Fxm2Y··2+cy2Y·1+ky1Y1=-Fym2Z··2+cz2Z·2+kz2Z2=-FzJ2θ2y··=-T2+Fzr2]]>式中m1、m2为主、被动轮的质量，其中m1＝4.65kg，m2＝9.65kg；J1、J2为两齿轮的转动惯量，其中J1＝0.009kg·m2，J2＝0.166kg·m2；cij、kij(i＝x,y,z；j＝1,2)为主、被动轮沿x、y、z方向的阻尼系数和刚度系数，其中cij＝5e5N/mm/s，kij＝1e7N·mm；Fx、Fy、Fz为齿轮间的动态啮合力Fn在x、y、z方向的分力：Fn=kh(t)λn+chλ·nFx=-Fn(sinαncosδ2+cosαnsinβ2sinδ2)Fy=Fn(sinαnsinδ2-cosαnsinβ2cosδ2)Fz=Fncosαncosβ2]]>式中kh(t)为弧齿锥齿轮的时变啮合刚度，kh(t)＝1e7(1+0.2cos(4569.6·t))；ch为啮合阻尼，αn为法面压力角，20°；β为中点螺旋角，35°；δ2为被动轮的节锥角，77°；λn为齿轮啮合点间由振动和误差产生的位移：λn＝(x1-x2)cosδ1sinαn-(y1-y2)cosδ1cosαnsinβ1-[(z1+r1θ2y)-(z2+r2θ1x)]cosαncosβ1-en(t)式中en(t)为齿轮副的法向传动静态误差，en(t)＝3.6e-2·cos(4569.6·t)，单位mm。弧齿锥齿轮动力学模型的方程组为非线性微分方程组，编制五阶变步长自适应龙格库塔求解程序，设定动态响应求解时间t，进行计算得到λn和进而可以得到Fn，通过kv和Fn的关系，可以得到工况1的弧齿锥齿轮的动载系数kv如图2所示。以同样的方法计算其他工况下的齿轮动载系数。步骤(3)：建立弧齿锥齿轮加载接触ABAQUS有限元分析模型，添加载荷，偏移量和动载系数边界条件，进行有限元模型计算；弧齿锥齿轮加载接触ABAQUS有限元分析模型为五齿分析模型，网格单元为六面体单元C3D8R，载荷边界条件为各个工况下的扭矩，幅值曲线设置为动载系数的周期变化曲线，分析步设为稳态分析步。ABAQUS有限元分析模型如图3所示。步骤(4)：从各个工况的有限元计算结果中提取齿面接触应力值和齿根弯曲应力值；从有12个工况的有限元分析结果中提取动载系数分别处于最大值、最小值、平均值的齿面点的一个啮合周期内的接触应力和齿根弯曲应力变化值，分别进行叠加得到齿面接触应力和齿根弯曲应力3个啮合周期过程中的应力变化值曲线。工况1的3个啮合周期过程中的接触应力值和齿根弯曲应力值变化曲线如图4所示。步骤(5)：根据功率谱中的各个工况时间比，将接触应力和齿根弯曲应力变化值分别扩展得到齿面整个工作过程的随机应力谱。各工况时间比η是各工况时间占总时间的百分比，将12个工况的接触应力和齿根弯曲应力值扩展η×10000次，最后将扩展后的12个工况下的齿面接触应力值和齿根弯曲应力值分别进行叠加，得到弧齿锥齿轮整个工作过程中的齿面接触随机应力谱和齿根弯曲随机应力谱。步骤(6)：对齿面随机应力谱进行分析，得到齿面应力谱分布规律。对齿面随机应力谱的分布拟合采用混合正态分布参数估计法进行估计，设定系统参数：l＝4为需要估计的正态分布的个数，r_w＝5为单峰值极小值位置求解半径，espi＝0.01为精度容差系数。计算得出齿面随机应力谱概率分布估计公式为式中：w_li为各正态分布的权重系数，μi为各正态分布的期望，σi为各正态分布的标准差。齿面接触应力谱的分布估计为：w_l1＝0.335，μ1＝5.2，σ1＝2.8；w_l2＝0.079，μ2＝67.0，σ2＝2.1；w_l3＝0.532，μ3＝77.0，σ3＝53.8；w_l4＝0.054，μ4＝154.4，σ4＝1.6。齿面接触应力谱的分布估计如图5a所示。齿面弯曲应力谱的分布估计为：w_l1＝0.247，μ1＝29.6，σ1＝11.5；w_l2＝0.183，μ2＝125.3，σ2＝13.2；w_l3＝0.165，μ3＝291.1，σ3＝13.5；w_l4＝0.405，μ4＝155.2，σ4＝55.3。齿根弯曲应力谱的分布估计如图5b所示。总之，本发明针弧齿锥齿轮传动系统由轴上功率谱向齿面应力谱精确转换的问题，利用轴上功率谱计算得到各工况下的扭矩，建立弧齿锥齿轮动力学模型得到各工况下的齿轮动载系数，通过弧齿锥齿轮ABAQUS有限模型，设定扭矩和动载系数边界条件，计算得到各个工况下的齿面点在啮合周期内的齿根弯曲应力和齿面接触应力变化值，根据各工况的时间比，对应力变化值进行扩展最后叠加得到齿面整个工作过程中的随机应力谱，通过混合正态分布参数估计法对齿面应力谱分布规律进行估计，从而为弧齿锥齿轮的寿命预测和设计工作提供重要的依据。