[发明专利]一种改进的轮齿加载接触分析方法

摘要

本发明公开了一种改进的轮齿加载接触分析方法。该方法基于变形协调条件，通过几何计算得到齿轮的齿面方程和法向量方程；并在此基础上进行无载荷的轮齿接触分析，得到齿轮啮合的接触印痕、传动误差，齿面的主曲率、主方向，相切齿面的诱导法曲率，从而建立起接触椭圆大小、齿面压缩量与载荷三者之间的关系；结合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是μ的函数，对μ求导得$(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理论公式计算小大齿轮各对轮齿副在小齿轮转动某转角φ₁时的齿面压力分布和接触应力分布；根据上述建立小、大齿轮三维有限元模型，将计算好的齿面压力分布等效施加在接触区域内的节点上，采用现有有限元软件进行弯曲强度计算，则可得到轮齿接触产生的弯曲应力分布。