CN109145484B - 基于双曲面壳单元模型的数值载荷齿面接触分析方法 - Google Patents
基于双曲面壳单元模型的数值载荷齿面接触分析方法 Download PDFInfo
- Publication number
- CN109145484B CN109145484B CN201811026119.0A CN201811026119A CN109145484B CN 109145484 B CN109145484 B CN 109145484B CN 201811026119 A CN201811026119 A CN 201811026119A CN 109145484 B CN109145484 B CN 109145484B
- Authority
- CN
- China
- Prior art keywords
- contact
- tooth
- tooth surface
- point
- load
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 37
- 238000000034 method Methods 0.000 claims abstract description 100
- 238000004364 calculation method Methods 0.000 claims abstract description 27
- 238000005452 bending Methods 0.000 claims abstract description 26
- 230000008569 process Effects 0.000 claims abstract description 19
- 238000011156 evaluation Methods 0.000 claims abstract description 11
- 238000006073 displacement reaction Methods 0.000 claims description 30
- 239000013598 vector Substances 0.000 claims description 28
- 239000011159 matrix material Substances 0.000 claims description 14
- 238000005381 potential energy Methods 0.000 claims description 14
- 239000000463 material Substances 0.000 claims description 10
- 230000008859 change Effects 0.000 claims description 9
- 230000014509 gene expression Effects 0.000 claims description 8
- 230000005540 biological transmission Effects 0.000 claims description 7
- 238000009826 distribution Methods 0.000 claims description 7
- 238000005520 cutting process Methods 0.000 claims description 6
- 230000007812 deficiency Effects 0.000 claims description 6
- 150000001875 compounds Chemical class 0.000 claims description 4
- 239000013256 coordination polymer Substances 0.000 claims description 4
- 238000001069 Raman spectroscopy Methods 0.000 claims description 3
- 230000009351 contact transmission Effects 0.000 claims description 3
- 230000008878 coupling Effects 0.000 claims description 3
- 238000010168 coupling process Methods 0.000 claims description 3
- 238000005859 coupling reaction Methods 0.000 claims description 3
- 238000009795 derivation Methods 0.000 claims description 3
- 230000005489 elastic deformation Effects 0.000 claims description 3
- 238000000926 separation method Methods 0.000 claims description 3
- 238000005259 measurement Methods 0.000 claims description 2
- 238000013461 design Methods 0.000 abstract description 10
- 238000003754 machining Methods 0.000 abstract description 8
- 230000000704 physical effect Effects 0.000 abstract description 7
- 230000006399 behavior Effects 0.000 description 14
- 238000012545 processing Methods 0.000 description 9
- 238000010586 diagram Methods 0.000 description 8
- 238000003801 milling Methods 0.000 description 7
- 238000004519 manufacturing process Methods 0.000 description 4
- 230000001413 cellular effect Effects 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 244000309464 bull Species 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 238000005094 computer simulation Methods 0.000 description 1
- 238000013016 damping Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 230000013011 mating Effects 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 229910052751 metal Inorganic materials 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000010008 shearing Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
- 239000000758 substrate Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Gears, Cams (AREA)
Abstract
本发明公开了一种基于双曲面壳单元模型的数值载荷齿面接触分析(NLTCA)方法,从几何层面分析螺旋锥齿轮齿面基本的结构特性,进行齿面弯曲行为分析;结合以往齿轮有限元建模的基本理论和方法,提出螺旋锥齿轮的双曲面壳单元的建模方法,既考虑齿轮啮合刚度又考虑齿轮接触柔性,以保证NLTCA的数值结果更真实有效;利用Rayleigh‑Ritz法来确定螺旋锥齿轮啮合刚度和接触柔性,为后续的NLTCA数值结果的求解及相关的齿面接触力学性能分析提供基础;采用有限元法和线性规划法相结合的形式求解NLTCA的目标函数,建立加工参数与齿面物理性能评价项的直接关联规律,提供实际加工参数与齿面接触性能评价之间的直接参数驱动关系。求解过程的计算量减少且能保证设计的柔性化和实用性。
Description
技术领域
本发明属于螺旋锥齿轮加工领域,具体为一种基于双曲面壳单元模型的数值载荷齿面接触分析方法。
背景技术
载荷齿面接触分析(LTCA)作为螺旋锥齿轮制造领域一种主要的齿面传动性能设计方法,目前主要的手段还是依托于有限元仿真软件平台进行求解,该仿真载荷齿面接触分析方法过程复杂,对齿轮有限元单元模型精度要求高,整个求解过程计算量大、耗时长,且一旦计算结果出错可修改性差,设计的柔性化和实用性大大降低。
发明内容
本发明的目的在于提供一种能用来构建形性协同制造中加工参数与物理性能的直接函数关系,探究二者之间的关联规律,为MOO加工参数反调及形性协同***建模提供基础的数值载齿面荷接触分析(NLTCA)方法。
本发明提供的这种基于双曲面壳单元模型的数值载荷齿面接触分析方法,包括以下步骤:
i)齿面弯曲特性分析
从几何层面分析其基本的结构特性;
ii)精确的有限元建模
ii.1)建立双曲面壳单元模型;
ii.2)建立双曲面壳单元方程;
ii.3)确定双曲面壳单元模型的基本几何定义和边界约束条件;
iii)齿面啮合刚度和接触柔性的确定
利用基于壳理论的Rayleigh-Ritz法来确定螺旋锥齿轮啮合刚度和接触柔性;
iiii)数值载荷齿面接触分析的精确计算
iiii.1)确定齿面接触点的精确匹配策略;
iiii.2)确定变形协调条件;
iiii.3)确定受力平衡条件;
iiii.4)确定求解数值载荷齿面接触分析问题的线性规划方法;
iiii.5)确定数值载荷齿面接触分析评价项。
上述技术方案的一种实施方式中,针对步骤i),为了建立以精确的有限元模型,首先需要从几何层面分析其基本的结构特性,分别从齿高hG、齿厚tG及齿宽方向wG对齿面进行三维尺寸定义;建立螺旋锥齿轮关于节锥母线上任一点P*的齿线方程:
式中,τG表示齿线的弧长,Sr为径向刀位,rc为刀盘半径,R1为外端锥距,R2为内端锥距。
上述技术方案的一种实施方式中,针对步骤ii.1),将壳单元定位为拥有不同厚度和高度的厚圆锥壳的一部分,其齿宽方向、齿高方向和齿厚方向近似于齿轮的三维尺寸方向;由于较小的半径与厚度及长度与厚度之比,螺旋锥齿轮的齿面弯曲行为可以利用剪切变形理论进行建模;建立局部主***(α,β,z),α和β分别表示中面的正交曲线主坐标,z表示法向厚度根据壳单元中面的不同几何形状,完成双曲面壳单元模型建模,双曲面壳单元模型同时兼容了圆柱壳与圆锥壳模型的特点。
上述技术方案的一种实施方式中,针对步骤ii.2),具体如下:
在坐标系(α,β,z)中,U、V和W分别表示壳单元对应的三个主方向的变形位移量;在全局坐标系O(i1,i2,i3)中,3D壳模型的一个任意点Ps可通过其位移矢量R(α,β,z)定义为
R(α,β,z):=r(α,β)+ziu(α,β)∈R3 (2)
式中,r(α,β)表示中面投影点Ps*的位移矢量,iu(α,β)为向外的单位法向量,Z表示为点Ps到投影点Ps*的距离;由于目前的壳单元大小是有限的,壳单元被限制在以下的尺寸范围内
式中,hG用来表示双曲面壳的全齿厚;
基于不同的几何表达,参考齿面的位移矢量可通过其沿全局坐标系的分量表示
r(α,β)=r1(α,β)i1+r2(α,β)i2+r3(α,β)i3 (4)
式中,i1,i2,i3分别是全局坐标系的单元矢量;此处,控制方程所需的拉曼参数通过参考曲面r(α,β)的第一类基本形式可以表示为
曲面的度量被定义为含有物理含义的数学项;第一类基本系数E、F和G,第二类基本系数L、M和N,通过在局部正交主坐标系(α,β,z)中曲率的参数线α和β的表达,目前存在关系F=M=0,则公式可简化为
式中,Rα(α,β),Rβ(α,β)为齿面的主曲率半径;表明了双曲面壳单元模型的意义,并说明了该模型可以足够精确表达齿面弯曲特性。
在双曲面壳单元的局部局部坐标系中,法向应变和剪切应变关系表示为
在该双曲面壳单元方程中,具有恒定曲率半径的参数化系数可写为
上述技术方案的一种实施方式中,针对步骤ii.3),具体如下:
主曲率半径Rβ表示刀盘切削轨迹上的一段圆弧
Rβ:=Rβ(α,β)=rc (10)
齿高hG和齿厚tG分别假定线性变化且通过线性拟合可以表示其完整形状,则运用关于大端h1和小端h2的一个线性拟合来表示齿高
描述了沿齿宽方向的高度变化,其中τG可通过公式(1)确定。
曲率半径Rα表示从齿根圆Rα1到齿根圆的Rα2的不断变化的曲率半径,根据几何关系,任一点的锥角是恒定的,可以表示为
式中,X:=G和X:=P分别表示大轮和小轮,Σ表示两个齿轮轴的夹角,Rα的线性拟合为
从Rα1到Rα2逐渐减小,和齿厚拟合所呈现的线性关系一致。
在大端齿顶和齿根的齿厚参数t1和t2,在小端齿顶和齿根的齿厚参数t3和t4,可用来进行齿厚的线性拟合,有
或者,也可以表示为
这是由于在齿顶和齿根处具有同样的弯曲行为变化趋势所致。
对于厚度变化的双曲面壳单元模型,双曲面壳底部沿边缘固定而齿顶面、凸面、凹面、大端面和小断面都自由,可以表现螺旋锥齿轮齿面的弯曲行为特性;对于剪切旋转βx和βθG,函数Πm(x)被选用来限制自由。假设将调和函数导入上述方程,使得系列通解的每一项表达成两种方式,一种是关于轴位移x的函数,
式中,多项式函数Πm(x)和Γn(θG)需要通过满足下列条件达到收敛:i)所有必要的边界条件与W,βx和βθG相关;ii)连续完整且独立线性。
上述技术方案的一种实施方式中,针对步骤iii),关于基于壳理论的Rayleigh-Ritz法,基于Bhimaraddi壳的高阶剪切理论的位移假设被用来反映双曲面壳单元模型的弯曲行为特性;位移假设、应变表达包括关于空间几何坐标的齿厚的推导,其必须在应变能中表达出来;横向偏移量W和剪切旋转量βx及βθG,则可利用Rayleigh-Ritz法将势能的一阶变化变为0来计算;此时,忽略了平面应变σz的弹性应变能QSE,σx,、σθG和τxθG与均匀各向同质材料属性相关,其中E为杨氏模量,υ为泊松比,存在势能QPE为弹性应变能QSE和虚功QWF之差。假设QPE=QSE-QWF,则假定势能的一阶差分为0,有
可采用成含未知系数的代数多项式的有限线性组合形式来表示横向偏移量W和剪切旋转量βx及βθG,必须满足预设的边界条件;
针对步骤iii),关于确定螺旋锥齿轮啮合刚度,具体如下:
将边界约束条件代入位移假设(16)-(19)中,通过评法向和剪切应力来计算应变能。而势能QPE和虚功QWF则需要利用高斯-正交法来进行计算,获得线性方程组
式中,子矩阵Kmn(m,n∈[1,3])表示齿轮啮合刚度单元,可以根据材料属性相关的体积积分和预设的调和函数进行确定。而在螺旋锥齿轮接触传动过程中,整个啮合刚度KTOTAL除了由(21)中所确定的大小轮的齿面啮合刚度KP和KG外,还包括Hertz接触刚度KH,可得
KTOTAL=KG+KP+KH (22)
而KH根据Hetrz接触理论进行计算
式中κ1和κ2分别表示对应小轮和大轮的材料参数,D1和D2分别表示在啮合接触点位置两个齿面的曲率半径,而μ为泊松比,E为弹性模量;
在螺旋锥齿轮的载荷接触过程中,单对齿轮的啮合刚度确定如下:在接触坐标系OG(X,Y,Z),K1和K2分别表示有载荷也没有载荷的接触点,其位移矢量和接触力分别为μ1、μ2和F1、F2,则可得到下列关系
当多齿同时发生接触时,由于耦合作用的复杂啮合刚度可表示为
式中,p表示接触过程中发生接触的齿轮对的数目。
针对步骤iii),关于螺旋锥齿轮接触柔性的确定,具体如下:
在公式(59)中,Fm(m∈[1,3])表示外力矢量,来自于无论是点接触还是线接触状态下的载荷;求解公式(59)中的系数矩阵χ中的Amn,Bmn和Cmn,即这些有接触点变形的相互影响系数,就可以确定齿轮接触柔性;假定整个接触线分成Nc部分,每一部分可计算其分配的载荷,则大轮或小轮的解除柔性矩阵可写为
式中,Cis,js(is,js∈[1,Nc])为在某部分is的变形量,由于在某部分js施加载荷所致。
上述技术方案的一种实施方式中,针对步骤iiii.1),在求解齿面接触点过程中,为了保证足够的齿面点的匹配精度,分别考虑各自的齿形误差εi(i=1,2);首先,大小齿轮齿面点M2和M1各自先接触切平面投影,则在可行域内可以按照匹配规则进行目标点的匹配;此时,选择大轮齿面某离散点M2 (ij)(i∈[1,n],j∈[1,m])作为主动点,来匹配另一个小轮齿面的目标点;在匹配过程中,对于任一点M2 (ij),先根据精度要求设置一定的匹配可行域范围半径(FRR)r*,然后来搜寻可行域内另一齿面上的投影点,其目标投影点与主动投影点的距离作为匹配半径,可表示为
r=|rf2(φ2,θ2)-ε2nf2(φ2,θ2)-rf1(φ1,θ1)-ε1nf1(φ1,θ1)| (28)
式中,投影平面即齿轮接触切平面坐标系为Of(if,jf,kf),齿面点的参数表示其点矢rfi和法矢nfi(i=1表示小轮,i=2表示大轮)。
齿面参数化离散点的精确匹配策略为:
1)对于任一点的匹配半径MR即rk>r*,如果不匹配,则调整FRR再次搜寻目标点直至满足条件;
2)对于在可行域内的任意两个候选点,比较其MR,当出现rq<rp≤r*情况时,选取点M1 (q)为匹配目标点;最终,点M2 (ij)和点M1 (q)匹配成功,成为合格的匹配接触对M2 (ij)-M1 (q);式中的参数d表示实际接触过程中目标匹配点的距离。
上述技术方案的一种实施方式中,针对步骤iiii.2),通过齿面接触点精确匹配策略得到的最终齿面接触对,可标记为M1 K-M2 K,而它们之间的初始距离假定为DK,它们之间的弹性变形量ωK(1)和ωK(2);它们在载荷变形情况下的相对距离为Θ;则在接触点位置,有如下的变形协调关系存在
ωK(1)+ωK(1)+DK-Θ≥0(K=1,2,…,N) (29)
式中,“=”表示已经接触;而“>”表示尚未接触,其中ωK(1)和ωK(2)和前文所求的接触单元柔性CKj有关,可得
式中FK表示目标匹配接触对M1 K-M2 K的接触。
上述技术方案的一种实施方式中,针对步骤iiii.3),作用于接触区所有接触点的力之和FK(K=1,2,…,N)需要等于所有外部载荷,则存在以下受力平衡关系
{e}T{F}=P (32)
方程(32)和(29)是用来判定齿面是否接触的约束条件。两个齿面的弹性接触问题可以看成如何实现接触力FK(K=1,2,…,N)在一致接触柔性、外部载荷和变性相对距离的情况下,同时满足变形协调条件和受力平衡条件。
上述技术方案的一种实施方式中,针对步骤iiii.4),具体如下:
由于方程(29)是一个大于或等于0的不等式约束,引入松弛变量Y,则可以变成等式约束,为
[C]{F}+{D}-Θ{e}-[I]{Y}={0} (33)
联立方程(32)和(33),共得到N+1个标量方程;要求解未知参数(F,D,Y),FK≥0,Θ≥0,YK≥0,并要求满足接触条件(29),即存在如下条件
YK=0,FK≥0;YK=0,FK≥0 (34)
利用改进单纯形法来求解满足接触条件(29)的方程(32)和(33),则该齿面弹性接触问题的目标函数表示为一般线性规划:
式中,对于N+1个标量方程需要人为引入一系列非负的变量参数Z=[Z1,Z2,…,ZN,ZN+1],使得等式约束两边的常数项都为非负值;
考虑接触条件(29),采用改进单纯形法,假定目标函数(34)中Z0=-ZN+2,有
将上式作为一个约束条件,ZN+2作为固定的基变量,用分离系数法将目标函数(34)写为单纯形表,得到变换后的标准形式。
上述技术方案的一种实施方式中,针对步骤iiii.5),具体如下:
首先,通过基于双曲面壳模型和Rayleigh-Ritz法确定齿面接触柔性系数CP,G;
然后,利用改进单纯形法的线性规划方法,求解齿面载荷分配和载荷传动误差及其最大载荷传动误差LTEMAX;
其次,根据载荷分配设置载荷边界条件,求解最大齿面接触应力CPMAX;
最后,根据啮合周期内的接触应力发生的时间历程计算出齿面重合度εγ。
本发明为了建立精确的有限元模型,从几何层面分析螺旋锥齿轮齿面基本的结构特性,进行齿面弯曲行为分析;根据齿面弯曲行为分析,结合以往齿轮有限元建模的基本理论和方法,提出了螺旋锥齿轮的双曲面壳单元的建模方法,既考虑齿轮啮合刚度又考虑齿轮接触柔性,以保证NLTCA的数值结果更加真实有效;利用Rayleigh-Ritz法来确定螺旋锥齿轮啮合刚度和接触柔性,同时为后续的NLTCA数值结果的求解及相关的齿面接触力学性能分析提供基础;进行NLTCA的数值精确求解时包括几个关键步骤:确定精确的齿面接触对;建立考虑载荷变形情况下的变形协调关系;建立考虑接触约束条件的受力平衡方程;基于改进单纯形法的NLTCA的线性规划。
总之,本发明充分考虑螺旋锥齿轮的齿面弯曲特性和高使役性能要求,提出了基于精确的双曲面壳有限元模型的求解分析。采用有限元法和线性规划法相结合的形式求解NLTCA的目标函数,建立加工参数与齿面物理性能评价项的直接关联规律,使NLTCA的计算方法可以有效的确定加工参数与齿面物理性能评价项之间的精确参数驱动关系,可以有效弥补现有技术的缺陷,同时提供实际加工参数与齿面接触性能评价之间的直接参数驱动关系,又为考虑齿面几何与物理性能协同优化的设计及今后的协同制造乃至智能制造创造有利条件。
附图说明
图1为本发明螺旋锥齿轮的齿面几何形状示意图。
图2为图1所示齿面的两种主要加工方式示意图。
图3为图1所示齿面的齿线弧长计算示意图。
图4为图1所示齿面的有限元双曲面壳单元模型示意图。
图5为双曲面壳单元的参考坐标系示意图。
图6为双曲面可单元模型的基本几何定义及边界约束条件示意图。
图7为在全局坐标系中单对齿的接触啮合刚度示意图。
图8为齿面匹配点的参数定义示意图。
图9为参数化齿面点的精确匹配策略示意图。
图10为目标匹配点的排序示意图。
图11为确定NLTCA评价项的基本流程示意图。
图12为算例进行15×20齿面离散化后的5×5细分图。
图13为单齿啮合刚度结果比较图。
图14为FEA多齿啮合刚度图。
图15Rayleigh-Ritz法与FEA法的多齿啮合刚度比较结果图。
图16为不同接触单元密度情况下的计算时间图。
图17载荷对齿面啮合刚度的影响图。
具体实施方式
下面结合附图及算例,对上述技术方案进行详细说明。
关于i)齿面弯曲特性分析,具体如下:
螺旋锥齿轮的齿形十分复杂,且具有自身的齿面弯曲行为特性。为了建立精确的有限元模型,首先需要从几何层面分析其基本的结构特性。图1给出了螺旋锥齿轮齿面形状的基本示意图,分别从齿高hG、齿厚tG及齿宽方向wG对齿面进行三维尺寸定义。尤其是对于端面铣削(face-milling)的收缩齿,齿高从小端到大端逐渐增加,齿根厚度从小端到大端的逐渐增加,齿厚在齿宽任意点处沿齿高逐渐减小,这意味着精确的齿轮模型,必须考虑以下弯曲特性:齿宽方向的刚度变化;齿面纵向的曲率变化;齿向螺旋且逐渐收缩变化(格里森齿制)。
众所周知,螺旋锥齿轮从在节锥母线方向看是一条曲线,而其相对于节锥母线的弯曲程度通常用在齿轮节锥曲线上的螺旋角βP来定义。
图2给出了螺旋锥齿轮最常见的两种加工方式,即单分度的端面铣削,主要是适用于Gleason收缩齿,和连续分度的端面滚切,主要适用于Klingelnberg和Oerlikon等高齿。对于端面铣削加工,齿高是变化的,齿面纵向曲率是一段圆弧,其曲率半径等于刀盘的曲率半径,齿轮被刀盘刀具的内外切削刃分别单独加工。对于端面滚切加工齿轮,齿高是恒定的,纵向齿线是一个延伸的外摆线,齿轮的凹面P1面和凸面P2同时被刀具内外切削刃加工。中点P的螺旋角定义为βP。在△OOcP,根据余弦定理,有以下关系
式中R为节锥,rc为刀盘半径,Sr为径向刀位。假设P*为节线任一点且其锥距为RX,其螺旋角βX根据关于△OOcP*的余弦定理可得
联立以上两个方程(1)和(1-1),则节锥线上任意一点的螺旋角表示为
因此,参数rc和βX则可参数定义螺旋锥齿轮的齿线。计算齿线的弧长τG,如图3所示。
在半径为rc的圆Oc上,PAOc=PBOc=rc,有以下关系
τG=∠PAOcPB·rc=(∠PAOcO-∠PBOcO)·rc (1-3)
在△PAOcO和△PBOcO,OcO=Sr,PAO=R2为内端锥距,PAO=R1为外端锥距。根据余弦定理可得:
联立以上两个方程(1-4,1-5)有
这表明齿线是恒定的,且可以被用来确定所提出的有限元壳单元模型的几何尺寸。由此可得,壳单元比一般的梁单元或者板单元更能精确的表达螺旋锥齿轮弯曲行为特性。
关于ii.1)建立精确的双曲面壳单元模型
精确的有限元模型是齿轮接触力学性能求解的关键步骤。有限元法方法就是把复杂的整体结构离散到有限个单元(Finite Element),再把这种理想化的假定和力学控制方程施加于结构内部的每一个单元,然后通过单元分析组装得到结构总刚度方程,再通过边界条件和其他约束解得结构总反应。
因此,数值载荷齿面接触分析(NLTCA)作为分析齿轮接触力学性能的重要手段,其齿轮有限元精确模型的建立也十分重要。本发明根据前述螺旋锥齿轮齿面弯曲行为特性分析,结合以往齿轮有限元建模的基本理论和方法,提出了螺旋锥齿轮双曲面壳单元的建模方法,具体如下:
ii.1.1)建立精确的双曲面壳模型
目前,螺旋主齿轮的有限元单元建模研究,经历了一个逐渐由浅入深的过程:从简化的悬臂梁单元到收缩板单元到目前的圆柱壳单元。区别以往有限元模型结构,本发明提出更为精确的双曲面壳单元模型,如图4所示。壳单元可以看成拥有不同厚度和高度的厚圆锥壳的一部分,其齿宽方向、齿高方向和齿厚方向近似于齿轮的三维尺寸方向。由于较小的半径与厚度及长度与厚度之比,螺旋锥齿轮的齿面弯曲行为可以利用剪切变形理论进行建模。特别地,在主要的集中壳单元模型中,双曲面壳单元模型被选用来完成螺旋锥齿轮的有限元建模,因为它能更加精确的反映出齿面的弯曲特性。
在图4所示的模型中,双曲面壳单元模型同时兼容了圆柱壳单元与圆锥壳单元模型的特点。建立局部主***(α,β,z),α和β分别表示中面的正交曲线主坐标,z表示法向厚度。由于端面铣削的齿面具有更复杂的弯曲特性,且厚度可变,本发明以此为例进行论述。而对于端面滚切等高齿面模型,只需要设定法向厚度恒定即可。值得注意的是,为了探讨给定齿面上逐点变化的曲率半径的双曲面几何特性,本发明给出壳单元中面的不同几何形状,以完成快速可靠且简便的有限元结构建模。
ii.2)建立双曲面壳单元方程
前面已经从螺旋锥齿轮几何层面分析其基本的结构特性,它在齿高、齿宽和齿厚方向的三维尺寸都是变化的,这也符合壳单元模型中的双曲面壳单元模型的几何特性。因此,考虑齿面弯曲特性的螺旋锥齿轮双曲面壳模型的几何特点,图5给出了其基本参数化描述。在坐标系(α,β,z)中,U、V和W分别表示了壳单元对应的三个主方向的变形位移量。在全局坐标系O(i1,i2,i3)中,3D壳模型的一个任意点Ps可通过其位移矢量R(α,β,z)定义为
R(α,β,z):=r(α,β)+ziu(α,β)∈R3 (2)
式中,r(α,β)表示壳单元中面投影点Ps*的位矢,iu(α,β)为向外的单位法向量,Z表示为点Ps到投影点Ps*的距离。由于目前的壳单元大小是有限的,被限制在以下的尺寸范围内
式中,hG表示双曲面壳的全齿厚。
基于不同的几何表达,参考齿面的位移矢量可通过其沿全局坐标系的分量表示
r(α,β)=r1(α,β)i1+r2(α,β)i2+r3(α,β)i3 (4)
式中,i1,i2,i3分别是全局坐标系的单元矢量。此处,控制方程所需的拉曼参数通过参考曲面r(α,β)的第一类基本形式可以表示为
曲面的度量被定义为含有物理含义的数学项。第一类基本系数E、F和G可构成一个矩阵CP
可以进行倒置,由于其行列式不为0。通过公式(5)和(5-1)可得,
由此,易于计算向外单位法矢为
同理,包含第二类基本系数L、M和N矩阵CCP为
联立(3)和(5),有一个诸如下三阶处理ΘP
包含齿面的曲率,其确定允许定义高斯曲率ΩG和平均曲率ΩM。故有以下关系
通过求解,在α和β方向的主曲率kα(α,β)和kβ(α,β)为
式中,Rα(α,β),Rβ(α,β)为齿面的主曲率半径。通过在局部正交主坐标系(α,β,z)中曲率的参数线α和β的表达,目前存在关系F=M=0,则公式(5-7)可简化为
表明了双曲面壳单元模型的意义,并说明了该模型可以足够精确表达螺旋锥齿面的弯曲特性。
在双曲面壳的局部坐标系中,三维应变分量表示为
另外,便于论述,引入如下表达
ε:=ε(α,β,z,t)=[εαα εββ εzz γαβ γαz γβz]T (6-3)
U=[U(α,β,t)V(α,β,t)W(α,β,t)]T (6-4)
联立以上(6-1)-(6-5),法向应变和剪切应变关系表示为
在该双曲面壳方程中,具有恒定曲率半径的参数化系数可写为
关于ii.3)确定双曲面壳单元模型的基本几何定义和边界约束条件,具体如下:
图6给出了建立的双曲面壳模型的基本几何定义和边界约束条件。其中,β=θG表示任意齿面接触点从刀盘中心到小端面之间的角度,α=x则表示在特定
接触点的齿高。参考图3,主曲率半径Rβ表示刀盘切削轨迹上的一段圆弧,
Rβ:=Rβ(α,β)=rc (10)
在该模型中,齿高hG和齿厚tG分别假定线性变化且通过线性拟合可以表示其完整形状。则运用关于大端h1和小端h2的一个线性拟合来表示齿高为,
描述了沿齿宽方向的高度变化,其中τG可通过公式(1)确定。
而另一个表示曲率半径Rα,表示从齿根圆Rα1到齿根圆的Rα2的不断变化的曲率半径。根据几何关系,任一点的锥角是恒定的,可以表示为
式中,X:=G和X:=P分别表示大轮和小轮,Σ表示两个齿轮轴的夹角。由此,Rα的线性拟合表示为
从Rα1到Rα2逐渐减小,和齿厚拟合所呈现的线性关系一致。
根据图6的几何定义,在大端齿顶和齿根的齿厚参数t1和t2,在小端齿顶和齿根的齿厚参数t3和t4,可用来进行齿厚的线性拟合,有
或者,也可以表示为
这是由于在齿顶和齿根处具有同样的弯曲行为变化趋势所致。
对于厚度变化的双曲面壳单元模型,图6设定了相应的边界条件。这样,双曲面壳底部沿边缘固定而齿顶面、凸面、凹面、大端面和小断面都自由,可以表现螺旋锥齿轮齿面的弯曲行为特性。剪切旋转βx和βθG,函数Πm(x)被选用来限制自由,该约束条件为
式中,H表示齿高。对于横向偏移量(W),函数Πm(x)被选用来固定自由,该约束为
同理,函数Γm(θG)被选用来释放自由,该约束条件为
假设将调和函数导入上述方程,使得系列通解的每一项表达成两种方式,一种是关于轴位移x的函数,
式中,多项式函数Πm(x)和Γn(θG)需要通过满足下列条件达到收敛:1)所有必要的边界条件与W,βx和βθG相关;2)连续完整且独立线性。
iii)关于齿面啮合刚度和接触柔性的确定
基于双曲面壳单元模型的螺旋锥齿轮及设定的边界条件,本文既考虑齿轮啮合刚度也考虑齿轮解除柔性,以保证提出的NLTCA的数值结果更加真实有效。根据最小势能原理,如果能够列出所有的几何可能位移,那么使总势能取最小值的那一组位移就是真实位移。为此,本发明利用Rayleigh-Ritz法来确定螺旋锥齿轮啮合刚度和接触柔性,同时为后续的NLTCA数值结果的求解及相关的齿面接触力学性能分析提供基础。具体步骤如下:
iii.1)基于壳理论的Rayleigh-Ritz法
在本文中,基于Bhimaraddi壳的高阶剪切理论的位移假设被用来反映壳模型的弯曲行为特性。考虑从齿顶面到齿底面的齿厚位置横向剪切应力的抛物线变化,详细的表达为
基于Bhimaraddi壳理论,位移假设表示为
应变表达包括关于空间几何坐标的齿厚的推导,其必须在应变能中表达出来。横向偏移量(W)和剪切旋转量βx及βθG,则可利用Rayleigh-Ritz法将势能的一阶变化变为0来计算。此时,忽略了平面应变σz的弹性应变能QSE为
σx,、σθG和τxθG与均匀各向同质材料属性相关,可表示为
式中E为杨氏模量,υ为泊松比。根据公式(19-3)-(19-12),存在
由于导致可变性的外力Fp所做的功QWF为
势能QPE为弹性应变能QSE和虚功QWF之差。假设QPE=QSE-QWF,则假定势能的一阶差分为0,有
可采用成含未知系数的代数多项式的有限线性组合形式来表示横向偏移量(W)和剪切旋转量βx及βθG,必须满足预设的边界条件。
iii.2)关于啮合刚度的确定
将边界约束条件(16)-(19)代入位移假设(19-6)-(19-8)中,通过评法向和剪切应力来计算应变能(56)。而势能QPE和虚功QWF则需要利用高斯-正交法来进行计算,获得线性方程组
式中,子矩阵Kmn(m,n∈[1,3])表示齿轮啮合刚度单元,可以根据材料属性相关的体积积分和预设的调和函数进行确定。而在螺旋锥齿轮接触传动过程中,整个啮合刚度KTOTAL除了由(59)中所确定的大小轮的齿面啮合刚度KP和KG外,还包括Hertz接触刚度KH,可得
KTOTAL=KG+KP+KH (22)
而KH根据Hetrz接触理论进行计算
式中κ1和κ2分别表示对应小轮和大轮的材料参数,D1和D2分别表示在啮合接触点位置两个齿面的曲率半径,而μ为泊松比,E为弹性模量。
在螺旋锥齿轮的载荷接触过程中,单对齿的啮合刚度如图7所示。在接触坐标系OG(X,Y,Z),K1和K2分别表示有载荷也没有载荷的接触点,其位移矢量和接触力分别为μ1、μ2和F1、F2,则可得到下列关系
当多齿同时发生接触时,由于耦合作用的复杂啮合刚度可表示为
式中,p表示接触过程中发生接触的齿轮对的数目。
iii.2)关于齿轮接触柔性的确定
在公式(21)中,Fm(m∈[1,3])表示外力矢量,来自于无论是点接触还是线接触状态下的载荷。求解公式(21)中的系数矩阵χ中的Amn,Bmn和Cmn,即这些有接触点变形的相互影响系数,就可以确定齿轮接触柔性。
由于每一对接触点位置都要进行计算,关于这些变形函数的未知系数的求解是十分耗时的。假定整个接触线分成Nc部分,每一部分可计算其分配的载荷,则大轮或小轮的接触柔性矩阵可写为
式中,Cis,js(is,js∈[1,Nc])为在某部分is的变形量,由于在某部分js施加载荷所致。
关于iiii)NLTCA的精确计算
基于精确有限元建模和啮合刚度与接触柔性矩阵的确定,就可以进行NLTCA的数值精确求解。其过程主要包括以下几个关键步骤:1)确定精确的齿面接触对;2)建立考虑载荷变形情况下的变形协调关系;3)建立考虑接触约束条件的受力平衡方程;4)基于改进单纯形法的NLTCA的线性规划。改进单纯形法是1953年美国数学家G.B.丹齐克提出了改进单纯形法,其基本步骤和单纯形法大致相同,主要区别是在逐次迭代中不再以高斯消去法为基础,而是由旧基阵的逆去直接计算新基阵的逆,再由此确定检验数。这样做可以减少迭代中的累积误差,提高计算精度,同时也减少了在计算机上的存储量。
步骤iiii)包括以下分步:
iiii.1)齿面接触点的精确匹配策略
本申请人已经就含误差齿面接触分析(eTCA)的求解申请了专利,其中公开了相应的TCA流程和接触点确定方法。本发明只重点论述NLTCA中的齿面接触点的精确匹配问题。在求解齿面接触点过程中,为了保证足够的齿面点的匹配精度,考虑了分别考虑各自的齿形误差εi(i=1,2)。首先,大小齿轮齿面点M2和M1各自先接触切平面投影,则在可行域内可以按照一定的匹配规则进行目标点的匹配。此时,选择大轮齿面某离散点M2 (ij)(i∈[1,n],j∈[1,m])作为主动点,来匹配另一个小轮齿面的目标点。在匹配过程中,对于任一点M2 (ij),先根据精度要求设置一定的匹配可行域范围半径(FRR)r*,然后来搜寻可行域内另一齿面上的投影点。如图8所示,其目标投影点与主动投影点的距离作为匹配半径,可表示为
r=|rf2(φ2,θ2)-ε2nf2(φ2,θ2)-rf1(φ1,θ1)-ε1nf1(φ1,θ1)| (28)
式中,投影平面即齿轮接触切平面坐标系为Of(if,jf,kf),齿面点的参数表示其点矢rfi和法矢nfi(i=1表示小轮,i=2表示大轮)。
因此,如图9所示,齿面参数化离散点的精确匹配策略为:
i)对于任一点的MR即rk>r*,如果不匹配,则调整FRR再次搜寻目标点直至满足条件;
ii)对于在可行域内的任意两个候选点,比较其MR,当出现rq<rp≤r*情况时,选取点M1 (q)为匹配目标点。最终,点M2 (ij)和点M1 (q)匹配成功,成为合格的匹配接触对M2 (ij)-M1 (q)。
很明显,对于两个候选点,存在完全重合关系nf1 (ij)=nf2 (ij)即r*=0这种情况是最优最有效解,但是这只是特例。另外,如果当前匹配精度不符合要求或者需要更高精度的匹配,本文提供了更为优化的匹配策略。当在可行域内不存在一个匹配点时,FRR可以增加进行再次匹配。只需要在上次未成功匹配点附近区域以较小的MR进行搜索。但是,为了提高匹配精度,搜索半径可以逐次增大。因此,可以再上一次目标区域进行细分,搜寻满足条件的目标点。基本的细分规则为,以上次候选点为中心进行一个简单的5×5网格等分;如果失败,则再次进行(4n+1)×(4n+1)等分的细分(n为细分次数),直至满足条件。
最后,所有的目标匹配点被确定,可进行精确编号排序,如图10所示。首先选取大小轮齿面点作为输入,并从点Mi (1-1)到Mi (n-m)进行精确编号,尤其是某一行的末尾点和下一行的起始点要注意衔接,诸如Mi (1-m)和Mi (2-1);然后,进行有效的匹配;每一对目标匹配点按照原来的序列号进行保存。通过匹配和排序后,下列的优化关系可用来求解对应的齿面接触点
mind=|rf1(φ1,θ1)-rf2(φ2,θ2)| (28-1)
式中的参数d表示实际接触过程中目标匹配点的距离。
iiii.2)变形协调条件
通过齿面接触点的精确匹配策略得到最终的齿面接触对,可标记为M1 K-M2 K,而它们之间的初始距离假定为DK,它们之间的弹性变形量ωK(1)和ωK(2);它们在载荷变形情况下的相对距离为Θ。则在接触点位置,有如下的变形协调关系存在
ωK(1)+ωK(1)+DK-Θ≥0(K=1,2,…,N) (29)
式中,“=”表示已经接触;而“>”表示尚未接触,其中ωK(1)和ωK(2)和前文所求的接触单元柔性CKj有关,可得
式中FK表示目标匹配接触对M1 K-M2 K的接触力。联立(68)-(69)可以得到:
[C]{F}+{D}-Θ{e}≥{0} (30-2)
式中,有
[C]=[CKj]=[CK(1)j+CK(2)j],{F}={F1,F2,…,FK,…,FN}T
{D}={D1,D2,…,DK,…,DN}T,{e}={1,1,…,1,…,1}T,{0}={0,0,…,0,…,0}T
iiii.3)受力平衡条件
作用于接触区所有接触点的力之和FK(K=1,2,…,N)需要等于所有外部载荷,则存在以下受力平衡关系
{e}T{F}=P (32)
方程(32)和(30-2)是用来判定齿面是否接触的约束条件。两个齿面的弹性接触问题可以看成如何实现接触力FK(K=1,2,…,N)在一致接触柔性、外部载荷和变性相对距离的情况下,同时满足变形协调条件和受力平衡条件。
iiii.4)求解NLTCA问题的线性规划方法
方程(30-2)是一个大于或等于0的不等式约束,引入松弛变量Y,则可变成等式约束为
[C]{F}+{D}-Θ{e}-[I]{Y}={0} (33)
联立方程(32)和(33),共得到N+1个标量方程。要求解未知参数(F,D,Y),FK≥0,Θ≥0,YK≥0,并要求满足接触条件(29),即存在如下条件
YK=0,FK≥0;YK=0,FK≥0 (34)
本发明利用改进单纯形法来求解满足接触条件(29)的方程(32)和(33),则该齿面弹性接触问题的目标函数表示为一般线性规划:
式中,对于N+1个标量方程需要人为引入一系列非负的变量参数Z=[Z1,Z2,…,ZN,ZN+1],使得等式约束两边的常数项都为非负值。
在本发明所提出的螺旋锥齿轮NLTCA求解过程中,与一般的线性规划不同,还有考虑接触条件(29),采用改进的单纯形法,假定目标函数(34)中Z0=-ZN+2,有
将上式作为一个约束条件,ZN+2作为固定的基变量,用分离系数法将目标函数(34)写为单纯形表形式,采用改进单纯形法进行求解。
iiii.5)NLTCA评价项的确定
图11给出了确定NLTCA评价项的基本流程,其关键点为:1)有限元建模;2)啮合刚度、解除柔性求解;3)NLTCA方程设置及其求解。
首先,通过基于双曲面壳模型和Rayleigh-Ritz法可以确定齿面接触柔性系数CP,G;
然后利用改进单纯形法的线性规划方法,就可以求解齿面载荷分配和载荷传动误差及其最大载荷传动误差LTEMAX;
其次,根据载荷分配设置载荷边界条件,就可以求解最大齿面接触应力CPMAX;
最后,根据啮合周期内的接触应力发生的时间历程就可以计算出齿面重合度εγ;
很显然,该NLTCA的计算方法可以有效的确定加工参数与齿面物理性能评价项之间的精确参数驱动关系。
算例
关于齿面接触点匹配,具体如下:
表3给出了齿坯设计参数。表4给出了属于端面铣削加工(face-milling)范畴的弧齿锥齿轮SGM加工调整卡参数及本发明所需的加工参数。根据以上设计参数就可以完成精确的螺旋锥建模并建立齿面TCA方程。根据本申请人已经就含误差齿面接触分析(eTCA)的求解申请了专利,其中公开的相应的TCA流程和接触点确定方法结合本发明给出的齿面接触点匹配策略就可以完成齿面接触点的精确求解。表5给出了15×20离散化齿面初始接触点的精确匹配结果。而表6则给出了更精确的15×20离散化齿面经过5×5细分的初始接触点的匹配结果。值得注意的是,此处的ease-off采用的是它的原始定义,表示大小齿轮面与它们共轭接触位置的偏差。
表3齿坯设计基本参数
表4弧齿锥齿轮SGM加工调整卡参数
表5 15×20离散化齿面初始接触点的精确匹配
表6 15×20离散化齿面经过5×5细分的初始接触点的精确匹配
表7以上两种初始接触点精确匹配的比较
其中,在表7的匹配结果中,目标匹配点(PTs)为大轮M2 (8,10)和小轮M1 (7,11),而FRR为0.25mm,计算的MR为0.173654。在表7更高精度匹配中,图12给出了具体的5×5细分示意图。设定一个新的FRR,以上一次的目标匹配点M1 (7,11)为中心进行5×5细分,搜寻到最后的目标匹配点为M1 (7,11) -1 (4,4)。通过以上两种匹配,比较了其匹配之后求解最小距离,一次来验证其正确性。其中,目标匹配点M2 (8,10)-M1 (7,11)的为0.1683314mm,而M2 (8,10)-M1 (7,11) -1 (4,4)的为0.0267345mm。因此,根据给出的齿面接触点匹配策略,最终会选取接触对M2 (8,10)-M1 (7 ,11) -1 (4,4)作为结果输出来求解精确初始点。
关于双曲面壳单元和Rayleigh-Ritz法,具体如下:
一副端面铣削加工的准双曲面齿轮齿面基本参数如表8和表9所示。其中,准双曲面齿轮的HGT调整卡中,大轮采用展成法加工,小轮采用刀倾法加工,选择机械式的NO.116机床采用较为普遍的五刀法加工。基于提供的设计参数就可以齿轮副建模,用于进行有限元单元设定和双曲面壳方程求解。这里,为了来验证本文提出的基于双曲面壳理论的Rayleigh-Ritz方法的可行性,选取了NLTCA中的重要的齿面接触刚度值作为基本指标,以成熟的商用有限元分析(FEA)软件例如的LTCA计算结果为参考,进行比较评价。
表8齿坯设计参数
表9 HGT加工调整卡参数
在NLTCA设置中,大轮凸面和凹面被设定为接触齿面,基本的材料属性为杨氏模量E=2.09×105Mpa,泊松比υ=0.3,切向接触方向的摩察系数设为0.05,输入载荷为200Nm。采用六面体一次减缩积分单元C3D8R来定义单元属性。图13给出了单齿啮合刚度的比较结果。在分析类型的定义中,采用隐式静态分析算法,计算步长设为0.1,即需要10步完成,以便在迭代过程中缓慢增加载荷并获得收敛。首先,齿轮围绕其自身轴线的自由度(DOF)给出少量的旋转,并且小齿轮上施加固定的边界条件。其次,在初始负载的情况下,取消齿轮旋转并且围绕其轴线增加扭矩来完成接触条件。第三,小齿轮的旋转角度围绕其轴线设置,以便它能够驱动齿轮旋转,并且采用递增方法来获得结果。在边界条件的定义中,除参考点处轴线旋转自由度以外的所有自由度都受到齿轮和小齿轮的所有自由度的限制。
通过齿面建模,考虑齿面弯曲行为特性的变厚度双曲面壳模型建立。建立的有限元壳模型的几何参数及边界条件分别确定。基于壳理论的Rayleigh-Ritz就可以确定齿面的齿面接触柔性。值得强调的是,由于齿面接触点的数目较少,齿面接触柔性和齿面啮合刚度是一对相反数,可直接将接触柔性求逆矩阵,以获得齿面啮合刚度。而刚度值又可以通过FEA软件仿真直接求取,故可选择这两种啮合刚度确定方法作为比较。单齿啮合刚度的比较结果,如图15所示。其中,FEA结果中,最大值为800000N/mm,最小值为180000N/mm。而本章提出的Rayleigh-Ritz相对于FEA法的残差值中,RMSE为3.02%,最大值为5.45%,最小值为-4.43%。而值得注意的是,用于FEA计算需要31.45个小时的CPU时间,基本配置为2.5GHz处理器和1Gbyte RAM,而Rayleigh-Ritz法只需要19.36分钟。
多齿啮合刚度可以根据单齿啮合刚度叠加而得,图14给出了FEA法中确定的多齿啮合刚度。其中,单齿接触主要发生在啮合周期的第I和III部分,即蓝线所表示的曲线;而第II部分则为多齿接触,其中最大值达到了1.1135×106N/mm,最小值为6.869×105N/mm。
图15给出了本发明提出的Rayleigh-Ritz法和FEA方法的多齿啮合刚度的比较结果。有11个接触网格节点的数值能分别能提取出来拟合成为多齿啮合曲线。有图可知,两种方法具有相同的变化趋势。比较结果显示:在第I部分,残差的RMSE为3.1672%,最大值为5.71%,最小值为-2.43%;在第II部分,残差的RMSE为2.3746%,最大值为4.01%,最小值为-3.61%;在第III部分,残差的RMSE为3.5083%,最大值为6.38%,最小值为-3.57%。但是,通过比较发现,在同样表示单齿啮合刚度的第I和III部分,尽管数值结果基本近似,但依然有所出入,可能是NLTCA计算过程中接触区域的单元网格分布密度稍微不同所致。
在基于有限元双曲面壳模型的NLTCA计算过程中,接触单元密度会直接影响其计算结果。表10给出了不同接触单元密度的计算结果比较。当接触网格单元为102时,最大接触啮合刚度为108629N/mm,参考图17中的最大值,计算误差为0.00732%,然而计算速度显著下降。图16给出了不同接触单元密度情况下的计算时间。结果表明,考虑到计算精度和效率,合理设置有限接触单元数目很有必要。图17表示了接触载荷对齿面啮合刚度的影响。随着载荷的增加,啮合周期变长且啮合刚度变大,则也意味着齿面接触柔性变小。毋庸置疑,在齿轮接触过程中,载荷对齿面接触性能会产生重大的影响。当然,在接触过程中,还有很多重要因素,例如齿面摩擦和阻尼等,由于其复杂性不再赘述。
表10不同接触单元密度的计算结果比较
Claims (7)
1.一种基于双曲面壳单元模型的数值载荷齿面接触分析方法,包括以下步骤:
i)齿面弯曲特性分析
从几何层面分析其基本的结构特性;
ii)精确的有限元建模
ii.1)建立双曲面壳单元模型;
ii.2)建立双曲面壳单元方程;
ii.3)确定双曲面壳单元模型的基本几何定义和边界约束条件;
iii)齿面啮合刚度和接触柔性的确定
利用基于壳理论的Rayleigh-Ritz法来确定螺旋锥齿轮啮合刚度和接触柔性;
iiii)数值载荷齿面接触分析的精确计算
iiii.1)确定齿面接触点的精确匹配策略;
iiii.2)确定变形协调条件
通过齿面接触点精确匹配策略得到的最终齿面接触对,可标记为M1 K-M2 K,而它们之间的初始距离假定为DK,它们之间的弹性变形量ωK(1)和ωK(2);它们在载荷变形情况下的相对距离为Θ;则在接触点位置,有如下的变形协调关系存在
ωK(1)+ωK(1)+DK-Θ≥0 (K=1,2,…,N) (29)
式中,“=”表示已经接触;而“>”表示尚未接触,其中ωK(1)和ωK(2)和前文所求的接触单元柔性CKj有关,可得
式中FK表示目标匹配接触对M1 K-M2 K的接触;
iiii.3)确定受力平衡条件
作用于接触区所有接触点的力之和FK(K=1,2,…,N)需要等于所有外部载荷,则存在以下受力平衡关系
{e}T{F}=P (32)
方程(32)和(29)是用来判定齿面是否接触的约束条件; 两个齿面的弹性接触问题可以看成如何实现接触力FK(K=1,2,…,N)在一致接触柔性、外部载荷和变性相对距离的情况下,同时满足变形协调条件和受力平衡条件;
iiii.4)确定求解数值载荷齿面接触分析问题的线性规划方法;
由于方程(29)是一个大于或等于0的不等式约束,引入松弛变量Y,则可以变成等式约束,为
[C]{F}+{D}-Θ{e}-[I]{Y}={0} (33)
式中,
{Y}={Y1,Y2,…,YK,…,YN}T
联立方程(32)和(33),共得到N+1个标量方程;要求解未知参数(F,D,Y),FK≥0,Θ≥0,YK≥0,并要求满足接触条件(29),即存在如下条件
YK=0,FK≥0;YK=0,FK≥0 (34)
利用改进单纯形法来求解满足接触条件(29)的方程(32)和(33),则该齿面弹性接触问题的目标函数表示为一般线性规划:
式中,对于N+1个标量方程需要人为引入一系列非负的变量参数Z=[Z1,Z2,…,ZN,ZN+1],使得等式约束两边的常数项都为非负值;
考虑接触条件(29),采用改进单纯形法,假定目标函数(34)中Z0=-ZN+2,有
将上式作为一个约束条件,ZN+2作为固定的基变量,用分离系数法将目标函数(34)写为单纯形表;在横线下面的每一行代表一个约束条件,共有N+2个,竖线右边为原等约束条件(34)的右边常数项;将第N+2行减去前N+1行的和,并表示在第N+2行中,得到变换后的标准形式;
iiii.5)确定数值载荷齿面接触分析评价项:
首先,通过基于双曲面壳模型和Rayleigh-Ritz法确定齿面接触柔性系数CP,G;
然后,利用改进单纯形法的线性规划方法,求解齿面载荷分配和载荷传动误差及其最大载荷传动误差LTEMAX;
其次,根据载荷分配设置载荷边界条件,求解最大齿面接触应力CPMAX;
最后,根据啮合周期内的接触应力发生的时间历程计算出齿面重合度εγ。
3.如权利要求2所述的基于双曲面壳单元模型的数值载荷齿面接触分析方法,其特征在于:针对步骤ii.1),将壳单元定位为拥有不同厚度和高度的厚圆锥壳的一部分,其齿宽方向、齿高方向和齿厚方向近似于齿轮的三维尺寸方向;由于较小的半径与厚度及长度与厚度之比,螺旋锥齿轮的齿面弯曲行为利用剪切变形理论进行建模;建立局部主***(α,β,z),α和β分别表示中面的正交曲线主坐标,z表示法向厚度根据壳单元中面的不同几何形状,完成双曲面壳单元模型建模,双曲面壳单元模型同时兼容了圆柱壳与圆锥壳模型的特点。
4.如权利要求3所述的基于双曲面壳单元模型的数值载荷齿面接触分析方法,其特征在于,针对步骤ii.2),具体如下:
在坐标系(α,β,z)中,U、V和W分别表示壳单元对应的三个主方向的变形位移量;在全局坐标系O(i1,i2,i3)中,3D壳模型的一个任意点Ps可通过其位移矢量R(α,β,z)定义为
R(α,β,z):=r(α,β)+ziu(α,β)∈R3 (2)
式中,r(α,β)表示中面投影点Ps*的位移矢量,iu(α,β)为向外的单位法向量,Z表示为点Ps到投影点Ps*的距离;由于目前的壳单元大小是有限的,壳单元被限制在以下的尺寸范围内
式中,hG用来表示双曲面壳的全齿厚;
基于不同的几何表达,参考齿面的位移矢量可通过其沿全局坐标系的分量表示
r(α,β)=r1(α,β)i1+r2(α,β)i2+r3(α,β)i3 (4)
式中,i1,i2,i3分别是全局坐标系的单元矢量;此处,控制方程所需的拉曼参数通过参考曲面r(α,β)的第一类基本形式可以表示为
曲面的度量被定义为含有物理含义的数学项;第一类基本系数E、F和G,第二类基本系数L、M和N,通过在局部正交主坐标系(α,β,z)中曲率的参数线α和β的表达,目前存在关系F=M=0,则公式可简化为
式中,Rα(α,β),Rβ(α,β)为齿面的主曲率半径;表明了双曲面壳单元模型的意义,并说明了该模型可以足够精确表达齿面弯曲特性;
在双曲面壳单元的局部坐标系中,法向应变和剪切应变关系表示为
在该双曲面壳单元方程中,具有恒定曲率半径的参数化系数可写为
5.如权利要求1所述的基于双曲面壳单元模型的数值载荷齿面接触分析方法,其特征在于,针对步骤ii.3),具体如下:
主曲率半径Rβ表示刀盘切削轨迹上的一段圆弧
Rβ:=Rβ(α,β)=rc (10)
齿高hG和齿厚tG分别假定线性变化且通过线性拟合可以表示其完整形状,则运用关于大端h1和小端h2的一个线性拟合来表示齿高
描述了沿齿宽方向的高度变化,其中τG可通过公式(1)确定;
曲率半径Rα表示从齿根圆Rα1到齿根圆的Rα2的不断变化的曲率半径,根据几何关系,任一点的锥角是恒定的,可以表示为
式中,X:=G和X:=P分别表示大轮和小轮,Σ表示两个齿轮轴的夹角;由此,Rα的线性拟合表示为
从Rα1到Rα2逐渐减小,和齿厚拟合所呈现的线性关系一致;
在大端齿顶和齿根的齿厚参数t1和t2,在小端齿顶和齿根的齿厚参数t3和t4,可用来进行齿厚的线性拟合,有
或者,也可以表示为
这是由于在齿顶和齿根处具有同样的弯曲行为变化趋势所致;
对于厚度变化的双曲面壳单元模型,双曲面壳底部沿边缘固定而齿顶面、凸面、凹面、大端面和小断面都自由,可以表现螺旋锥齿轮齿面的弯曲行为特性;对于剪切旋转βx和βθG,函数Πm(x)被选用来限制自由; 假设将调和函数导入上述方程,使得系列通解的每一项表达成两种方式,一种是关于轴位移x的函数,
式中,多项式函数Πm(x)和Γn(θG)需要通过满足下列条件达到收敛:i)所有必要的边界条件与W,βx和βθG相关;ii)连续完整且独立线性。
6.如权利要求5所述的基于双曲面壳单元模型的数值载荷齿面接触分析方法,其特征在于,针对步骤iii),关于基于壳理论的Rayleigh-Ritz法,基于Bhimaraddi壳的高阶剪切理论的位移假设被用来反映双曲面壳单元模型的弯曲行为特性;位移假设、应变表达包括关于空间几何坐标的齿厚的推导,其必须在应变能中表达出来;横向偏移量W和剪切旋转量βx及βθG,则可利用Rayleigh-Ritz法将势能的一阶变化变为0来计算;此时,忽略了平面应变σz的弹性应变能QSE,σx,、σθG和τxθG与均匀各向同质材料属性相关,其中E为杨氏模量,υ为泊松比,存在势能QPE为弹性应变能QSE和虚功QWF之差; 假设QPE=QSE-QWF,则假定势能的一阶差分为0,有
可采用成含未知系数的代数多项式的有限线性组合形式来表示横向偏移量W和剪切旋转量βx及βθG,必须满足预设的边界条件;
针对步骤iii),关于确定螺旋锥齿轮啮合刚度,具体如下:
将边界约束条件代入位移假设(16)-(19)中,通过评法向和剪切应力来计算应变能;而势能QPE和虚功QWF则需要利用高斯-正交法来进行计算,获得线性方程组
式中,子矩阵Kmn(m,n∈[1,3])表示齿轮啮合刚度单元,可以根据材料属性相关的体积积分和预设的调和函数进行确定; 而在螺旋锥齿轮接触传动过程中,整个啮合刚度KTOTAL除了由(21)中所确定的大小轮的齿面啮合刚度KP和KG外,还包括Hertz接触刚度KH,可得
KTOTAL=KG+KP+KH (22)
而KH根据Hetrz接触理论进行计算
式中κ1和κ2分别表示对应小轮和大轮的材料参数,D1和D2分别表示在啮合接触点位置两个齿面的曲率半径,而μ为泊松比,E为弹性模量;
在螺旋锥齿轮的载荷接触过程中,单对齿轮的啮合刚度确定如下:在接触坐标系OG(X,Y,Z),K1和K2分别表示有载荷也没有载荷的接触点,其位移矢量和接触力分别为μ1、μ2和F1、F2,则可得到下列关系
当多齿同时发生接触时,由于耦合作用的复杂啮合刚度可表示为
式中,p表示接触过程中发生接触的齿轮对的数目;
针对步骤iii),关于螺旋锥齿轮接触柔性的确定,具体如下:
在公式(21)中,Fm(m∈[1,3])表示外力矢量,来自于无论是点接触还是线接触状态下的载荷;求解公式(21)中的系数矩阵χ中的Amn,Bmn和Cmn,即这些有接触点变形的相互影响系数,就可以确定齿轮接触柔性;
假定整个接触线分成Nc部分,每一部分可计算其分配的载荷,则大轮或小轮的解除柔性矩阵可写为
式中,Cis,js(is,js∈[1,Nc])为在某部分is的变形量,由于在某部分js施加载荷所致。
7.如权利要求6所述的基于双曲面壳单元模型的数值载荷齿面接触分析方法,其特征在于,针对步骤iiii.1),在求解齿面接触点过程中,为了保证足够的齿面点的匹配精度,分别考虑各自的齿形误差εi(i=1,2);首先,大小齿轮齿面点M2和M1各自先接触切平面投影,则在可行域内可以按照匹配规则进行目标点的匹配;此时,选择大轮齿面某离散点M2 (ij)(i∈[1,n],j∈[1,m])作为主动点,来匹配另一个小轮齿面的目标点;在匹配过程中,对于任一点M2 (ij),先根据精度要求设置一定的匹配可行域范围半径FRRr*,然后来搜寻可行域内另一齿面上的投影点,其目标投影点与主动投影点的距离作为匹配半径,可表示为
r=|rf2(φ2,θ2)-ε2nf2(φ2,θ2)-rf1(φ1,θ1)-ε1nf1(φ1,θ1)| (28)
式中,投影平面即齿轮接触切平面坐标系为Of(if,jf,kf),齿面点的参数表示其点矢rfi和法矢nfi,其中i=1时表示小轮,i=2时表示大轮;
齿面参数化离散点的精确匹配策略为:
1)对于任一点的匹配半径MR rk>r*,如果不匹配,则调整FRR再次搜寻目标点直至满足条件;
2)对于在可行域内的任意两个候选点,比较其MR,当出现rq<rp≤r*情况时,选取点M1 (q)为匹配目标点;最终,点M2 (ij)和点M1 (q)匹配成功,成为合格的匹配接触对M2 (ij)-M1 (q);式中的参数d表示实际接触过程中目标匹配点的距离。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811026119.0A CN109145484B (zh) | 2018-09-04 | 2018-09-04 | 基于双曲面壳单元模型的数值载荷齿面接触分析方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811026119.0A CN109145484B (zh) | 2018-09-04 | 2018-09-04 | 基于双曲面壳单元模型的数值载荷齿面接触分析方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109145484A CN109145484A (zh) | 2019-01-04 |
CN109145484B true CN109145484B (zh) | 2022-03-18 |
Family
ID=64826691
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811026119.0A Active CN109145484B (zh) | 2018-09-04 | 2018-09-04 | 基于双曲面壳单元模型的数值载荷齿面接触分析方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109145484B (zh) |
Families Citing this family (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109902377B (zh) * | 2019-02-25 | 2021-05-04 | 华中科技大学 | 一种间隙转动副接触应力分析的方法 |
CN109948279B (zh) * | 2019-03-29 | 2022-12-20 | 江苏精研科技股份有限公司 | 一种金属件整形的仿真设计方法 |
CN110096756B (zh) * | 2019-04-08 | 2020-01-17 | 河海大学 | 一种考虑荷载不确定性的自由曲面结构形态创建方法 |
CN110750922A (zh) * | 2019-09-11 | 2020-02-04 | 中南大学 | 含腹板结构的直齿轮有限元接触模型快速建模方法 |
CN110826158B (zh) * | 2019-10-28 | 2024-02-02 | 长安大学 | 一种基于啮入冲击最小的螺旋锥齿轮齿面Ease-off修形设计方法 |
CN111382503B (zh) * | 2020-02-27 | 2022-05-20 | 中南大学 | 弹性支承下旋转的柔性圆环的振动分析方法及*** |
CN112417593B (zh) * | 2020-11-18 | 2024-02-23 | 西北工业大学 | 一种航空发动机圆弧端齿连接装配安装角度的优化方法 |
CN117634057A (zh) * | 2023-10-18 | 2024-03-01 | 南京航空航天大学 | 一种含剥落故障的弧齿锥齿轮时变啮合刚度计算方法 |
CN117973151B (zh) * | 2024-03-27 | 2024-05-31 | 北京航空航天大学 | 一种压电层合板壳结构振动控制方法、***及设备 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106570278A (zh) * | 2016-11-08 | 2017-04-19 | 江苏大学 | 融合齿距偏差的圆柱齿轮啮合刚度计算方法 |
CN107577876A (zh) * | 2017-09-07 | 2018-01-12 | 清华大学 | 一种螺旋锥齿轮齿面加载性能多目标优化方法 |
CN108153981A (zh) * | 2017-12-26 | 2018-06-12 | 中航沈飞民用飞机有限责任公司 | 一种基于有限元分析的复合材料机身加筋壁板结构后屈曲分析方法 |
-
2018
- 2018-09-04 CN CN201811026119.0A patent/CN109145484B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106570278A (zh) * | 2016-11-08 | 2017-04-19 | 江苏大学 | 融合齿距偏差的圆柱齿轮啮合刚度计算方法 |
CN107577876A (zh) * | 2017-09-07 | 2018-01-12 | 清华大学 | 一种螺旋锥齿轮齿面加载性能多目标优化方法 |
CN108153981A (zh) * | 2017-12-26 | 2018-06-12 | 中航沈飞民用飞机有限责任公司 | 一种基于有限元分析的复合材料机身加筋壁板结构后屈曲分析方法 |
Non-Patent Citations (3)
Title |
---|
A novel operation approach to determine initial contact point for tooth contact analysis with errors of spiral bevel and hypoid gears;Han Ding et al.;《Mechanism and Machine Theory》;20161204;第155-170页 * |
An innovative determination approach to tooth compliance for spiral bevel and hypoid gears by using double-curved shell model and Rayleigh–Ritz approach;Han Ding et al;《Mechanism and Machine Theory》;20180827;第27-46页 * |
准双曲面齿轮动态啮合性能的有限元分析研究;唐进元 等;《振动与冲击》;20110725;第30卷(第7期);第101-105页 * |
Also Published As
Publication number | Publication date |
---|---|
CN109145484A (zh) | 2019-01-04 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109145484B (zh) | 基于双曲面壳单元模型的数值载荷齿面接触分析方法 | |
Vivet et al. | An analytical model for accurate and numerically efficient tooth contact analysis under load, applied to face-milled spiral bevel gears | |
CN109408857B (zh) | 螺旋锥齿轮形性协同制造的智能参数驱动模块化设计方法 | |
Tiwari et al. | Stress analysis of mating involute spur gear teeth | |
Sánchez et al. | Calculation of tooth bending strength and surface durability of internal spur gear drives | |
Barbieri et al. | Adaptive grid-size finite element modeling of helical gear pairs | |
Litvin et al. | New version of Novikov–Wildhaber helical gears: computerized design, simulation of meshing and stress analysis | |
CN102314534A (zh) | 一种基于振动可靠性和遗传算法的齿轮齿廓修形方法 | |
Zhu et al. | Pitch cone design and influence of misalignments on tooth contact behaviors of crossed beveloid gears | |
Dong et al. | Optimum design of the tooth root profile for improving bending capacity | |
Gonzalez-Perez et al. | Implementation of a finite element model for gear stress analysis based on tie-surface constraints and its validation through the Hertz's theory | |
Sanchez-Marin et al. | Numerical tooth contact analysis of gear transmissions through the discretization and adaptive refinement of the contact surfaces | |
Litvin et al. | Modified approach for tooth contact analysis of gear drives and automatic determination of guess values | |
Fuentes-Aznar et al. | Computational approach to design face-milled spiral bevel gear drives with favorable conditions of meshing and contact | |
Liu et al. | Parameterized High‐Precision Finite Element Modelling Method of 3D Helical Gears with Contact Zone Refinement | |
Ding et al. | Life cycle assessment-driven collaborative optimization model of power dry cutting for face-hobbing hypoid gear production | |
Litvin et al. | Generalized concept of meshing and contact of involute crossed helical gears and its application | |
Fang et al. | Interference-based technique for designing cutter flank using multiple radial infeed in gear skiving | |
Zhai et al. | A mathematical model for parametric tooth profile of spur gears | |
Cao et al. | Robust geometric parameter optimization of a crossed beveloid gear pair with approximate line contact | |
CN109033669B (zh) | 基于万能运动参数驱动的螺旋锥齿轮仿真加工建模方法 | |
Aziz et al. | Knowledge-based geometry generation for spur and helical gears | |
Zeng et al. | Optimized design of straight bevel gear tooth root transition surface | |
Lealin | Comparison and evaluation of classical methods of dimensional chains theory and their modern analogues | |
Kawalec et al. | Simulation of generation and tooth contact analysis of helical gears with crowned flanks |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CB03 | Change of inventor or designer information | ||
CB03 | Change of inventor or designer information |
Inventor after: Tang Jinyuan Inventor after: Ding Han Inventor after: Hu Zehua Inventor after: Wang Zhiwei Inventor before: Tang Jinyuan Inventor before: Ding Han |