齿轮接触有限元分析.pdf

第20卷第2期2003年4月 计 算 力 学 学 报 Chinese Journal of ComputationalM echan icsVol.20,No.2April 2003文章编号:100724708(2003)0220189206齿轮接触有限元分析杨生华(煤炭科学研究总院 上海分院,上海200030)摘要:通过接触仿真分析研究了通用接触单元在轮齿变形和接触应力计算中的应用。建立了一对齿轮接触仿真分析的模型,并使用新的接触单元法计算了轮齿变形和接触应力,与赫兹理论比较,同时也计算了摩擦力对接触应力的影响。计算分析了单元离散、几何、边界范围与加载或约束处理方式的误差,建立了一个计算轮齿变形和接触应力的标准,说明了新的接触单元法的精确性、有效性和可靠性。关键词:接触单元;轮齿变形;接触应力;计算标准;仿真分析中图分类号:TP391文献标识码:A收稿日期:2001204228;修改稿收到日期:20022062241基金项目:上海自然科学基金资助项目 1作者简介:杨生华(19632),男,硕士生,工程师 11引言计算接触非线性问题有许多方法,例如罚函数法、拉格朗日乘子法等,其中罚函数法由于其经济和方便而得到广泛使用。过去使用点2点接触单元,求解接触问题,对于象齿轮类接触,模型构造很麻烦,计算结果精度和准确性很难保证。随着计算机和有限元法的发展,新的接触单元法产生精确的几何模型,自动划分网格,自适应求解。新的单元计算精度更高,更有效,功能更强大。其中接触单元能非常有效地求解接触非线性问题,新的通用接触单元(包括点2面和面2面单元)特别适合于计算齿轮接触问题。在微机上能实现齿轮接触仿真分析,大大地促进了齿轮CA E的形成和发展。轮齿变形的有限元分析20世纪70年代已开始,但仅仅计算挠曲变形。接触变形和接触应力的有限元分析在20世纪90年代才真正开始。总之,过去的计算是基于试验的计算方法,计算方法是简化的、近似的,不够精确更不够可靠;没有使用有限元法研究轮齿接触变形和应力,并说明与赫兹变形和应力之间的差别,没有分析计算误差,没有考虑齿轮本体变形对轮齿变形的影响,没有计算摩擦力对接触应力的影响。文中使用AN SYS大型通用有限元分析软件,在个人计算机上建立齿轮接触仿真分析模型。通过两圆柱赫兹接触变形和应力验证其有效性和精度,分析计算了一对直齿轮的轮齿变形和接触应力,说明了新的接触单元法的精确性、有效性和可靠性。建立了一个计算轮齿变形和接触应力的标准或基准,给力学研究和机械设计人员一个参考。2通用接触单元的赫兹计算为了检验通用接触单元的有效性和精确性,赫兹计算验证是必要的。两无限长圆柱有限元计算网格模型如图1所示。结构单元是具有附加形状函数的四节点等参单元(一次单元)。图中接触处网格边长为二十分之一接触半宽,该模型节点为7444,单元为7280(其中接触单元为80个点2面单元)。计算参数和结果如表1所示,理论结果按公式(1)2(4)计算1。计算结果表明:有限元计算结果和理论计算结果一致,两圆柱变形计算误差仅分别为0.08%和0.045%。注意到公式(2)、(4)是按赫兹接触半无限空间推导的公式,因而是理论近似的(变形误差为1.7%、0.6%,应力误差为0.6%、0.4%),在接触点不远处一点的变形和应力与有限元计算结果基本一致,有限元计算结果略大于公式(2)和(4)与理论一致1。1995-2006 Tsinghua Tongfang Optical Disc Co.,Ltd.All rights reserved.表1两无限长圆柱接触分析Tab.1Two cylinder contact analysis参数半径距离理论计算有限元计算R(mm)d(mm)H(m)dH(m)max(N?mm2)H(m)dH(m)max(N?mm2)圆柱138.0737.1986.1513.898229.5786.1463.966230.85圆柱263.6706.0056.7103.701229.5786.7073.722230.28 注:弹性模量:E1=E2=2.06105N?mm2,泊松比:1=2=0.3,kn=10E1表2 有摩擦接触应力分析(单位:N?mm2)Tab.2Contact stress analysis w ith friction(U nit:N?mm2)摩擦系数00(0.001)1(10)20.1(0.001)1(42)20.2(0.05)1(93)20.2(0.001)1(14)2计算应力赫兹应力max(误差)3误差max(误差)3误差max(误差)3误差max(误差)3误差圆柱1229.578231.07(3.88)1.0%234.58(4.04)2.2%245.19(3.32)6.8%245.55(0.42)7.0%圆柱2229.578230.40(3.93)0.5%233.82(4.15)1.8%244.53(3.22)6.5%244.93(0.36)6.7%注:上标1为力的计算收敛误差,上标2为迭代次数,上标3为误差估计4 罚参数大小与计算效率和精度有关,罚参数越小计算误差越大2,但罚参数太大计算效率降低,而且由于单元离散本身有误差,计算精度不会有明显提高。因此罚参数有最佳范围,通常取1210倍接触体的弹性模量。网格密度也与计算效率和精度有关,网格越密计算精度越高而效率降低。使用一次单元时摩擦力使得计算效率明显减低,需要更多的迭代次数,摩擦系数越大效率和精度越低。表2是摩擦力对接触应力影响计算,计算模型网格密度接触处只有图1中一半,但需要另一半的对称模型。无摩擦时,迭代次数为10,摩擦系数为0.1时,力的收敛误差为0.001,迭代次数为42;摩擦系数0.2时,收敛误差为0.05,迭代次数为93。摩擦系数大于0.2时,计算难于进行。而使用同样网格二次8节点等参单元和面 面接触单元,能有效计算有摩擦接触问题。当摩擦系数为0.2时,收敛误差为0.001,迭代次数为14。a2=4PRE3(1)p0=2Pa其中1E3=1-v21E1+1-v22E1,1R=1R1+1R2diH=P(1-v2i)E32ln(2dia)-vi1-vi(2)H=P(1-v2i)E32ln(4Ria)-1(3)max=0.3003p0(4)a是赫兹接触半宽,p0是最大赫兹压力,H是接触变形,max是最大赫兹应力,d、R如图2。3 轮齿变形和接触应力的计算实例311 齿轮参数表1是计算轮齿变形和接触应力的一对齿轮参数,所有参数通过精确计算,根据功率和转速计算出单齿啮合时额定载荷沿齿向分布力大小为F=386.5N?mm。表3 计算一对齿形参数Tab.3Parameter of a pair of gear摸数压力角齿数变位系数顶圆半径齿数变位系数顶圆半径齿宽中心距功率转速mZ1X1ra1Z2X2ra2bAP(kw)n820290.1672125.25400.1687169.2634278.6060400 注:E=2.06105N?mm2,=0.3,刀具圆角半径=0.38m,刀具齿顶高hao=1.25m.091计算力学学报第20卷 1995-2006 Tsinghua Tongfang Optical Disc Co.,Ltd.All rights reserved.312 轮齿变形和接触应力的有限元计算模型轮齿变形包括挠曲变形和接触变形及基础变形:轮齿挠曲变形计算模型的边界范围通常取一齿宽,计算的变形是轮齿对称中心点G的载荷方向(齿面法向)的变形;轮齿接触变形是载荷作用点至轮齿对称点之间的变形;轮齿基础变形为轮齿根部的弹性倾斜对轮齿变形的影响。为了计算方便,把基础变形包括在挠曲变形中。轮齿挠曲变形的单齿有限元模型如图2所示,边界范围为PQRS(图中为二齿宽,轮缘厚度1.5m)。过去轮齿接触变形用赫兹接触理论公式近似计算,但轮齿接触变形和赫兹接触变形之间存在多大误差,考虑轮齿基础变形影响的轮齿挠曲变形有限元计算的边界范围应该取多大,需要建立一对啮合齿轮接触有限元仿真分析模型,图3是三齿啮合接触计算模型。接触分析还能计算接触应力和应力分布,并能考虑摩擦力的影响,计算齿轮接触应力与赫兹接触应力之间的误差。313 计算模型的网格对齿轮接触分析来说,为了有效地生成网格,模型划分为接触区域,接触轮齿和非接触轮齿三个部分,接触单元最后再产生。由于自动生成的接触单元较多,需要控制接触面和目标面范围,接触范围一般不超过两倍的赫兹接触长度。图4是三齿接触一个模型网格,为了得到精确的变形,接触轮齿的相邻轮齿的网格也需要适当加密。为了得到精确的接触应力,接触处网格要更密些,通常单元边长为赫兹半宽的十分之一或更小。图4中接触处单元边长为赫兹半宽的十分之一,接触区半径为赫兹半宽的1.5倍,图中右下角为接触区网格放大图,该模型总节点数5632,单元数5325,其中接触单元为60个。对整轮接触仿真模型,接触模型的其余部分(即两齿轮本体)可用超单元表示3。314 轮齿变形分开计算和仿真分析结果比较轮齿变形的计算方法有两种:一种是分开计算,即轮齿的挠曲变形按图2模型计算,接触变形按公式(2)计算;第二种方法,建立一对齿轮的啮合接触仿真分析模型,进行接触分析而得出轮齿变形。表4为分开计算时和仿真分析计算结果,仿真分析时三齿接触网格模型如图4,整轮啮合接触计算模型两轮本体都为实心,两轮本体内圆直径(轴径)都为90mm。分开计算时轮齿挠曲变形单齿模型的边界范围PQRS相对两个仿真模型分别取二齿宽和三齿宽。仿真计算结果表明:接触变形按赫兹变形公式计算有误差,由于齿轮接触已经是非赫兹接触,按照公式(2)计算有-7%左右误差。单齿挠曲变形计算的误差来源于边界范围,轮齿挠曲变形边界范围PS、QR应在23齿宽之间,只要适当调整可以和仿真分析取得一致结果。轮齿变形受到齿轮本体变形的影响,局部和整轮仿真分析结果误差达-9.1%。315 轮齿接触应力仿真分析结果与赫兹应力计算比较轮齿接触应力计算方法也有两种:赫兹接触应力公式计算和有限元接触仿真分析计算。由于齿轮是渐开线轮齿接触,赫兹接触应力肯定是近似的,特别在有摩擦时,必然存在误差,而接触仿真分析能计算其误差大小。表5是有无摩擦接触时整轮仿真分析计算结果和赫兹接触应力比较,齿轮啮合时由于摩擦力造成接触力增加。由表中看出:齿轮实际接触应力比赫兹接触应力大,均超过5%。当摩擦系数从0提高到0.2时,赫兹接触应力误差达10%,而齿轮接触应力也提高5%以上,当载荷增加时赫兹接触应力误差也增加,3倍载荷时达10%,而且接触应力分布计算结果最大应力深度大于赫兹理论0.786a(a为赫兹接触半宽)。191第2期杨生华:齿轮接触有限元分析 1995-2006 Tsinghua Tongfang Optical Disc Co.,Ltd.All rights reserved.表4 一对齿轮变形分析比较(单位:m)Tab.4Comparing the tooth deformation analysis of a pair of gear(U nit:m)计算方法和误差轮1挠曲变形接触变形轮2挠曲变形接触变形总变形分开计算仿真分析误差%分开计算仿真分析误差%三齿模型整轮模型4.3823.8987.0473.70118.9913.7924.1747.0883.99119.045(4.130)14.1767.0883.993(19.387)114.8(-2.2)1-6.3-1.4-6.9-0.4(-3.7)15.1353.8987.8183.70120.5524.9134.1767.8863.99220.9674.5-6.7-0.9-7.3-2.0 注:上标1为力分布作用在轮1的三个边界上,如图3所示。表5 有摩擦接触仿真分析和赫兹接触应力比较(单位:N?mm2)Tab.5Comparing the contact si m ulation analysis w ith friction w ith Hertz stress(U nit:N?mm2)摩擦系数00.1(0.001)10.2(0.05)10(3倍载荷)计算应力max赫兹应力误差max赫兹应力误差max赫兹应力误差max赫兹应力误差轮1242.05229.585.4%248.18233.706.2%259.23238.048.9%433.68397.699.1%轮2245.01229.586.7%254.39233.708.9%260.63238.049.5%438.91397.6910.4%注:上标1为接触力收敛误差4轮齿变形和接触应力的计算 误差分析不考虑物理模型本身的误差,包括计算参数不确定性和随机性,轮齿变形和接触应力计算模型的误差主要包括单元本身、单元离散、几何、边界范围和计算模型处理方法、二维与三维之间的误差。由于个人计算机现在都是32位,正向64位转变,数值计算误差通常很小,一般能达5位精度以上。单元本身、单元离散、几何误差可以使用更多的或更精确的单元、更大的接触刚度3、更精确的几何实体模型解决。在轮齿变形和接触应力计算中,这些误差现在可以肯定地减少到1%以下,而边界范围、二维与三维之间对轮齿变形的误差较大,需要建立仿真分析模型,即二维整轮接触仿真和三维整轮仿真分析模型,通过仿真分析计算误差也能达到1%以内。三维计算需要更多的花费,而且对计算机有更高的要求,通常不进行,三维模型产生的误差参考文献2。一次等参单元模型结构时由于刚度变大5,计算变形略小,赫兹变形实际计算误差可小于1(参考表1)。赫兹应力误差较大,表2中有1%。表2中的误差估计4说明了它的计算精度和可靠性。表6是不同接触刚度和网格及边界约束条件时 三齿接触仿真分析模型计算结果,由接触刚度和几何误差产生的变形计算误差都很小(约0.2左右),三次样条曲线拟合渐开线(齿侧用18节点三次样条拟合)使计算变形略变小。接触处网格加密一倍计算结果与原来的比较相差0.9。使用不同边界约束(不同加载方式)产生的误差为1.8%。表7是使用一次和二次单元仿真分析计算变形的比较,一次单元和单轮仿真计算变形的误差。由表7可知表4中轮1挠曲变形有0.8计算误差,这个误差是因为轮1旋转非线性产生。单轮仿真时计算轮齿挠曲变形使用集中载荷有0.1%的误差,单轮仿真网格与接触仿真网格相同,由p单元计算出单轮仿真网格离散误差。齿轮接触应力计算使用三对齿接触模型(图4)已足够,接触区网格边长为赫兹半宽的二十分之一(在赫兹长度上划分40个单元),如图4所示,接触处几何误差小于6.5微米。子模型是以接触点为原点,半径为1毫米的圆,网格密度为原网格的两倍,而粗网格仅为原网格密度一半。二次单元使用面接触单元,网格与一次单元相同。表8为不同加载和几何精度及计算方法下的接触应力,周向和切向载荷分别为作用在轮1底边,周向为轮1底边切线方向,切向是沿啮合线方向。对轮1旋转载荷非线性及几何非线性也进行了计算,几何精度按渐开线修正和子模型方法2计算得到更精确的值。表8计算表明,通用接触单元和其它单元一样对小变形分析对计算结果影响很小,子模型方法能与二次单元计算结果比较,所有计算模型相对误差都在1%内。根据表1赫兹应力计算比较,齿轮接触应力有限元计算结果误差估计在1%内。291计算力学学报第20卷 1995-2006 Tsinghua Tongfang Optical Disc Co.,Ltd.All rights reserved.表6三对齿轮接触仿真变形分析(计算轮1齿顶干涉量)(单位:m)Tab.6Contact si mulation analysis of three pair of teeth(U nit:m)加载位置轮1变形对称点接触点干涉点轮2变形对称点接触点干涉点总变形干涉量底边115.99011.81619.7827.80811.7990.72019.045底边215.98311.80719.7757.80811.8030.72019.049底边315.99411.81819.7907.81611.8110.72119.062底边415.99311.81719.7857.80811.8010.72019.048左右底边15.99311.81720.1247.80811.8010.72019.387注:上标1是接触刚度为10e,上标2是接触刚度为100e,上标3是更密网格,4是渐开线精确几何接触。表7轮齿变形仿真分析一次单元和二次单元比较(单位:m)Tab.7Comparing si mulation analyses of the tooth deformation(U nit:m)计算变形接触变形一次单元二次单元误差挠曲变形(接触仿真)一次单元(误差)二次单元挠曲变形(单轮仿真)二次单元(误差1)p单元(误差1)轮1变形4.174044.182390.20%4.90931(0.4%)4.927824.92847(0.13)4.92854(0.15)轮2变形3.991453.998790.18%7.88586(0.5%)7.925337.93385(0.1%)7.94308(0.2%)注:上标1与接触仿真比表8几个计算模型的接触应力比较(单位:N?mm2)Tab.8Comparison of the contact stresses of severalmodel computations(U nit:N?mm2)计算模型一次单元整轮仿真周向载荷切向载荷旋转渐开线子模型粗网格Kn=100E二次单元周向载荷齿轮1244.17241.74241.68242.25241.89241.66242.92241.95241.92齿轮2244.10246.42246.51246.94246.59243.68243.72246.64243.75通过接触仿真分析,接触单元法和计算技术能使轮齿变形和接触应力的计算误差从过去10%减少到1%以内,而且人为误差大大减少,模型的可靠性大大提高。通过随机和不确定分析将进一步提高模型的可靠性和分析精度,分析和仿真真实世界。5结语本文检验了通用接触单元的有效性。通过接触仿真分析,定义了轮齿变形(挠曲、接触和基础变形)和接触应力的有限元计算方法,说明了轮齿接触变形和应力传统赫兹理论存在的误差和齿轮本体变形对轮齿变形与摩擦力对轮齿接触应力的影响,证明了本文计算模型的有效性。分析了轮齿变形和接触应力计算的误差来源,建立了计算轮齿变形和接触应力的标准。为齿轮修缘和接触强度评价提供更可靠、更有价值的信息,齿轮接触过程的仿真分析,将为齿轮的动态设计、优化设计和可靠性设计打下新的基础。这样不仅能优化齿轮结构、齿形和齿廓,还能优化齿轮材料和工艺,能实现齿轮结构、材料和工艺的创新设计。随着计算机和有限元软件的发展,新的更有效、更精确、更通用的接触单元将继续产生。二次接触单元能更精确更可靠计算接触应力和变形,更有效地计算摩擦接触问题。求解方法多样化,使用罚函数法和拉格朗日乘子法两种方法求解具有实际意义,能使求解接触问题可靠性更高和稳健性更强,且无需与赫兹理论比较。齿轮变形和应力的仿真分析是发展的必然趋势,仿真分析进入三维领域。分析内容除了静态外,还包括动态和啮合过程。计算模型将更真实、更精确、更全面,通过个人图形工作站,既能快速计算又能更加直观、仔细、迅速、精确地观察到计算结果。误差可控制在1%内,不仅是实验法无法相比而且也是实验法无法做到的,为齿轮CA E奠定了基础。参考文献(References):1Johnson K L.接触力学M.徐秉业,等译.北京:高等教育出版社,1992.(Johnson K L.Contact M e2chanicsM.Cambridge U niversity Press,1985.(inChinese)2杨生华.子结构、子模型方法在齿轮应力和变形计算中的应用A.工程与科学中的计算力学C.北京:北京大学出版社,2001:7012707.(Yang Shenghua.Applications of Substructuring and Submodeling in391第2期杨生华:齿轮接触有限元分析 1995-2006 Tsinghua Tongfang Optical Disc Co.,Ltd.All rights reserved.Calculations of the Gear Stress and DeformationA.Computational M echanics in Engineering andScienceC.Beijing U niversity Press,2001:7012707.(in Chinese)3Barlam D,Zahavi E.The reliability of solutions incontact problem s J.Comp&S truct,1999,70:35245.4Zienkiew iczO C,Zhu J Z.A si mple error esti matorand adaptive procedure for practical engineering anal2ysis J.Int N um ber M ethods Eng,1987,24:3372357.5库克R D.有限元分析的概念和应用M.何穷,程耿东,译.北京:科学出版社,1981.(Cook R D.Concep ts and A pp lications of F inite E lem ent A nalysisM.JohnW iley,1974.(in Chinese)Fin ite element analysis of gear contactYang Shenghua(ShanghaiBranch of China Coal Research Institute,Shanghai 200030,China)Abstract:By the contact si mulation analysis the applicationsof the general contact elements in the calcu2lations of the tooth deformation and contact stress of gears have been investigated.The tooth contactmodel of a pair of mesh gear has been created,new finite element methods are used for calculating thetooth deformation and contact stress of gear,comparing w ith the hertz theory,and the contact stressw ith friction has been calculated.The errors in mesh discretization,geometry,boundary area and load2ing methods have been analyzed by calculation,the computation benchmark of the tooth deformation andcontact stress is created,the accurateness,effectiveness and reliability of new finite elementmethods aredescribed.Key words:general contact element;toth deformation;contact stress;computation benchmark;si mula2tion analysis.491计算力学学报第20卷 1995-2006 Tsinghua Tongfang Optical Disc Co.,Ltd.All rights reserved.