CN104408220A - 一种改进的轮齿加载接触分析方法 - Google Patents

Abstract

本发明公开了一种改进的轮齿加载接触分析方法。该方法基于变形协调条件，通过几何计算得到齿轮的齿面方程和法向量方程；并在此基础上进行无载荷的轮齿接触分析，得到齿轮啮合的接触印痕、传动误差，齿面的主曲率、主方向，相切齿面的诱导法曲率，从而建立起接触椭圆大小、齿面压缩量与载荷三者之间的关系；结合Hertz理论公式和柔度系数法，考虑安装误差、支撑误差的影响和接触柔度系数的非线性，求出齿轮的载荷分布；根据齿间载荷分配结果，采用Hertz公式计算小大齿轮某时刻下的齿面压力分布和接触应力分布，利用三维有限元模型进行受力计算得到可靠的弯曲应力。轮齿加载接触分析方法对齿轮学的设计、计算和应用有重要意义。

Description

一种改进的轮齿加载接触分析方法

技术领域

本发明属于齿轮结构设计分析领域，具体地说，涉及一种改进的轮齿加载接触分析方法。

背景技术

轮齿加载接触分析(LTCA)方法在国内外已引起了广泛的重视，但大量的研究还停留在轮齿接触分析的基础上，如何在齿轮研究领域中建立好几何设计与力学分析的桥梁，特别是获得齿轮在拟真实工况下的工作性能，对齿轮学的设计、计算和发展具有重要意义。

现有公开的技术文献“齿轮轮齿承载接触分析(LTCA)的模型和方法”(机械传动，1998，22(2)，p1-3)中,基于柔度系数法和变形协调条件提出了一种轮齿加载接触分析(LTCA)方法，用于计算轮齿接触应力和弯曲应力，但该方法忽略了接触柔度系数的非线性影响，以及将接触区域载荷近似处理为线接触载荷来替代实际的区域载荷，当实际接触区域不是细长形区域，该方法的计算精度将变得更差，故不具有普遍性和较好的计算精度；其次，该方法将线接触载荷离散为一些点上的集中载荷，当离散点取得较多时，计算变量较多，计算复杂，对算法的收敛性要求高。

发明内容

为了避免传统TCA方法忽略轮齿承载接触变形和弯曲变形的影响，能够真实有效地反应轮齿在加载情况下的齿间载荷分布情况，齿面载荷分布情况，传动误差以及精确的计算轮齿接触应力和弯曲应力，本发明提出一种改进的轮齿加载接触分析方法。该方法基于变形协调条件，结合Hertz理论公式和柔度系数法，考虑安装误差、支撑误差的影响和接触柔度系数的非线性，求出齿轮的载荷分布；在此基础上，利用三维有限元模型进行受力计算得到可靠的弯曲应力，利用Hertz理论进行了接触应力的计算。

本发明解决其技术问题所采用的技术方案是:一种改进的轮齿加载接触分析方法，其特征在于包括以下步骤:

步骤1.推导及获得齿面方程、法向量；

根据微分几何与啮合原理，求出一对小大齿轮的齿面方程及法向量，其表示为:在与小齿轮固联的坐标系S₁下，小齿轮相邻三个轮齿的齿面方程及法向量为r₁₁(u₁,v₁),n₁₁(u₁,v₁),r₁₂(u₁,v₁)，n₁₂(u₁,v₁)和r₁₃(u₁,v₁)，n₁₃(u₁,v₁)；在与大齿轮固联的坐标系S₂下，与小齿轮各轮齿相啮合的大齿轮齿面的齿面方程及法向量为r₂₁(u₂,v₂),n₂₁(u₂,v₂),r₂₂(u₂,v₂),n₂₂(u₂,v₂)和r₂₃(u₂,v₂),n₂₃(u₂,v₂)；

步骤2.进行无载荷的轮齿接触分析；

①计算接触印痕和传动误差

以小齿轮三轮齿的中间齿为基准面，分别计算进行多对齿面的无载荷轮齿接触分析，得到在小齿轮取不同转角φ₁时，多对齿面副的接触印痕和传动误差，其具体如下:

建立一个装配坐标系S_f，分别将小齿轮、大齿轮的某对啮合轮齿的齿面方程r_1j(u₁,v₁),r_2j(u₂,v₂)及其法向量n_1j(u₁,v₁),n_2j(u₂,v₂)转换到装配坐标系S_f下，建立轮齿接触分析方程组；

$\left\{\begin{array}{c}{r}_{f1j}(u_1,v_1,{\mathrm{\φ}}_1)=r_{f2j}(u_2,v_2,{\mathrm{\φ}}_2)\\ n_{f1j}(u_1,v_1,{\mathrm{\φ}}_1)=n_{f2j}(u_2,v_2,{\mathrm{\φ}}_2)\end{array}\right.---\left(1\right)$

给定一系列的小齿轮转角φ₁，通过上述方程组可求出另外五个未知量:u₁,v₁,u₂,v₂和大齿轮转角φ₂；

将上述求得变量分别代入小大齿轮的齿面方程r_1j(u₁,v₁),r_2j(u₂,v₂)，则可得到某对齿面上一系列的接触点，即接触路径；

将(1)式求解参数代入传动误差公式则可得到接触点处的传动误差，并绘制不同轮齿接触对的传动误差曲线；

②推导齿面的主曲率、主方向

已知，一个齿面方程为r(u,v)，则单位法矢n＝(r_u×r_v)/|r_u×r_v|，式中， $r_v=\∂r/\∂v;$

曲面r(u,v)上某点P(x,y)的法曲率表示如下:

$k_n=\frac{L{(\mathrm{du}/\mathrm{dv})}^2+2M(\mathrm{du}/\mathrm{dv})+N}{E{(\mathrm{du}/\mathrm{dv})}^2+2F(\mathrm{du}/\mathrm{dv})+G}---\left(2\right)$

式中，E,F,G为曲面的第一基本量，L,M,N为曲面的第二基本量，du/dv为点P(x，y)处的切线方向，其中，F=r_u·r_v，L=-n·r_uu，M=n·r_uv，N＝n·r_vv；

由(2)式知，k_n随μ＝du/dv而变化，其最大值和最小值成为点P(x,y)的主曲率，对应切线方向即为主方向；

k_n是μ的函数，对μ求导得

$({\mathrm{E\μ}}^2+2\mathrm{F\μ}+G)\frac{{\mathrm{dk}}_n}{\mathrm{d\μ}}+2[(k_{n}E-L)\mathrm{\μ}+(k_{n}F-M)]=0---\left(3\right)$

当dk_n/dμ＝0时，k_n取极值，有

(k_nE-L)μ+(k_nF-M)＝0     (4)

联立式(2)与(4)得到

(k_nF-M)μ+(k_nG-N)＝0     (5)

联立式(4)与(5)，消除μ，则有

$(\mathrm{EG}-F^2)k_n^2-(\mathrm{EN}-2\mathrm{FM}+\mathrm{GL})k_n+(\mathrm{LN}-M^2)=0---\left(6\right)$

求解式(6)则可得到主曲率；

联立式(4)与(5)，消除k_n，则有

(EM-FL)μ²+(EN-GL)μ+(FN-GM)＝0     (7)

求解式(7)则可得到主方向；

根据罗德里克定理，判断与主曲率对应的主方向，其关系式为:

$k_{\mathrm{ni}}=-\frac{\mathrm{dn}}{\mathrm{dr}}=-\frac{n_u{\mathrm{\μ}}_i+n_v}{r_u{\mathrm{\μ}}_i+r_v}---\left(8\right)$

根据式(7)求出的主方向参数μ，可确定主方向的单位矢量:

$e_i=\frac{r_u\mathrm{du}+r_v\mathrm{dv}}{|r_u\mathrm{du}+r_v\mathrm{dv}|}=\frac{r_u{u}_i+r_v}{|r_u{u}_i+r_v|}---\left(9\right)$

③相切齿面的诱导法曲率，求得接触椭圆的长短轴方向

根据步骤②求得的主曲率和主方向可得到相切齿面的诱导法曲率，求得接触椭圆的长短轴方向，其过程如下:

规定小齿轮齿面上的第一主方向到大齿轮齿面上的第一主方向的有向角为

${\mathrm{\σ}}^{\left(12\right)}=\mathrm{sign}(e_I^{\left(1\right)}\×e_I^{\left(2\right)}\·n)\·\arccos (e_I^{\left(1\right)}\·e_I^{\left(2\right)})---\left(10\right)$

已知曲面的主曲率及主方向，根据法曲率的欧拉公式表达出任意方向的法曲率

$k_n^{\left(i\right)}=k_I^{\left(i\right)}{\cos }^2{\mathrm{\σ}}^{\left(i\right)}+k_{\mathrm{II}}^{\left(i\right)}{\sin }^2{\mathrm{\σ}}^{\left(i\right)}=1/2(k_{\mathrm{\Σ}}^{\left(i\right)}+g_i\cos 2{\mathrm{\σ}}^{\left(i\right)})---\left(11\right)$

式中， $k_{\mathrm{\Σ}}^{\left(i\right)}=k_I^{\left(i\right)}+k_{\mathrm{II}}^{\left(i\right)},$ $g_i=k_I^{\left(i\right)}-k_{\mathrm{II}}^{\left(i\right)},$ i＝1,2；

两曲面在切点沿同一方向的法曲率之差，即该方向的诱导法曲率，表示为

$k_n^{\left(12\right)}=k_n^{\left(1\right)}-k_n^{\left(2\right)}=1/2(k_{\mathrm{\Σ}}^{\left(1\right)}-k_{\mathrm{\Σ}}^{\left(2\right)}+g_1\cos 2{\mathrm{\σ}}^{\left(1\right)}-g_2\cos 2{\mathrm{\σ}}^{\left(2\right)})---\left(12\right)$

将σ⁽²⁾＝σ⁽¹⁾-σ⁽¹²⁾代入式(12)，化简得

$k_n^{\left(12\right)}=1/2[(k_{\mathrm{\Σ}}^{\left(1\right)}-k_{\mathrm{\Σ}}^{\left(2\right)})+(g_1-g_2\cos 2{\mathrm{\σ}}^{\left(12\right)})\cos 2{\mathrm{\σ}}^{\left(1\right)}-g_2\sin 2{\mathrm{\σ}}^{\left(12\right)}\sin 2{\mathrm{\σ}}^{\left(1\right)}]---\left(13\right)$

是σ⁽¹⁾的函数，对式(13)求导，令得到诱导法曲率取极值时对应的σ⁽¹⁾，其化简后表示为:

$\left\{\begin{array}{c}{\cos }^2{\mathrm{\σ}}^{\left(1\right)}=\frac{1}{2}(1+\frac{g_1-g_2\cos 2{\mathrm{\σ}}^{\left(12\right)}}{R_{12}})\\ {\sin }^2{\mathrm{\σ}}^{\left(1\right)}=\frac{1}{2}(1-\frac{g_1\cos 2{\mathrm{\σ}}^{\left(12\right)}-g_2}{R_{12}})\end{array}\right.---\left(14\right)$

式中， $R_{12}=\sqrt{g_1^2-{2g}_1{g}_2\cos 2{\mathrm{\σ}}^{\left(12\right)}+g_2^2};$

由此，可以进一步得到接触椭圆的长短轴方向:

$\left\{\begin{array}{c}\stackrel{\→}{\mathrm{\η}}=e_I^{\left(1\right)}\cos {\mathrm{\σ}}^{\left(1\right)}+e_{\mathrm{II}}^{\left(1\right)}\sin {\mathrm{\σ}}^{\left(1\right)}\\ \mathrm{\ζ}=-e_I^{\left(1\right)}\sin {\mathrm{\σ}}^{\left(1\right)}+e_{\mathrm{II}}^{\left(1\right)}\cos {\mathrm{\σ}}^{\left(1\right)}\end{array}\right.---\left(15\right)$

④建立单个接触位置上接触椭圆大小、齿面压缩量与载荷的关系

由上述步骤①、②、③，求得某时刻接触点位置、法向量、主曲率、主方向、接触椭圆长短轴方向，建立某时刻接触点位置上接触椭圆大小、齿面压缩量与载荷的关系，具体如下:

根据Hertz接触理论，可得到两轮齿压缩量δ_c与接触椭圆大小之间的关系

$\left\{\begin{array}{c}{\mathrm{\δ}}_c={\mathrm{Aa}}^2\\ {\mathrm{\δ}}_c={\mathrm{Bb}}^2\end{array}\right.---\left(16\right)$

式中 $A=1/4(k_{\mathrm{\Σ}}^{\left(1\right)}-k_{\mathrm{\Σ}}^{\left(2\right)}-R_{12}),$ $B=1/4(k_{\mathrm{\Σ}}^{\left(1\right)}-k_{\mathrm{\Σ}}^{\left(2\right)}+R_{12});$

同时，可得到两轮齿压缩量δ_c与载荷P_i的关系

${\mathrm{\δ}}_c={\left(\frac{{9P}_i^2}{{16E}^{*2}{R}_e}\right)}^{1/3}{F}_2\left(e\right)---\left(17\right)$

式中，v₁,v₂,E₁,E₂分别为小大齿轮的泊松比和弹性模量；R_e为等效半径，R_e＝(R′R″)^1/2，1/R′＝1/R₁′+1/R₂′，1/R″＝1/R₁″+1/R₂″，R_i′和R_i″分别为小大齿轮在椭圆中心位置的主曲率半径；F₂(e)函数可看作是对于椭圆偏心率的修正因子，其相对于相对曲率比值(R′/R″)的变化曲线已知；

接触椭圆大小、齿面压缩量与载荷的关系由式(16)、(17)确定；

步骤3.计算载荷分配；

通过步骤1，步骤2,得到某时刻小大齿轮副的多齿面接触情况，即某时刻各个齿面的接触位置、传动误差，同时建立单个接触位置上接触椭圆大小、齿面压缩量与载荷的关系，下面依据变形协调条件可求出各个接触位置上的载荷分量；

首先，根据轮齿无载荷接触分析得到的多对齿面的传动误差，可判断小齿轮在转动某转角φ₁时，小大齿轮的接触点对数m；

设各接触齿面对上分担的载荷分别为P₁,P₂,……P_m，则根据变形协调条件和总载荷不变得到方程组:

$\left\{\begin{array}{c}\mathrm{\Δ}{\mathrm{\θ}}_1\left(P_1\right)=\mathrm{\Δ}{\mathrm{\θ}}_2\left(P_2\right)=\·\·\·={\mathrm{\Δ\θ}}_i\left(P_i\right)\·\·\·={\mathrm{\Δ\θ}}_m\left(P_m\right)\\ R_{G1}{P}_1\cos {\mathrm{\α}}_{t1}+R_{G2}{P}_2\cos {\mathrm{\α}}_{t2}+\·\·\·\·\·\·+R_{\mathrm{Gm}}{P}_m\cos {\mathrm{\α}}_{\mathrm{tm}}=T\end{array}\right.---\left(18\right)$

式中，R_Gi为第i对齿轮接触点距转动中心的转动半径；α_ti表示载荷与切平面之间的夹角；Δθ_i为加载传动误差，T为小大齿轮传递总扭矩；

已知，加载传动误差Δθ_i由三部分组成，表示为

Δθ_i(P_i)＝Δθ_io(P_i)+Δθ_ic(P_i)+Δθ_ib(P_i)     (19)

式中，Δθ_io(P_i)为无加载传动误差，Δθ_ic(P_i)为接触变形引起的加载传动误差，Δθ_ib(P_i)为轮齿弯曲变形引起的加载传动误差；

下面对三种传动误差分量进行求解:

①求解Δθ_io(P_i),根据上述无载荷轮齿接触分析，可得到

${\mathrm{\Δ\θ}}_{\mathrm{io}}\left(P_i\right)=({\mathrm{\φ}}_{i2}-{\mathrm{\φ}}_{i2}^{*})-({\mathrm{\φ}}_{i1}-{\mathrm{\φ}}_{i1}^{*})\frac{z_1}{z_2}---\left(20\right)$

②求解Δθ_ic(P_i),上述按Hertz接触理论建立起两轮齿压缩量δ_ic与载荷P_i的关系即δ_ic＝f(P_i)，则有

${\mathrm{\Δ\θ}}_{\mathrm{ic}}\left(P_i\right)=\frac{{2\mathrm{\δ}}_{\mathrm{ic}}}{R_{\mathrm{Gi}}}---\left(21\right)$

③求解Δθ_ib(P_i),首先，分别建立大小齿轮的三齿三维有限元网格模型，进行齿轮各个齿面的弯曲柔度系数矩阵提取，其方法如下:

对齿轮的三齿三维有限元模型添加两侧面的位移约束和底面位移约束，下面，分别给工作齿面上的各个节点施加1N的法向力，提取齿面上所有节点的位移，若某个工作齿面上有n个节点，则提取出的柔度系数矩阵为[C_f1]_n×n×3；

对齿轮的三齿三维有限元模型添加两侧面的位移约束和底面位移约束，同时将各个非工作齿面添加位移约束，下面分别给工作齿面上的各个节点施加1N的法向力，提取齿面上所有节点的位移，若某个工作齿面上有n个节点，则提取出的柔度系数矩阵为[C_f2]_n×n×3；

则齿轮的弯曲柔度系数矩阵为

[C_f]＝[C_f1]-[C_f2]       (22)

根据Hertz接触理论可得，一般外形曲面的表面压力为

$p=\frac{3P}{2\mathrm{\πab}}\{1-{(x/a)}^2-{(y/b)}^2{\}}^{1/2}---\left(23\right)$

将接触椭圆区域离散化，建立分割区域内载荷对接触点的弯曲弹性变形量的微分表达式

dδ_bi＝p(ρ,θ)·c_f(ρ,θ)·ρdρdθ      (24)

由此，建立一个分割区域内载荷对接触点的弯曲弹性变形量的近似表达式

Δδ_bi＝p(ρ,θ)·c_f(ρ,θ)·ρ·Δρ·Δθ     (25)

其中，ρ,θ是通过将接触椭圆在极坐标下的表示，且有

$\left\{\begin{array}{c}x=a\·\mathrm{\ρ}\·\cos \mathrm{\θ}\\ y=b\·\mathrm{\ρ}\·\sin \mathrm{\θ}\end{array}\right.---\left(26\right)$

故ρ的积分区间为[0,1]，θ的积分区间为[0,2π]；c_f(ρ,θ)为分割区域类的柔度系数，通过对区域范围内或周围的节点柔度系数插值而言，插值选取节点越多，c_f(ρ,θ)越精确；

最后，对接触椭圆内所有分割区域内载荷对接触点的弯曲弹性变形量进行求和，则可得到一个齿面上由弯曲引起的接触点位置的位移量:

${\mathrm{\δ}}_{\mathrm{bi}}=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\Δ\δ}}_{\mathrm{bij}}---\left(27\right)$

由弯曲变形引起的传动误差分量为

${\mathrm{\Δ\θ}}_{\mathrm{bi}}\left(P_i\right)=\frac{{\mathrm{\δ}}_{\mathrm{bi}1}+{\mathrm{\δ}}_{\mathrm{bi}2}}{R_{\mathrm{Gi}}}---\left(28\right)$

式中，δ_bi1,δ_bi2分别代表小大齿轮一对齿轮副在接触点位置的轮齿弯曲变形量；

联立式(18)～(28)，进行求解，则可得到齿间载荷分布P₁,P₂,……P_m；

再次计算计入弯曲变形和接触变形后的传动误差，判断小齿轮在转动某转角φ₁时，小大齿轮的接触点对数k，若k>m，则将k赋予m，重新进行步骤3，若k＝m，则进行下一步；

步骤4.进行强度计算；

按照上述分析计算得到的齿间载荷分配结果，则可根据Hertz理论公式计算小大齿轮各对轮齿副在小齿轮转动某转角φ₁时的齿面压力分布和接触应力分布；

根据上述建立小、大齿轮三维有限元模型，将计算好的齿面压力分布等效施加在接触区域内的节点上，采用现有有限元软件进行弯曲强度计算，则可得到轮齿接触产生的弯曲应力分布。

有益效果

本发明提出的一种改进的轮齿加载接触分析方法，通过几何计算得到齿轮的齿面方程和法向量方程；在此基础上进行无载荷的轮齿接触分析，得到齿轮啮合的接触印痕、传动误差，齿面的主曲率、主方向，相切齿面的诱导法曲率，从而建立接触椭圆大小、齿面压缩量与载荷三者之间的关系；通过上述得到的传动误差、啮合印痕、以及载荷与接触椭圆大小、齿面压缩量之间的关系，联系变形协调方程、Hertz方程，则可求解一对齿轮副的齿间载荷分布；最终，根据齿间载荷分配结果，采用Hertz公式计算小大齿轮某时刻下的齿面压力分布和接触应力分布，采用现有有限元软件计算得到轮齿接触产生的弯曲应力分布。该方法考虑了安装误差、支撑误差的影响和接触柔度系数的非线性，能够精确和简便的计算轮齿接触应力和弯曲应力，对齿轮学的设计、计算和发展有重要意义。

附图说明

下面结合附图和实施方式对本发明一种改进的轮齿加载接触分析方法作进一步详细说明。

图1为相切齿面的诱导法曲率示意图。

图2a为非工作齿面无位移约束的有限元计算模型示意图。

图2b为添加非工作齿面位移约束的有限元计算模型示意图。

图3为接触椭圆区域的分割示意图。

具体实施方式

本实施例轮齿加载接触分析方法是以一对三齿几何模型进行计算分析。

第一步,推导及获得齿面方程、法向量；

根据微分几何与啮合原理，求出一对小大齿轮的齿面方程，其表示方式为:

小齿轮中间啮合齿面在与小齿轮固联的坐标系S₁下的齿面方程及法向量为r₁₁(u₁,v₁)和n₁₁(u₁,v₁)，根据坐标旋转变换可得到小齿轮中间齿左右轮齿的齿面方程及其法向量为r₁₂(u₁,v₁),n₁₂(u₁,v₁)和r₁₃(u₁,v₁),n₁₃(u₁,v₁)；

同理，在与大齿轮固联的坐标系S₂下，与小齿轮各轮齿相啮合的大齿轮齿面的齿面方程及法向量为r₂₁(u₂,v₂),n₂₁(u₂,v₂)、r₂₂(u₂,v₂),n₂₂(u₂,v₂)和r₂₃(u₂,v₂),n₂₃(u₂,v₂)；

第二步,进行无载荷的轮齿接触分析；

①计算接触印痕和传动误差

以小齿轮三轮齿的中间齿为基准轮齿，分别计算进行多对齿面的无载荷轮齿接触分析，得到在小齿轮取不同转角φ₁时，多对齿面副的接触印痕和传动误差，其具体如下:

建立一个装配坐标系S_f，分别将小、大齿轮的某对啮合轮齿的齿面方程r_1j(u₁,v₁),r_2j(u₂,v₂)及其法向量n_1j(u₁,v₁),n_2j(u₂,v₂)转换到装配坐标系S_f下，建立轮齿接触分析方程组；

$\left\{\begin{array}{c}{r}_{f1j}(u_1,v_1,{\mathrm{\φ}}_1)=r_{f2j}(u_2,v_2,{\mathrm{\φ}}_2)\\ n_{f1j}(u_1,v_1,{\mathrm{\φ}}_1)=n_{f2j}(u_2,v_2,{\mathrm{\φ}}_2)\end{array}\right.---\left(1\right)$

给定一系列的小齿轮转角φ₁，通过上述方程组可求出另外五个未知量:u₁,v₁,u₂,v₂和大齿轮转角φ₂；

将上述求得变量分别代入小大齿轮的齿面方程r_1j(u₁,v₁),r_2j(u₂,v₂)，则可得到某对齿面上一系列的接触点，即接触路径；将上述求得变量分别代入小大齿轮的法向量方程n_1j(u₁,v₁),n_2j(u₂,v₂)，则可得到某对齿面上对应接触点位置的法向量；

通过传动误差公式则可得到接触点处的传动误差，并绘制不同轮齿接触对的传动误差曲线；

②推导齿面的主曲率、主方向

已知，一个齿面方程表示为r(u,v)，则单位法矢n＝(r_u×r_v)/|r_u×r_v|，式中， $r_u=\∂r/\∂u,$ $r_v=\∂r/\∂v;$

曲面r(u,v)上某点P(x,y)的法曲率表示如下:

$k_n=\frac{L{(\mathrm{du}/\mathrm{dv})}^2+2M(\mathrm{du}/\mathrm{dv})+N}{E{(\mathrm{du}/\mathrm{dv})}^2+2F(\mathrm{du}/\mathrm{dv})+G}---\left(2\right)$

式中，E,F,G为曲面的第一基本量，L,M,N为曲面的第二基本量，du/dv为点P(x,y)处的切线方向；其中，F＝r_u·r_v,L＝-n·r_uu,M＝n·r_uv,N＝n·r_vv；

由(2)式知，k_n随μ＝du/dv而变化，其最大值和最小值成为点P(x,y)的主曲率，对应切线方向即为主方向；

k_n是μ的函数，对μ求导得

$({\mathrm{E\μ}}^2+2\mathrm{F\μ}+G)\frac{{\mathrm{dk}}_n}{\mathrm{d\μ}}+2[(k_{n}E-L)\mathrm{\μ}+(k_{n}F-M)]=0---\left(3\right)$

当dk_n/dμ＝0时，k_n取极值，有

(k_nE-L)μ+(k_nF-M)＝0     (4)

联立式(2)与(4)得到

(k_nF-M)μ+(k_nG-N)＝0     (5)

联立式(4)与(5)，消除μ，则有

$(\mathrm{EG}-F^2)k_n^2-(\mathrm{EN}-2\mathrm{FM}+\mathrm{GL})k_n+(\mathrm{LN}-M^2)=0---\left(6\right)$

求解式(6)则可得到主曲率；

联立式(4)与(5)，消除k_n，则有

(EM-FL)μ²+(EN-GL)μ+(FN-GM)＝0     (7)

求解式(7)则可得到主方向；

根据罗德里克定理，判断与主曲率对应的主方向，其关系式为：

$k_{\mathrm{ni}}=-\frac{\mathrm{dn}}{\mathrm{dr}}=-\frac{n_u{\mathrm{\μ}}_i+n_v}{r_u{\mathrm{\μ}}_i+r_v}---\left(8\right)$

根据式(7)求出的主方向参数μ，可确定主方向的单位矢量:

$e_i=\frac{r_u\mathrm{du}+r_v\mathrm{dv}}{|r_u\mathrm{du}+r_v\mathrm{dv}|}=\frac{r_u{u}_i+r_v}{|r_u{u}_i+r_v|}---\left(9\right)$

下面，用分别代表主方向的单位矢量、主曲率，上标(i)＝(1)、(2)分别代表小、大齿轮，下标j＝I、II分别代表第一主曲率、第二主曲率及其对应的主方向，下面计算并规定小齿轮的第一主曲率为法曲率的最大值，大齿轮的第一主曲率为法曲率的最小值；

③相切齿面的诱导法曲率，求得接触椭圆的长短轴方向

根据步骤②求得的主曲率和主方向可得到相切齿面的诱导法曲率，求得接触椭圆的长短轴方向，其过程如下:

参阅图1，Σ₁、Σ₂为空间相切的小、大齿面，切点为P，在切点处两齿面的单位法矢n的方向相同，σ⁽ⁱ⁾分别为主方向到任意方向的有向角；

规定小齿轮齿面上的第一主方向到大齿轮齿面上的第一主方向的有向角为

${\mathrm{\σ}}^{\left(12\right)}=\mathrm{sign}(e_I^{\left(1\right)}\×e_I^{\left(2\right)}\·n)\·\arccos (e_I^{\left(1\right)}\·e_I^{\left(2\right)})---\left(10\right)$

已知曲面的主曲率及主方向，根据法曲率的欧拉公式表达出任意方向的法曲率

$k_n^{\left(i\right)}=k_I^{\left(i\right)}{\cos }^2{\mathrm{\σ}}^{\left(i\right)}+k_{\mathrm{II}}^{\left(i\right)}{\sin }^2{\mathrm{\σ}}^{\left(i\right)}=1/2(k_{\mathrm{\Σ}}^{\left(i\right)}+g_i\cos 2{\mathrm{\σ}}^{\left(i\right)})---\left(11\right)$

式中， $k_{\mathrm{\Σ}}^{\left(i\right)}=k_I^{\left(i\right)}+k_{\mathrm{II}}^{\left(i\right)},$ $g_i=k_I^{\left(i\right)}-k_{\mathrm{II}}^{\left(i\right)},$ i＝1,2；

两曲面在切点沿同一方向的法曲率之差称为该方向的诱导法曲率，表示为

$k_n^{\left(12\right)}=k_n^{\left(1\right)}-k_n^{\left(2\right)}=1/2(k_{\mathrm{\Σ}}^{\left(1\right)}-k_{\mathrm{\Σ}}^{\left(2\right)}+g_1\cos 2{\mathrm{\σ}}^{\left(1\right)}-g_2\cos 2{\mathrm{\σ}}^{\left(2\right)})---\left(12\right)$

将σ⁽²⁾＝σ⁽¹⁾-σ⁽¹²⁾代入式(12)，化简得

$k_n^{\left(12\right)}=1/2[(k_{\mathrm{\Σ}}^{\left(1\right)}-k_{\mathrm{\Σ}}^{\left(2\right)})+(g_1-g_2\cos 2{\mathrm{\σ}}^{\left(12\right)})\cos 2{\mathrm{\σ}}^{\left(1\right)}-g_2\sin 2{\mathrm{\σ}}^{\left(12\right)}\sin 2{\mathrm{\σ}}^{\left(1\right)}]---\left(13\right)$

是σ⁽¹⁾的函数，对式(13)求导，令得到诱导法曲率取极值时对应的σ⁽¹⁾，其化简后表示为:

$\left\{\begin{array}{c}{\cos }^2{\mathrm{\σ}}^{\left(1\right)}=\frac{1}{2}(1+\frac{g_1-g_2\cos 2{\mathrm{\σ}}^{\left(12\right)}}{R_{12}})\\ {\sin }^2{\mathrm{\σ}}^{\left(1\right)}=\frac{1}{2}(1-\frac{g_1\cos 2{\mathrm{\σ}}^{\left(12\right)}-g_2}{R_{12}})\end{array}\right.---\left(14\right)$

式中， $R_{12}=\sqrt{g_1^2-{2g}_1{g}_2\cos 2{\mathrm{\σ}}^{\left(12\right)}+g_2^2};$

由此，可进一步得到接触椭圆的长短轴方向:

$\left\{\begin{array}{c}\stackrel{\→}{\mathrm{\η}}=e_I^{\left(1\right)}\cos {\mathrm{\σ}}^{\left(1\right)}+e_{\mathrm{II}}^{\left(1\right)}\sin {\mathrm{\σ}}^{\left(1\right)}\\ \mathrm{\ζ}=-e_I^{\left(1\right)}\sin {\mathrm{\σ}}^{\left(1\right)}+e_{\mathrm{II}}^{\left(1\right)}\cos {\mathrm{\σ}}^{\left(1\right)}\end{array}\right.---\left(15\right)$

④建立单个接触位置上接触椭圆大小、齿面压缩量与载荷的关系

由上述步骤①、②、③,求得某时刻接触点位置、法向量、主曲率、主方向、接触椭圆长短轴方向，建立某时刻接触点位置上接触椭圆大小、齿面压缩量与载荷的关系，具体如下:

根据Hertz接触理论，可得到两轮齿压缩量δ_c与接触椭圆大小之间的关系

$\left\{\begin{array}{c}{\mathrm{\δ}}_c={\mathrm{Aa}}^2\\ {\mathrm{\δ}}_c={\mathrm{Bb}}^2\end{array}\right.---\left(16\right)$

式中 $A=1/4(k_{\mathrm{\Σ}}^{\left(1\right)}-k_{\mathrm{\Σ}}^{\left(2\right)}-R_{12}),$ $B=1/4(k_{\mathrm{\Σ}}^{\left(1\right)}-k_{\mathrm{\Σ}}^{\left(2\right)}+R_{12});$

同时，可得到两轮齿压缩量δ_c与载荷P_i的关系

${\mathrm{\δ}}_c={\left(\frac{{9P}_i^2}{{16E}^{*2}{R}_e}\right)}^{1/3}{F}_2\left(e\right)---\left(17\right)$

式中，v₁,v₂,E₁,E₂分别为小大齿轮的泊松比和弹性模量；R_e为等效半径，R_e＝(R′R″)^1/2，1/R′＝1/R₁′+1/R₂′，1/R″＝1/R₁″+1/R₂″，R_i′和R_i″分别为小大齿轮在椭圆中心位置的主曲率半径；F₂(e)函数可作为是对于椭圆偏心率的修正因子，其相对于相对曲率比值(R′/R″)的变化曲线已知；

至此，接触椭圆大小、齿面压缩量与载荷的关系由式(16)、(17)确定；

第三步,计算载荷分配；

通过第一步,第二步,得到某时刻小大齿轮副的多齿面接触情况，即某时刻各个齿面的接触位置、传动误差，同时建立单个接触位置上接触椭圆大小、齿面压缩量与载荷的关系，下面依据变形协调条件则可求出各个接触位置上的载荷分量；

首先，根据上述轮齿无载荷接触分析得到多对齿面的传动误差，可判断小齿轮在转动某转角φ₁时，小大齿轮的接触点对数m；

假设，各接触齿面对上分担的载荷分别为P₁,P₂,……P_m，则根据变形协调条件和总载荷不变得到方程组:

$\left\{\begin{array}{c}\mathrm{\Δ}{\mathrm{\θ}}_1\left(P_1\right)=\mathrm{\Δ}{\mathrm{\θ}}_2\left(P_2\right)=\·\·\·={\mathrm{\Δ\θ}}_i\left(P_i\right)\·\·\·={\mathrm{\Δ\θ}}_m\left(P_m\right)\\ R_{G1}{P}_1\cos {\mathrm{\α}}_{t1}+R_{G2}{P}_2\cos {\mathrm{\α}}_{t2}+\·\·\·\·\·\·+R_{\mathrm{Gm}}{P}_m\cos {\mathrm{\α}}_{\mathrm{tm}}=T\end{array}\right.---\left(18\right)$

式中，R_Gi为第i对齿轮接触点距转动中心的转动半径，α_ti为载荷与切平面之间的夹角，Δθ_i为加载传动误差，T为小大齿轮传递总扭矩；

已知，加载传动误差Δθ_i由三部分组成，表示为

Δθ_i(P_i)＝Δθ_io(P_i)+Δθ_ic(P_i)+Δθ_ib(P_i)      (19)

式中，Δθ_io(P_i)为无加载传动误差，Δθ_ic(P_i)为接触变形引起的加载传动误差，Δθ_ib(P_i)为轮齿弯曲变形引起的加载传动误差；

下面对这三种传动误差分量进行求解:

①求解Δθ_io(P_i),根据上述无载荷轮齿接触分析，可得到

${\mathrm{\Δ\θ}}_{\mathrm{io}}\left(P_i\right)=({\mathrm{\φ}}_{i2}-{\mathrm{\φ}}_{i2}^{*})-({\mathrm{\φ}}_{i1}-{\mathrm{\φ}}_{i1}^{*})\frac{z_1}{z_2}---\left(20\right)$

②求解Δθ_ic(P_i),上述按Hertz接触理论建立起两轮齿压缩量δ_ic与载荷P_i的关系，即δ_ic＝f(P_i)，则有

${\mathrm{\Δ\θ}}_{\mathrm{ic}}\left(P_i\right)=\frac{{2\mathrm{\δ}}_{\mathrm{ic}}}{R_{\mathrm{Gi}}}---\left(21\right)$

③求解Δθ_ib(P_i),首先，分别建立大小齿轮的三齿三维有限元网格模型，进行齿轮各个齿面的弯曲柔度系数矩阵提取，其方法如下:

参阅图2a，对齿轮的三齿三维有限元模型添加两侧面的位移约束和底面位移约束，分别给工作齿面上的各个节点施加1N的法向力，提取齿面上所有节点的位移，若某个工作齿面上有n个节点，则提取出的柔度系数矩阵为[C_f1]_n×n×3；

参阅图2b，对齿轮的三齿三维有限元模型添加两侧面的位移约束和底面位移约束，同时将各个非工作齿面添加位移约束，分别给工作齿面上的各个节点施加1N的法向力，提取齿面上所有节点的位移，若某个工作齿面上有n个节点，则提取出的柔度系数矩阵为[C_f2]_n×n×3；

则齿轮的弯曲柔度系数矩阵为

[C_f]＝[C_f1]-[C_f2]      (22)

根据Hertz接触理论可得，一般外形曲面的表面压力为

$p=\frac{3P}{2\mathrm{\πab}}\{1-{(x/a)}^2-{(y/b)}^2{\}}^{1/2}---\left(23\right)$

参阅图3，将接触椭圆区域离散化，建立分割区域内载荷对接触点的弯曲弹性变形量的微分表达式

dδ_bi＝p(ρ,θ)·c_f(ρ,θ)·ρdρdθ      (24)

由此，建立一个分割区域内载荷对接触点的弯曲弹性变形量的近似表达式

Δδ_bi＝p(ρ,θ)·c_f(ρ,θ)·ρ·Δρ·Δθ     (25)

其中，ρ,θ是通过将接触椭圆在极坐标下的表示，且有

$\left\{\begin{array}{c}x=a\·\mathrm{\ρ}\·\cos \mathrm{\θ}\\ y=b\·\mathrm{\ρ}\·\sin \mathrm{\θ}\end{array}\right.---\left(26\right)$

故ρ的积分区间为[0,1]，θ的积分区间为[0,2π]；c_f(ρ,θ)为分割区域类的柔度系数，通过对区域范围内或周围的节点柔度系数插值而言，插值选取节点越多，c_f(ρ,θ)越精确；

最后，对接触椭圆内所有分割区域内载荷对接触点的弯曲弹性变形量进行求和，则可得到一个齿面上由弯曲引起的接触点位置的位移量:

${\mathrm{\δ}}_{\mathrm{bi}}=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\Δ\δ}}_{\mathrm{bij}}---\left(27\right)$

由弯曲变形引起的传动误差分量为

${\mathrm{\Δ\θ}}_{\mathrm{bi}}\left(P_i\right)=\frac{{\mathrm{\δ}}_{\mathrm{bi}1}+{\mathrm{\δ}}_{\mathrm{bi}2}}{R_{\mathrm{Gi}}}---\left(28\right)$

式中，δ_bi1,δ_bi2分别代表小大齿轮一对齿轮副在接触点位置的轮齿弯曲变形量；

联立式(18)～(28)，进行求解，则可得到齿间载荷分布P₁,P₂,……P_m；

再次计算计入弯曲变形和接触变形后的传动误差，判断小齿轮在转动某转角φ₁时，小大齿轮的接触点对数k，若k>m，则将k赋予m，重新进行第三步，若k＝m，则进行下一步；

第四步,进行强度计算；

按照上述分析计算得到的齿间载荷分配结果，则根据Hertz理论公式计算小大齿轮各对轮齿副在小齿轮转动某转角φ₁时的齿面压力分布和接触应力分布；

根据上述建立小、大齿轮三维有限元模型，将计算好的齿面压力分布等效施加在接触区域内或者含接触区域附近的节点上，采用现有有限元软件ANSYS或ABAQUS进行弯曲强度计算，则可分别得到小大轮齿承载接触产生的弯曲应力分布。

本实施例所计算轮齿弯曲应力和接触应力的方法，考虑安装误差、支撑误差的影响和接触柔度系数的非线性，同时又尽量保证接触区域的几何形状，能较好的获得齿轮在拟真实工况下的工作性能。

Claims (1)

1.一种改进的轮齿加载接触分析方法，其特征在于包括以下步骤:

步骤1.推导及获得齿面方程、法向量；

根据微分几何与啮合原理，求出一对小大齿轮的齿面方程及法向量，其表示为:在与小齿轮固联的坐标系S₁下，小齿轮相邻三个轮齿的齿面方程及法向量为r₁₁(u₁,v₁),n₁₁(u₁,v₁),r₁₂(u₁,v₁)，n₁₂(u₁,v₁)和r₁₃(u₁,v₁)，n₁₃(u₁,v₁)；在与大齿轮固联的坐标系S₂下，与小齿轮各轮齿相啮合的大齿轮齿面的齿面方程及法向量为r₂₁(u₂,v₂),n₂₁(u₂,v₂),r₂₂(u₂,v₂),n₂₂(u₂,v₂)和r₂₃(u₂,v₂),n₂₃(u₂,v₂)；

步骤2.进行无载荷的轮齿接触分析；

①计算接触印痕和传动误差

以小齿轮三轮齿的中间齿为基准面，分别计算进行多对齿面的无载荷轮齿接触分析，得到在小齿轮取不同转角φ₁时，多对齿面副的接触印痕和传动误差，其具体如下:

建立一个装配坐标系S_f，分别将小齿轮、大齿轮的某对啮合轮齿的齿面方程r_1j(u₁,v₁),r_2j(u₂,v₂)及其法向量n_1j(u₁,v₁),n_2j(u₂,v₂)转换到装配坐标系S_f下，建立轮齿接触分析方程组；

$\left\{\begin{array}{c}{r}_{f1j}(u_1,v_1,{\mathrm{\φ}}_1)=r_{f2j}(u_2,v_2,{\mathrm{\φ}}_2)\\ n_{f1j}(u_1,v_1,{\mathrm{\φ}}_1)=n_{f2j}(u_2,v_2,{\mathrm{\φ}}_2)\end{array}\right.---\left(1\right)$

给定一系列的小齿轮转角φ₁，通过上述方程组可求出另外五个未知量:u₁,v₁,u₂,v₂和大齿轮转角φ₂；

将上述求得变量分别代入小大齿轮的齿面方程r_1j(u₁,v₁),r_2j(u₂,v₂)，则可得到某对齿面上一系列的接触点，即接触路径；

将(1)式求解参数代入传动误差公式则可得到接触点处的传动误差，并绘制不同轮齿接触对的传动误差曲线；

②推导齿面的主曲率、主方向

已知，一个齿面方程为r(u,v)，则单位法矢n＝(r_u×r_v)/|r_u×r_v|，式中， $r_v=\∂r/\∂v;$

曲面r(u,v)上某点P(x,y)的法曲率表示如下:

$k_n=\frac{{L(\mathrm{du}/\mathrm{dv})}^2+2M(\mathrm{du}/\mathrm{dv})+N}{{E(\mathrm{du}/\mathrm{dv})}^2+2F(\mathrm{du}/\mathrm{dv})+G}---\left(2\right)$

式中，E,F,G为曲面的第一基本量，L,M,N为曲面的第二基本量，du/dv为点P(x,y)处的切线方向,其中，F＝r_u·r_v,L＝-n·r_uu,M＝n·r_uv,N＝n·r_vv；

由(2)式知，k_n随μ＝du/dv而变化，其最大值和最小值成为点P(x,y)的主曲率，对应切线方向即为主方向；

k_n是μ的函数，对μ求导得

$(E{\mathrm{\μ}}^2+2\mathrm{F\μ}+G)\frac{{\mathrm{dk}}_n}{\mathrm{d\μ}}+2[(k_{n}E-L)\mathrm{\μ}+(k_{n}F-M)]=0---\left(3\right)$

当dk_n/dμ＝0时，k_n取极值，有

(k_nE-L)μ+(k_nF-M)＝0       (4)

联立式(2)与(4)得到

(k_nF-M)μ+(k_nG-N)＝0        (5)

联立式(4)与(5)，消除μ，则有

$(\mathrm{EG}-F^2)k_n^2-(\mathrm{EN}-2\mathrm{FM}+\mathrm{GL})k_n+(\mathrm{LN}-M^2)=0---\left(6\right)$

求解式(6)则可得到主曲率；

联立式(4)与(5)，消除k_n，则有

(EM-FL)μ²+(EN-GL)μ+(FN-GM)＝0        (7)

求解式(7)则可得到主方向；

根据罗德里克定理，判断与主曲率对应的主方向，其关系式为:

$k_{\mathrm{ni}}=-\frac{\mathrm{dn}}{\mathrm{dr}}=-\frac{n_u{\mathrm{\μ}}_i+n_v}{r_u{\mathrm{\μ}}_i+r_v}---\left(8\right)$

根据式(7)求出的主方向参数μ，可确定主方向的单位矢量:

$e_i=\frac{r_u\mathrm{du}+r_v\mathrm{dv}}{|r_u\mathrm{du}+r_v\mathrm{dv}|}=\frac{r_u{u}_i+r_v}{|r_u{u}_i+r_v|}---\left(9\right)$

③相切齿面的诱导法曲率，求得接触椭圆的长短轴方向

根据步骤②求得的主曲率和主方向可得到相切齿面的诱导法曲率，求得接触椭圆的长短轴方向，其过程如下:

规定小齿轮齿面上的第一主方向到大齿轮齿面上的第一主方向的有向角为

${\mathrm{\σ}}^{\left(12\right)}=\mathrm{sign}(e_I^{\left(1\right)}\×e_I^{\left(2\right)}\·n)\·\arccos (e_I^{\left(1\right)}\·e_I^{\left(2\right)})---\left(10\right)$

已知曲面的主曲率及主方向，根据法曲率的欧拉公式表达出任意方向的法曲率

$k_n^{\left(i\right)}=k_I^{\left(i\right)}{\cos }^2{\mathrm{\σ}}^{\left(i\right)}+k_{\mathrm{II}}^{\left(i\right)}{\sin }^2{\mathrm{\σ}}^{\left(i\right)}=1/2(k_{\mathrm{\Σ}}^{\left(i\right)}+g_i\cos {2\mathrm{\σ}}^{\left(i\right)})---\left(11\right)$

式中， $k_{\mathrm{\Σ}}^{\left(i\right)}=k_I^{\left(i\right)}+k_{\mathrm{II}}^{\left(i\right)},g_i=k_I^{\left(i\right)}-k_{\mathrm{II}}^{\left(i\right)},$ i＝1,2；

两曲面在切点沿同一方向的法曲率之差，即该方向的诱导法曲率，表示为

$k_n^{\left(12\right)}=k_n^{\left(1\right)}-k_n^{\left(2\right)}=1/2(k_{\mathrm{\Σ}}^{\left(1\right)}-k_{\mathrm{\Σ}}^{\left(2\right)}+g_i\cos {2\mathrm{\σ}}^{\left(1\right)}-g_2\cos {2\mathrm{\σ}}^{\left(2\right)})---\left(12\right)$

将σ⁽²⁾＝σ⁽¹⁾-σ⁽¹²⁾代入式(12)，化简得

$k_n^{\left(12\right)}=1/2[(k_{\mathrm{\Σ}}^{\left(1\right)}-k_{\mathrm{\Σ}}^{\left(2\right)})+(g_1-g_2\cos 2{\mathrm{\σ}}^{\left(12\right)})\cos 2{\mathrm{\σ}}^{\left(1\right)}-g_2\sin 2{\mathrm{\σ}}^{\left(12\right)}\sin 2{\mathrm{\σ}}^{\left(1\right)}]---\left(13\right)$

是σ⁽¹⁾的函数，对式(13)求导，令得到诱导法曲率取极值时对应的σ⁽¹⁾，其化简后表示为:

$\left\{\begin{array}{c}{\cos }^2{\mathrm{\σ}}^{\left(1\right)}=\frac{1}{2}(1+\frac{g_1-g_2\cos {2\mathrm{\σ}}^{\left(12\right)}}{R_{12}})\\ {\sin }^2{\mathrm{\σ}}^{\left(1\right)}=\frac{1}{2}(1-\frac{g_1\cos 2{\mathrm{\σ}}^{\left(12\right)}-g_2}{R_{12}})\end{array}\right.---\left(14\right)$

式中， $R_{12}=\sqrt{g_1^2-2g_1{g}_2\cos {2\mathrm{\σ}}^{\left(12\right)}+g_2^2};$

由此，可以进一步得到接触椭圆的长短轴方向:

$\left\{\begin{array}{c}\stackrel{\→}{\mathrm{\η}}=e_I^{\left(1\right)}\cos {\mathrm{\σ}}^{\left(1\right)}+e_{\mathrm{II}}^{\left(1\right)}\sin {\mathrm{\σ}}^{\left(1\right)}\\ \mathrm{\ζ}=-e_I^{\left(1\right)}\sin {\mathrm{\σ}}^{\left(1\right)}+e_{\mathrm{II}}^{\left(1\right)}\cos {\mathrm{\σ}}^{\left(1\right)}\end{array}\right.---\left(15\right)$

④建立单个接触位置上接触椭圆大小、齿面压缩量与载荷的关系

由上述步骤①、②、③，求得某时刻接触点位置、法向量、主曲率、主方向、接触椭圆长短轴方向，建立某时刻接触点位置上接触椭圆大小、齿面压缩量与载荷的关系，具体如下:

根据Hertz接触理论，可得到两轮齿压缩量δ_c与接触椭圆大小之间的关系

$\left\{\begin{array}{c}{\mathrm{\δ}}_c={\mathrm{Aa}}^2\\ {\mathrm{\δ}}_c=B{b}^2\end{array}\right.---\left(6\right)$

式中 $A=1/4(k_{\mathrm{\Σ}}^{\left(1\right)}-k_{\mathrm{\Σ}}^{\left(2\right)}-R_{12}),B=1/4(k_{\mathrm{\Σ}}^{\left(1\right)}-k_{\mathrm{\Σ}}^{\left(2\right)}+R_{12});$

同时，可得到两轮齿压缩量δ_c与载荷P_i的关系

${\mathrm{\δ}}_c={\left(\frac{9{P_i}^2}{16E^{*2}{R}_e}\right)}^{1/3}{F}_2\left(e\right)---\left(17\right)$

式中，v₁,v₂,E₁,E₂分别为小大齿轮的泊松比和弹性模量；R_e为等效半径，R_e＝(R′R″)^1/2，1/R′＝1/R₁′+1/R₂′，1/R″＝1/R₁″+1/R₂″，R_i′和R_i″分别为小大齿轮在椭圆中心位置的主曲率半径；F₂(e)函数可看作是对于椭圆偏心率的修正因子，其相对于相对曲率比值(R′R″)的变化曲线已知；

接触椭圆大小、齿面压缩量与载荷的关系由式(16)、(17)确定；

步骤3.计算载荷分配；

通过步骤1，步骤2,得到某时刻小大齿轮副的多齿面接触情况，即某时刻各个齿面的接触位置、传动误差，同时建立单个接触位置上接触椭圆大小、齿面压缩量与载荷的关系，下面依据变形协调条件可求出各个接触位置上的载荷分量；

首先，根据轮齿无载荷接触分析得到的多对齿面的传动误差，可判断小齿轮在转动某转角φ₁时，小大齿轮的接触点对数m；

设各接触齿面对上分担的载荷分别为P₁,P₂,……P_m，则根据变形协调条件和总载荷不变得到方程组:

$\left\{\begin{array}{c}\mathrm{\Δ}{\mathrm{\θ}}_1\left(P_1\right)=\mathrm{\Δ}{\mathrm{\θ}}_2\left(P_2\right)=\·\·\·=\mathrm{\Δ}{\mathrm{\θ}}_i\left(P_i\right)\·\·\·=\mathrm{\Δ}{\mathrm{\θ}}_m\left(P_m\right)\\ R_{G1}{P}_1\cos {\mathrm{\α}}_{t1}+R_{G2}{P}_2\cos {\mathrm{\α}}_{t2}+\·\·\·\·\·\·+R_{\mathrm{Gm}}{P}_m\cos {\mathrm{\α}}_{\mathrm{tm}}=T\end{array}\right.----\left(18\right)$

式中，R_Gi为第i对齿轮接触点距转动中心的转动半径；α_ti表示载荷与切平面之间的夹角；△θ_i为加载传动误差，T为小大齿轮传递总扭矩；

已知，加载传动误差△θ_i由三部分组成，表示为

△θ_i(P_i)＝△θ_io(P_i)+△θ_ic(P_i)+△θ_ib(P_i)   (19)

式中，△θ_io(P_i)为无加载传动误差，△θ_ic(P_i)为接触变形引起的加载传动误差，△θ_ib(P_i)为轮齿弯曲变形引起的加载传动误差；

下面对三种传动误差分量进行求解:

①求解△θ_io(P_i),根据上述无载荷轮齿接触分析，可得到

$\mathrm{\Δ}{\mathrm{\θ}}_{\mathrm{io}}\left(P_i\right)=({\mathrm{\φ}}_{i2}-{\mathrm{\φ}}_{i2}^{*})-{(\mathrm{\φ}}_{i1}-{\mathrm{\φ}}_{i1}^{*})\frac{z_1}{z_2}---\left(20\right)$

②求解△θ_ic(P_i),上述按Hertz接触理论建立起两轮齿压缩量δ_ic与载荷P_i的关系即δ_ic＝f(P_i)，则有

$\mathrm{\Δ}{\mathrm{\θ}}_{\mathrm{ic}}\left(P_i\right)=\frac{2{\mathrm{\δ}}_{\mathrm{ic}}}{R_{\mathrm{Gi}}}---\left(21\right)$

③求解△θ_ib(P_i),首先，分别建立大小齿轮的三齿三维有限元网格模型，进行齿轮各个齿面的弯曲柔度系数矩阵提取，其方法如下:

对齿轮的三齿三维有限元模型添加两侧面的位移约束和底面位移约束，下面，分别给工作齿面上的各个节点施加1N的法向力，提取齿面上所有节点的位移，若某个工作齿面上有n个节点，则提取出的柔度系数矩阵为[C_f1]_n×n×3；

对齿轮的三齿三维有限元模型添加两侧面的位移约束和底面位移约束，同时将各个非工作齿面添加位移约束，下面分别给工作齿面上的各个节点施加1N的法向力，提取齿面上所有节点的位移，若某个工作齿面上有n个节点，则提取出的柔度系数矩阵为[C_f2]_n×n×3；

则齿轮的弯曲柔度系数矩阵为

[C_f]＝[C_f1]-[C_f2]        (22)

根据Hertz接触理论可得，一般外形曲面的表面压力为

$p=\frac{3P}{2\mathrm{\πab}}{\{1-{(x/a)}^2-{(y/b)}^2\}}^{1/2}---\left(23\right)$

将接触椭圆区域离散化，建立分割区域内载荷对接触点的弯曲弹性变形量的微分表达式

dδ_bi＝p(ρ,θ)·c_f(ρ,θ)·ρdρdθ   (24)

由此，建立一个分割区域内载荷对接触点的弯曲弹性变形量的近似表达式

△δ_bi＝p(ρ,θ)·c_f(ρ,θ)·ρ·△ρ·△θ   (25)

其中，ρ,θ是通过将接触椭圆在极坐标下的表示，且有

$\left\{\begin{array}{c}x=a\·\mathrm{\ρ}\·\cos \mathrm{\θ}\\ y=b\·\mathrm{\ρ}\·\sin \mathrm{\θ}\end{array}\right.---\left(26\right)$

故ρ的积分区间为[0,1]，θ的积分区间为[0,2π]；c_f(ρ,θ)为分割区域类的柔度系数，通过对区域范围内或周围的节点柔度系数插值而言，插值选取节点越多，c_f(ρ,θ)越精确；

最后，对接触椭圆内所有分割区域内载荷对接触点的弯曲弹性变形量进行求和，则可得到一个齿面上由弯曲引起的接触点位置的位移量:

${\mathrm{\δ}}_{\mathrm{bi}}=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\Δ}{\mathrm{\δ}}_{\mathrm{bij}}---\left(27\right)$

由弯曲变形引起的传动误差分量为

$\mathrm{\Δ}{\mathrm{\θ}}_{\mathrm{bi}}\left(P_i\right)=\frac{{\mathrm{\δ}}_{\mathrm{bi}1}+{\mathrm{\δ}}_{\mathrm{bi}2}}{R_{\mathrm{Gi}}}---\left(28\right)$

式中，δ_bi1,δ_bi2分别代表小大齿轮一对齿轮副在接触点位置的轮齿弯曲变形量；

联立式(18)～(28)，进行求解，则可得到齿间载荷分布P₁,P₂,……P_m；

再次计算计入弯曲变形和接触变形后的传动误差，判断小齿轮在转动某转角φ₁时，小大齿轮的接触点对数k，若k>m，则将k赋予m，重新进行步骤3，若k＝m，则进行下一步；

步骤4.进行强度计算；

按照上述分析计算得到的齿间载荷分配结果，则可根据Hertz理论公式计算小大齿轮各对轮齿副在小齿轮转动某转角φ₁时的齿面压力分布和接触应力分布；

根据上述建立小、大齿轮三维有限元模型，将计算好的齿面压力分布等效施加在接触区域内的节点上，采用现有有限元软件进行弯曲强度计算，则可得到轮齿接触产生的弯曲应力分布。