齿轮接触有限元分析_杨生华

第20卷第2期2003年4月

　计算力学学报　

C hinese Journal of Computational Mechanics

V ol. 20, N o. 2April 2003

文章编号:1007-4708(2003) 02-0189-06

齿轮接触有限元分析

杨生华

(煤炭科学研究总院上海分院, 上海200030)

摘　要:通过接触仿真分析研究了通用接触单元在轮齿变形和接触应力计算中的应用。建立了一对齿轮接触仿真分析的模型, 并使用新的接触单元法计算了轮齿变形和接触应力, 与赫兹理论比较, 同时也计算了摩擦力对接触应力的影响。计算分析了单元离散、几何、边界范围与加载或约束处理方式的误差, 建立了一个计算轮齿变形和接触应力的标准, 说明了新的接触单元法的精确性、有效性和可靠性。关键词:接触单元; 轮齿变形; 接触应力; 计算标准; 仿真分析中图分类号:T P 391　　　文献标识码:A

1　引　言

计算接触非线性问题有许多方法, 例如罚函数法、拉格朗日乘子法等, 其中罚函数法由于其经济和方便而得到广泛使用。过去使用点-点接触单元, 求解接触问题, 对于象齿轮类接触, 模型构造很麻烦, 计算结果精度和准确性很难保证。随着计算机和有限元法的发展, 新的接触单元法产生精确的几何模型, 自动划分网格, 自适应求解。新的单元计算精度更高, 更有效, 功能更强大。其中接触单元能非常有效地求解接触非线性问题, 新的通用接触单元(包括点-面和面-面单元) 特别适合于计算齿轮接触问题。在微机上能实现齿轮接触仿真分析, 大大地促进了齿轮C AE 的形成和发展。

轮齿变形的有限元分析20世纪70年代已开始, 但仅仅计算挠曲变形。接触变形和接触应力的有限元分析在20世纪90年代才真正开始。总之, 过去的计算是基于试验的计算方法, 计算方法是简化的、近似的, 不够精确更不够可靠; 没有使用有限元法研究轮齿接触变形和应力, 并说明与赫兹变形和应力之间的差别, 没有分析计算误差, 没有考虑齿轮本体变形对轮齿变形的影响, 没有计算摩擦力对接触应力的影响。

文中使用AN SYS 大型通用有限元分析软件, 在个人计算机上建立齿轮接触仿真分析模型。通过两圆柱赫兹接触变形和应力验证其有效性和精度, 分析计算了一对直齿轮的轮齿变形和接触应力, 说

收稿日期:2001-04-28; 修改稿收到日期:2002-06-24.

基金项目:上海自然科学基金资助项目.

作者简介:杨生华(1963-) , 男, 硕士生, 工程师.

明了新的接触单元法的精确性、有效性和可靠性。建立了一个计算轮齿变形和接触应力的标准或基准, 给力学研究和机械设计人员一个参考。

2　通用接触单元的赫兹计算

为了检验通用接触单元的有效性和精确性, 赫兹计算验证是必要的。两无限长圆柱有限元计算网格模型如图1所示。结构单元是具有附加形状函数的四节点等参单元(一次单元) 。图中接触处网格边长为二十分之一接触半宽, 该模型节点为7444, 单元为7280(其中接触单元为80个点-面单元) 。计算参数和结果如表1所示, 理论结果按公式(1) -(4) 计算[1]。计算结果表明:有限元计算结果和理论计算结果一致, 两圆柱变形计算误差仅分别为0. 08%和0. 045%。注意到公式(2) 、(4) 是按赫兹接触半无限空间推导的公式, 因而是理论近似的(变形误差为1. 7%、0. 6%, 应力误差为0. 6%、0. 4%) , 在接触点不远处一点的变形和应力与有限元计算结果基本一致, 有限元计算结果略大于公式(2) 和(4) 与理论一致[1]。

190

计算力学学报

　第20卷　

表1　两无限长圆柱接触分析Tab. 1　Tw o cylinder co ntact a nalysis

参数圆柱1圆柱2

半径R (mm) 38. 07363. 670

距离d (mm ) 7. 1986. 005

W H (m) 6. 1516. 710

理论计算

d

W H (

有限元计算

2f max (N /mm)

m) W H (m) 6. 1466. 707

d

W H (m) 2f ma x (N /mm )

3. 8983. 701

229. 578229. 578

3. 9663. 722

230. 85230. 28

　　　　注:弹性模量:E 1=E 2=2. 06×105N /mm2, 泊松比:ν1=ν2=0. 3, k n =10E 1

表2　有摩擦接触应力分析(单位:N /mm)

Tab. 2　Contact stress a nalysis with friction (Unit :N /mm2)

摩擦系数计算应力圆柱1圆柱2

0赫兹应力229. 578229. 578

0(0. 001) 1(10) 2

3f 误差max (误差) 231. 07(3. 88) 230. 40(3. 93)

1. 0%0. 5%

0. 1(0. 001) 1(42) 2

3f 误差max (误差) 234. 58(4. 04) 233. 82(4. 15)

2. 2%1. 8%

0. 2(0. 05) 1(93) 2

3f 误差max (误差) 245. 19(3. 32) 244. 53(3. 22)

6. 8%6. 5%

0. 2(0. 001) 1(14) 2

3f 误差m a x (误差) 245. 55(0. 42) 244. 93(0. 36)

7. 0%6. 7%

2

　　注:上标1为力的计算收敛误差, 上标2为迭代次数, 上标3为误差估计[4]

　　罚参数大小与计算效率和精度有关, 罚参数越小计算误差越大[2], 但罚参数太大计算效率降低, 而且由于单元离散本身有误差, 计算精度不会有明显提高。因此罚参数有最佳范围, 通常取1-10倍接触体的弹性模量。网格密度也与计算效率和精度有关, 网格越密计算精度越高而效率降低。使用一次单元时摩擦力使得计算效率明显减低, 需要更多的迭代次数, 摩擦系数越大效率和精度越低。表2是摩擦力对接触应力影响计算, 计算模型网格密度接触处只有图1中一半, 但需要另一半的对称模型。无摩擦时, 迭代次数为10, 摩擦系数为0. 1时, 力的收敛误差为0. 001, 迭代次数为42; 摩擦系数0. 2时, 收敛误差为0. 05, 迭代次数为93。摩擦系数大于0. 2时, 计算难于进行。而使用同样网格二次8节点等参单元和面—面接触单元, 能有效计算有摩擦接触问题。当摩擦系数为0. 2时, 收敛误差为0. 001, 迭代次数为14。

a 2=

πE *

(1)

2i i i

　　W=P 2ln(a ) -1-v i

πE

d i H

(2) (3) (4)

i i

W H =P 2ln() -*

πE a

f max =0. 3003p 0

2

a 是赫兹接触半宽, p 0是最大赫兹压力, W H 是接触变形, f max 是最大赫兹应力, d 、R 如图2。

3　轮齿变形和接触应力的计算实例

3. 1　齿轮参数

表1是计算轮齿变形和接触应力的一对齿轮参数, 所有参数通过精确计算, 根据功率和转速计算出单齿啮合时额定载荷沿齿向分布力大小为F =386. 5N /mm 。

p 0=πa

其中

22

=1+2, 　=+E 1E 1R R 1R 2E

表3　计算一对齿形参数Tab. 3　Parameter of a pair o f gear

摸数m 8

压力角T 20°

齿数Z 129

变位系数

X 10. 1672

顶圆半径

r a 1125. 25

齿数Z 240

变位系数

X 20. 1687

顶圆半径

r a 2169. 26

齿宽b 34

中心距A 278. 60

功率P (kw ) 60

转速n 400

　　注:E =2. 06×105N /mm2, ν=0. 3, 刀具圆角半径d =0. 38m, 刀具齿顶高h ao =1. 25m.

　第2期

杨生华:齿轮接触有限元分析

191

3. 2　轮齿变形和接触应力的有限元计算模型

轮齿变形包括挠曲变形和接触变形及基础变形:轮齿挠曲变形计算模型的边界范围通常取一齿宽, 计算的变形是轮齿对称中心点G 的载荷方向(齿面法向) 的变形; 轮齿接触变形是载荷作用点至轮齿对称点之间的变形; 轮齿基础变形为轮齿根部的弹性倾斜对轮齿变形的影响。为了计算方便, 把基础变形包括在挠曲变形中。轮齿挠曲变形的单齿有限元模型如图2所示, 边界范围为PQRS(图中为二齿宽, 轮缘厚度1. 5m ) 。

过去轮齿接触变形用赫兹接触理论公式近似计算, 但轮齿接触变形和赫兹接触变形之间存在多大误差, 考虑轮齿基础变形影响的轮齿挠曲变形有限元计算的边界范围应该取多大, 需要建立一对啮合齿轮接触有限元仿真分析模型, 图3是三齿啮合接触计算模型。接触分析还能计算接触应力和应力分布, 并能考虑摩擦力的影响, 计算齿轮接触应力与赫兹接触应力之间的误差。3. 3　计算模型的网格

对齿轮接触分析来说, 为了有效地生成网格, 模型划分为接触区域, 接触轮齿和非接触轮齿三个部分, 接触单元最后再产生。由于自动生成的接触单元较多, 需要控制接触面和目标面范围, 接触范围一般不超过两倍的赫兹接触长度。图4是三齿接触一个模型网格, 为了得到精确的变形, 接触轮齿的相邻轮齿的网格也需要适当加密。为了得到精确的接触应力, 接触处网格要更密些, 通常单元边长为赫兹半宽的十分之一或更小。图4中接触处单元边长为赫兹半宽的十分之一, 接触区半径为赫兹半宽的1. 5倍, 图中右下角为接触区网格放大图, 该模型总节点数5632, 单元数5325, 其中接触单元为60个。对整轮接触仿真模型, 接触模型的其余部分(即两齿轮本体) 可用超单元表示。

[3]

3. 4　轮齿变形分开计算和仿真分析结果比较

轮齿变形的计算方法有两种:一种是分开计算, 即轮齿的挠曲变形按图2模型计算, 接触变形

按公式(2) 计算; 第二种方法, 建立一对齿轮的啮合接触仿真分析模型, 进行接触分析而得出轮齿变形。

表4为分开计算时和仿真分析计算结果, 仿真分析时三齿接触网格模型如图4, 整轮啮合接触计算模型两轮本体都为实心, 两轮本体内圆直径(轴径) 都为90mm 。分开计算时轮齿挠曲变形单齿模型的边界范围PQ RS 相对两个仿真模型分别取二齿宽和三齿宽。

仿真计算结果表明:接触变形按赫兹变形公式计算有误差, 由于齿轮接触已经是非赫兹接触, 按照公式(2) 计算有-7%左右误差。单齿挠曲变形计算的误差来源于边界范围, 轮齿挠曲变形边界范围PS 、QR 应在2～3齿宽之间, 只要适当调整可以和仿真分析取得一致结果。轮齿变形受到齿轮本体变形的影响, 局部和整轮仿真分析结果误差达-9. 1%。

3. 5　轮齿接触应力仿真分析结果与赫兹

应力计算比较

轮齿接触应力计算方法也有两种:赫兹接触应力公式计算和有限元接触仿真分析计算。由于齿轮是渐开线轮齿接触, 赫兹接触应力肯定是近似的, 特别在有摩擦时, 必然存在误差, 而接触仿真分析能计算其误差大小。表5是有无摩擦接触时整轮仿真分析计算结果和赫兹接触应力比较, 齿轮啮合时由于摩擦力造成接触力增加。由表中看出:齿轮实际接触应力比赫兹接触应力大, 均超过5%。当摩擦系数从0提高到0. 2时, 赫兹接触应力误差达10%, 而齿轮接触应力也提高5%以上, 当载荷增加时赫兹接触应力误差也增加, 3倍载荷时达10%, 而且接触应力分布计算结果最大应力深度大于赫兹理论0. 786a(a 为赫兹接触半宽) 。

192

计算力学学报

　第20卷　

表4　一对齿轮变形分析比较(单位:μm )

Tab. 4　Co mpa ring the too th defo rmation analy sis of a pair o f gear(Unit:μm )

计算方法和误差

分开计算仿真分析误差%分开计算仿真分析误差%

整轮模型三齿模型

轮1

挠曲变形4. 3823. 792(4. 130) 114. 8(-2. 2) 1

5. 1354. 9134. 5

接触变形3. 8984. 1744. 176-6. 33. 8984. 176-6. 7

挠曲变形7. 0477. 0887. 088-1. 47. 8187. 886-0. 9

轮2

接触变形3. 7013. 9913. 993-6. 93. 7013. 992-7. 3

总变形18. 99119. 045(19. 387) 1-0. 4(-3. 7) 1

20. 55220. 967-2. 0

　　注:上标1为力分布作用在轮1的三个边界上, 如图3所示。

表5　有摩擦接触仿真分析和赫兹接触应力比较(单位:N /mm2)

Tab. 5　Comparing the co ntact simulatio n analy sis with frictio n with Hertz stress (Unit :N /mm2)

摩擦系数计算应力轮1轮2

f max 242. 05245. 01

0赫兹应力229. 58229. 58

误差5. 4%6. 7%

f ma x 248. 18254. 39

0. 1(0. 001) 1赫兹应力233. 70233. 70

误差6. 2%8. 9%

f max 259. 23260. 63

0. 2(0. 05) 1赫兹应力238. 04238. 04

误差8. 9%9. 5%

f max 433. 68438. 91

0(3倍载荷) 赫兹应力397. 69397. 69

误差9. 1%10. 4%

　　注:上标1为接触力收敛误差

4　轮齿变形和接触应力的计算

　误差分析

　　不考虑物理模型本身的误差, 包括计算参数不确定性和随机性, 轮齿变形和接触应力计算模型的误差主要包括单元本身、单元离散、几何、边界范围和计算模型处理方法、二维与三维之间的误差。由于个人计算机现在都是32位, 正向64位转变, 数值计算误差通常很小, 一般能达5位精度以上。

单元本身、单元离散、几何误差可以使用更多的或更精确的单元、更大的接触刚度、更精确的几何实体模型解决。在轮齿变形和接触应力计算中, 这些误差现在可以肯定地减少到1%以下, 而边界范围、二维与三维之间对轮齿变形的误差较大, 需要建立仿真分析模型, 即二维整轮接触仿真和三维整轮仿真分析模型, 通过仿真分析计算误差也能达到1%以内。三维计算需要更多的花费, 而且对计算机有更高的要求, 通常不进行, 三维模型产生的误差参考文献[2]。

一次等参单元模型结构时由于刚度变大[5], 计算变形略小, 赫兹变形实际计算误差可小于1‰(参考表1) 。赫兹应力误差较大, 表2中有1%。表2中的误差估计说明了它的计算精度和可靠性。表6是不同接触刚度和网格及边界约束条件时三齿接触仿真分析模型计算结果, 由接触刚度和几何误差产生的变形计算误差都很小(约0. 2‰左右) , [4]

[3]

条拟合) 使计算变形略变小。接触处网格加密一倍计算结果与原来的比较相差0. 9‰。使用不同边界约束(不同加载方式) 产生的误差为1. 8%。

表7是使用一次和二次单元仿真分析计算变形的比较, 一次单元和单轮仿真计算变形的误差。由表7可知表4中轮1挠曲变形有0. 8‰计算误差, 这个误差是因为轮1旋转非线性产生。单轮仿真时计算轮齿挠曲变形使用集中载荷有0. 1%的误差, 单轮仿真网格与接触仿真网格相同, 由p 单元计算出单轮仿真网格离散误差。

齿轮接触应力计算使用三对齿接触模型(图4) 已足够, 接触区网格边长为赫兹半宽的二十分之一(在赫兹长度上划分40个单元) , 如图4所示, 接触处几何误差小于6. 5微米。子模型是以接触点为原点, 半径为1毫米的圆, 网格密度为原网格的两倍, 而粗网格仅为原网格密度一半。二次单元使用面接触单元, 网格与一次单元相同。表8为不同加载和几何精度及计算方法下的接触应力, 周向和切向载荷分别为作用在轮1底边, 周向为轮1底边切线方向, 切向是沿啮合线方向。对轮1旋转载荷非线性及几何非线性也进行了计算, 几何精度按渐开线修正和子模型方法计算得到更精确的值。表8计算表明, 通用接触单元和其它单元一样对小变形分析对计算结果影响很小, 子模型方法能与二次单元计算结果比较, 所有计算模型相对误差都在1%内。根据表1赫兹应力计算比较, 齿轮接触应力有1。

[2]

　第2期

杨生华:齿轮接触有限元分析

表6　三对齿轮接触仿真变形分析(计算轮1齿顶干涉量) (单位:μm ) Ta b. 6　Contact sim ula tion analysis o f three pair of teeth (Unit :μm)

轮1变形

对称点15. 99015. 98315. 99415. 99315. 993

接触点11. 81611. 80711. 81811. 81711. 817

干涉点19. 78219. 77519. 79019. 78520. 124

对称点7. 8087. 8087. 8167. 8087. 808

轮2变形接触点11. 79911. 80311. 81111. 80111. 801

干涉点0. 7200. 7200. 7210. 7200. 720

193

加载位置底边1底边2底边3底边

4

总变形干涉量

19. 04519. 04919. 06219. 04819. 387

左右底边

　　　　注:上标1是接触刚度为10e , 上标2是接触刚度为100e , 上标3是更密网格, 4是渐开线精确几何接触。

表7　轮齿变形仿真分析一次单元和二次单元比较(单位:μm )

Tab. 7　Comparing sim ula tion analyses of the too th defo rmatio n (Unit :μm )

计算变形轮1变形轮2变形

接触变形

一次单元4. 174043. 99145

二次单元4. 182393. 99879

误差0. 20%0. 18%

挠曲变形(接触仿真) 一次单元(误差) 4. 90931(0. 4%) 7. 88586(0. 5%)

二次单元4. 927827. 92533

挠曲变形(单轮仿真)

二次单元(误差1) 4. 92847(0. 13‰) 7. 93385(0. 1%)

p 单元(误差1) 4. 92854(0. 15‰) 7. 94308(0. 2%)

　　注:上标1与接触仿真比

表8　几个计算模型的接触应力比较(单位:N /mm )

　Com pariso n of the co ntact stresses of sev eral mo del com puta tions (Unit :N /mm 2) Tab . 8

计算模型齿轮1齿轮2

一次单元

整轮仿真244. 17244. 10

周向载荷241. 74246. 42

切向载荷241. 68246. 51

旋转242. 25246. 94

渐开线241. 89246. 59

子模型241. 66243. 68

粗网格242. 92243. 72

K n =100E

241. 95246. 64

二次单元周向载荷241. 92243. 75

2

　　通过接触仿真分析, 接触单元法和计算技术能使轮齿变形和接触应力的计算误差从过去10%减少到1%以内, 而且人为误差大大减少, 模型的可靠性大大提高。通过随机和不确定分析将进一步提高模型的可靠性和分析精度, 分析和仿真真实世界。

效、更精确、更通用的接触单元将继续产生。二次接触单元能更精确更可靠计算接触应力和变形, 更有效地计算摩擦接触问题。求解方法多样化, 使用罚函数法和拉格朗日乘子法两种方法求解具有实际意义, 能使求解接触问题可靠性更高和稳健性更强, 且无需与赫兹理论比较。齿轮变形和应力的仿真分析是发展的必然趋势, 仿真分析进入三维领域。分析内容除了静态外, 还包括动态和啮合过程。计算模型将更真实、更精确、更全面, 通过个人图形工作站, 既能快速计算又能更加直观、仔细、迅速、精确地观察到计算结果。误差可控制在1%内, 不仅是实验法无法相比而且也是实验法无法做到的, 为齿轮C AE 奠定了基础。

5　结　语

本文检验了通用接触单元的有效性。通过接触仿真分析, 定义了轮齿变形(挠曲、接触和基础变形) 和接触应力的有限元计算方法, 说明了轮齿接触变形和应力传统赫兹理论存在的误差和齿轮本体变形对轮齿变形与摩擦力对轮齿接触应力的影响, 证明了本文计算模型的有效性。分析了轮齿变形和接触应力计算的误差来源, 建立了计算轮齿变形和接触应力的标准。为齿轮修缘和接触强度评价提供更可靠、更有价值的信息, 齿轮接触过程的仿真分析, 将为齿轮的动态设计、优化设计和可靠性设计打下新的基础。这样不仅能优化齿轮结构、齿形和齿廓, 还能优化齿轮材料和工艺, 能实现齿轮结构、材料和工艺的创新设计。

, 参考文献(References ):

[1]　J o hnso n K L. 接触力学[M ].徐秉业, 等译. 北京:高

等教育出版社, 1992. (J oh nso n K L. Contact Me -chanics [M].Cambridge Univ ersity Press, 1985. (in Chinese ) )

[2]　杨生华. 子结构、子模型方法在齿轮应力和变形计算

中的应用[A ].工程与科学中的计算力学[C ].北京:北京大学出版社, 2001:701-707. (Yang Shenghua. a in

194

计算力学学报

　第20卷　

Calculations o f the Gear Stress and Defo rma tio n [A ].Co mputa tio nal M echanics in Engineering and :701-707. Science [C ].Beijing U niv er sity Press , 2001(inChinese) )

[3]　Bar lam D, Zahavi E. The reliability o f so lutions in

co nta ct pr oblem s [J].Comp &Struct , 1999, 70:35-45.

[4]　Zienkiewicz O C , Zhu J Z . A sim ple er ro r estimato r

and adaptiv e pr ocedure fo r practica l eng ineering ana l-:337-ysis [J ].Int N umber Methods Eng , 1987, 24357.

[5]　库克R D. 有限元分析的概念和应用[M].何　穷,

程耿东, 译. 北京:科学出版社, 1981. (Co ok R D. Concepts and Applications of Finite Element Analysis [M].John Wiley, 1974. (in Chinese) )

Finite element analysis of gear contact

Ya ng Shenghua

(Sha ng hai Bra nch of China Coa l Resear ch Institute , Shanghai 200030, China )

:By the contact simulatio n analysis the a pplicatio ns of the g eneral co ntact elem ents in the calcu-Abstract

lations of th e too th deforma tion and co ntact stress o f g ears hav e been inv estig ated. The too th co ntact m odel of a pair of mesh g ea r has been crea ted, new finite elem ent m ethods are used for calcula ting the too th deforma tio n and co ntact stress of g ear , com paring with the hertz theo ry , a nd the co ntact stress w ith friction has been calcula ted. The errors in mesh discretization, g eom etry , boundary a rea and load-ing methods have been a naly zed by calculatio n, the com puta tion benchmark of the to oth defo rmatio n and contact stress is crea ted , the accurateness , effectiveness and reliability of new finite elem ent methods are described.

Key words :g eneral co ntact elem ent ; toth defo rmatio n ; co ntact stress ; computation benchmark ; sim ula-tion analysis .