CN109145484B - 基于双曲面壳单元模型的数值载荷齿面接触分析方法 - Google Patents

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),为了建立以精确的有限元模型，首先需要从几何层面分析其基本的结构特性，分别从齿高h_G、齿厚t_G及齿宽方向w_G对齿面进行三维尺寸定义；建立螺旋锥齿轮关于节锥母线上任一点P*的齿线方程：

式中，τ_G表示齿线的弧长，S_r为径向刀位，r_c为刀盘半径，R₁为外端锥距，R₂为内端锥距。

上述技术方案的一种实施方式中，针对步骤ii.1)，将壳单元定位为拥有不同厚度和高度的厚圆锥壳的一部分，其齿宽方向、齿高方向和齿厚方向近似于齿轮的三维尺寸方向；由于较小的半径与厚度及长度与厚度之比，螺旋锥齿轮的齿面弯曲行为可以利用剪切变形理论进行建模；建立局部主系统(α,β,z)，α和β分别表示中面的正交曲线主坐标，z表示法向厚度根据壳单元中面的不同几何形状，完成双曲面壳单元模型建模，双曲面壳单元模型同时兼容了圆柱壳与圆锥壳模型的特点。

上述技术方案的一种实施方式中，针对步骤ii.2)，具体如下：

在坐标系(α,β,z)中，U、V和W分别表示壳单元对应的三个主方向的变形位移量；在全局坐标系O(i₁,i₂,i₃)中，3D壳模型的一个任意点P_s可通过其位移矢量R(α,β,z)定义为

R(α,β,z):＝r(α,β)+zi_u(α,β)∈R³ (2)

式中，r(α,β)表示中面投影点P_s*的位移矢量，i_u(α,β)为向外的单位法向量，Z表示为点P_s到投影点P_s*的距离；由于目前的壳单元大小是有限的，壳单元被限制在以下的尺寸范围内

式中，h_G用来表示双曲面壳的全齿厚；

基于不同的几何表达，参考齿面的位移矢量可通过其沿全局坐标系的分量表示

r(α,β)＝r₁(α,β)i₁+r₂(α,β)i₂+r₃(α,β)i₃ (4)

式中，i₁，i₂，i₃分别是全局坐标系的单元矢量；此处，控制方程所需的拉曼参数通过参考曲面r(α,β)的第一类基本形式可以表示为

曲面的度量被定义为含有物理含义的数学项；第一类基本系数E、F和G，第二类基本系数L、M和N，通过在局部正交主坐标系(α,β,z)中曲率的参数线α和β的表达，目前存在关系F＝M＝0，则公式可简化为

式中，R_α(α,β)，R_β(α,β)为齿面的主曲率半径；表明了双曲面壳单元模型的意义，并说明了该模型可以足够精确表达齿面弯曲特性。

在双曲面壳单元的局部局部坐标系中，法向应变和剪切应变关系表示为

在该双曲面壳单元方程中，具有恒定曲率半径的参数化系数可写为

上述技术方案的一种实施方式中，针对步骤ii.3)，具体如下：

主曲率半径R_β表示刀盘切削轨迹上的一段圆弧

R_β:＝R_β(α,β)＝r_c (10)

齿高h_G和齿厚t_G分别假定线性变化且通过线性拟合可以表示其完整形状，则运用关于大端h₁和小端h₂的一个线性拟合来表示齿高

描述了沿齿宽方向的高度变化，其中τ_G可通过公式(1)确定。

曲率半径R_α表示从齿根圆R_α1到齿根圆的R_α2的不断变化的曲率半径，根据几何关系，任一点的锥角是恒定的，可以表示为

式中，X:＝G和X:＝P分别表示大轮和小轮，Σ表示两个齿轮轴的夹角，R_α的线性拟合为

从R_α1到R_α2逐渐减小，和齿厚拟合所呈现的线性关系一致。

在大端齿顶和齿根的齿厚参数t₁和t₂，在小端齿顶和齿根的齿厚参数t₃和t₄，可用来进行齿厚的线性拟合，有

或者，也可以表示为

这是由于在齿顶和齿根处具有同样的弯曲行为变化趋势所致。

对于厚度变化的双曲面壳单元模型，双曲面壳底部沿边缘固定而齿顶面、凸面、凹面、大端面和小断面都自由，可以表现螺旋锥齿轮齿面的弯曲行为特性；对于剪切旋转β_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的弹性应变能Q_SE，σ_x,、σ_θG和τ_xθG与均匀各向同质材料属性相关，其中E为杨氏模量，υ为泊松比，存在势能Q_PE为弹性应变能Q_SE和虚功Q_WF之差。假设Q_PE＝Q_SE-Q_WF，则假定势能的一阶差分为0，有

可采用成含未知系数的代数多项式的有限线性组合形式来表示横向偏移量W和剪切旋转量β_x及β_θG，必须满足预设的边界条件；

针对步骤iii)，关于确定螺旋锥齿轮啮合刚度，具体如下：

将边界约束条件代入位移假设(16)-(19)中，通过评法向和剪切应力来计算应变能。而势能Q_PE和虚功Q_WF则需要利用高斯-正交法来进行计算，获得线性方程组

式中，子矩阵K^mn(m,n∈[1,3])表示齿轮啮合刚度单元，可以根据材料属性相关的体积积分和预设的调和函数进行确定。而在螺旋锥齿轮接触传动过程中，整个啮合刚度K_TOTAL除了由(21)中所确定的大小轮的齿面啮合刚度K_P和K_G外，还包括Hertz接触刚度K_H，可得

K_TOTAL＝K_G+K_P+K_H (22)

而K_H根据Hetrz接触理论进行计算

式中κ₁和κ₂分别表示对应小轮和大轮的材料参数，D₁和D₂分别表示在啮合接触点位置两个齿面的曲率半径，而μ为泊松比，E为弹性模量；

在螺旋锥齿轮的载荷接触过程中，单对齿轮的啮合刚度确定如下：在接触坐标系O_G(X,Y,Z)，K₁和K₂分别表示有载荷也没有载荷的接触点，其位移矢量和接触力分别为μ₁、μ₂和F₁、F₂，则可得到下列关系

当多齿同时发生接触时，由于耦合作用的复杂啮合刚度可表示为

式中，p表示接触过程中发生接触的齿轮对的数目。

针对步骤iii)，关于螺旋锥齿轮接触柔性的确定，具体如下：

在公式(59)中，F^m(m∈[1，3])表示外力矢量，来自于无论是点接触还是线接触状态下的载荷；求解公式(59)中的系数矩阵χ中的A_mn,B_mn和C_mn，即这些有接触点变形的相互影响系数，就可以确定齿轮接触柔性；假定整个接触线分成N_c部分，每一部分可计算其分配的载荷，则大轮或小轮的解除柔性矩阵可写为

式中，C_is,js(i_s,j_s∈[1,N_c])为在某部分i_s的变形量，由于在某部分j_s施加载荷所致。

上述技术方案的一种实施方式中，针对步骤iiii.1)，在求解齿面接触点过程中，为了保证足够的齿面点的匹配精度，分别考虑各自的齿形误差ε_i(i＝1,2)；首先，大小齿轮齿面点M₂和M₁各自先接触切平面投影，则在可行域内可以按照匹配规则进行目标点的匹配；此时，选择大轮齿面某离散点M₂ ^(ij)(i∈[1,n],j∈[1,m])作为主动点，来匹配另一个小轮齿面的目标点；在匹配过程中，对于任一点M₂ ^(ij)，先根据精度要求设置一定的匹配可行域范围半径(FRR)r*，然后来搜寻可行域内另一齿面上的投影点，其目标投影点与主动投影点的距离作为匹配半径，可表示为

r＝|r_f2(φ₂,θ₂)-ε₂n_f2(φ₂,θ₂)-r_f1(φ₁,θ₁)-ε₁n_f1(φ₁,θ₁)| (28)

式中，投影平面即齿轮接触切平面坐标系为O_f(i_f，j_f，k_f)，齿面点的参数表示其点矢r_fi和法矢n_fi(i＝1表示小轮，i＝2表示大轮)。

齿面参数化离散点的精确匹配策略为：

1)对于任一点的匹配半径MR即r_k>r*，如果不匹配，则调整FRR再次搜寻目标点直至满足条件；

2)对于在可行域内的任意两个候选点，比较其MR，当出现r_q<r_p≤r*情况时，选取点M₁ ^(q)为匹配目标点；最终，点M₂ ^(ij)和点M₁ ^(q)匹配成功，成为合格的匹配接触对M₂ ^(ij)-M₁ ^(q)；式中的参数d表示实际接触过程中目标匹配点的距离。

上述技术方案的一种实施方式中，针对步骤iiii.2)，通过齿面接触点精确匹配策略得到的最终齿面接触对，可标记为M₁ ^K-M₂ ^K，而它们之间的初始距离假定为D_K，它们之间的弹性变形量ω_K(1)和ω_K(2)；它们在载荷变形情况下的相对距离为Θ；则在接触点位置，有如下的变形协调关系存在

ω_K(1)+ω_K(1)+D_K-Θ≥0(K＝1,2,…,N) (29)

式中，“＝”表示已经接触；而“>”表示尚未接触，其中ω_K(1)和ω_K(2)和前文所求的接触单元柔性C_Kj有关，可得

式中F_K表示目标匹配接触对M₁ ^K-M₂ ^K的接触。

上述技术方案的一种实施方式中，针对步骤iiii.3)，作用于接触区所有接触点的力之和F_K(K＝1,2,…,N)需要等于所有外部载荷，则存在以下受力平衡关系

{e}^T{F}＝P (32)

方程(32)和(29)是用来判定齿面是否接触的约束条件。两个齿面的弹性接触问题可以看成如何实现接触力F_K(K＝1,2,…,N)在一致接触柔性、外部载荷和变性相对距离的情况下，同时满足变形协调条件和受力平衡条件。

上述技术方案的一种实施方式中，针对步骤iiii.4)，具体如下：

由于方程(29)是一个大于或等于0的不等式约束，引入松弛变量Y，则可以变成等式约束，为

[C]{F}+{D}-Θ{e}-[I]{Y}＝{0} (33)

式中，

联立方程(32)和(33)，共得到N+1个标量方程；要求解未知参数(F，D，Y)，F_K≥0，Θ≥0，Y_K≥0，并要求满足接触条件(29)，即存在如下条件

Y_K＝0,F_K≥0；Y_K＝0,F_K≥0 (34)

利用改进单纯形法来求解满足接触条件(29)的方程(32)和(33)，则该齿面弹性接触问题的目标函数表示为一般线性规划：

式中，对于N+1个标量方程需要人为引入一系列非负的变量参数Z＝[Z₁，Z₂，…，Z_N，Z_N+1]，使得等式约束两边的常数项都为非负值；

考虑接触条件(29)，采用改进单纯形法，假定目标函数(34)中Z₀＝-Z_N+2，有

将上式作为一个约束条件，Z_N+2作为固定的基变量，用分离系数法将目标函数(34)写为单纯形表，得到变换后的标准形式。

上述技术方案的一种实施方式中，针对步骤iiii.5)，具体如下：

首先，通过基于双曲面壳模型和Rayleigh-Ritz法确定齿面接触柔性系数C_P,G；

然后，利用改进单纯形法的线性规划方法，求解齿面载荷分配和载荷传动误差及其最大载荷传动误差LTE_MAX；

其次，根据载荷分配设置载荷边界条件，求解最大齿面接触应力CP_MAX；

最后，根据啮合周期内的接触应力发生的时间历程计算出齿面重合度ε_γ。

本发明为了建立精确的有限元模型，从几何层面分析螺旋锥齿轮齿面基本的结构特性，进行齿面弯曲行为分析；根据齿面弯曲行为分析，结合以往齿轮有限元建模的基本理论和方法，提出了螺旋锥齿轮的双曲面壳单元的建模方法，既考虑齿轮啮合刚度又考虑齿轮接触柔性，以保证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给出了螺旋锥齿轮齿面形状的基本示意图，分别从齿高h_G、齿厚t_G及齿宽方向w_G对齿面进行三维尺寸定义。尤其是对于端面铣削(face-milling)的收缩齿，齿高从小端到大端逐渐增加，齿根厚度从小端到大端的逐渐增加，齿厚在齿宽任意点处沿齿高逐渐减小，这意味着精确的齿轮模型，必须考虑以下弯曲特性：齿宽方向的刚度变化；齿面纵向的曲率变化；齿向螺旋且逐渐收缩变化(格里森齿制)。

众所周知，螺旋锥齿轮从在节锥母线方向看是一条曲线，而其相对于节锥母线的弯曲程度通常用在齿轮节锥曲线上的螺旋角β_P来定义。

图2给出了螺旋锥齿轮最常见的两种加工方式，即单分度的端面铣削，主要是适用于Gleason收缩齿，和连续分度的端面滚切，主要适用于Klingelnberg和Oerlikon等高齿。对于端面铣削加工，齿高是变化的，齿面纵向曲率是一段圆弧，其曲率半径等于刀盘的曲率半径，齿轮被刀盘刀具的内外切削刃分别单独加工。对于端面滚切加工齿轮，齿高是恒定的，纵向齿线是一个延伸的外摆线，齿轮的凹面P₁面和凸面P₂同时被刀具内外切削刃加工。中点P的螺旋角定义为β_P。在△OO_cP，根据余弦定理，有以下关系

式中R为节锥，r_c为刀盘半径，S_r为径向刀位。假设P*为节线任一点且其锥距为R_X，其螺旋角β_X根据关于△OO_cP*的余弦定理可得

联立以上两个方程(1)和(1-1)，则节锥线上任意一点的螺旋角表示为

因此，参数r_c和β_X则可参数定义螺旋锥齿轮的齿线。计算齿线的弧长τ_G，如图3所示。

在半径为r_c的圆O_c上，P_AO_c＝P_BO_c＝r_c，有以下关系

τ_G＝∠P_AO_cP_B·r_c＝(∠P_AO_cO-∠P_BO_cO)·r_c (1-3)

在△P_AO_cO和△P_BO_cO，O_cO＝S_r，P_AO＝R₂为内端锥距，P_AO＝R₁为外端锥距。根据余弦定理可得：

联立以上两个方程(1-4,1-5)有

这表明齿线是恒定的，且可以被用来确定所提出的有限元壳单元模型的几何尺寸。由此可得，壳单元比一般的梁单元或者板单元更能精确的表达螺旋锥齿轮弯曲行为特性。

关于ii.1)建立精确的双曲面壳单元模型

精确的有限元模型是齿轮接触力学性能求解的关键步骤。有限元法方法就是把复杂的整体结构离散到有限个单元(Finite Element)，再把这种理想化的假定和力学控制方程施加于结构内部的每一个单元，然后通过单元分析组装得到结构总刚度方程，再通过边界条件和其他约束解得结构总反应。

因此，数值载荷齿面接触分析(NLTCA)作为分析齿轮接触力学性能的重要手段，其齿轮有限元精确模型的建立也十分重要。本发明根据前述螺旋锥齿轮齿面弯曲行为特性分析，结合以往齿轮有限元建模的基本理论和方法，提出了螺旋锥齿轮双曲面壳单元的建模方法，具体如下：

ii.1.1)建立精确的双曲面壳模型

目前，螺旋主齿轮的有限元单元建模研究，经历了一个逐渐由浅入深的过程：从简化的悬臂梁单元到收缩板单元到目前的圆柱壳单元。区别以往有限元模型结构，本发明提出更为精确的双曲面壳单元模型，如图4所示。壳单元可以看成拥有不同厚度和高度的厚圆锥壳的一部分，其齿宽方向、齿高方向和齿厚方向近似于齿轮的三维尺寸方向。由于较小的半径与厚度及长度与厚度之比，螺旋锥齿轮的齿面弯曲行为可以利用剪切变形理论进行建模。特别地，在主要的集中壳单元模型中，双曲面壳单元模型被选用来完成螺旋锥齿轮的有限元建模，因为它能更加精确的反映出齿面的弯曲特性。

在图4所示的模型中，双曲面壳单元模型同时兼容了圆柱壳单元与圆锥壳单元模型的特点。建立局部主系统(α,β,z)，α和β分别表示中面的正交曲线主坐标，z表示法向厚度。由于端面铣削的齿面具有更复杂的弯曲特性，且厚度可变，本发明以此为例进行论述。而对于端面滚切等高齿面模型，只需要设定法向厚度恒定即可。值得注意的是，为了探讨给定齿面上逐点变化的曲率半径的双曲面几何特性，本发明给出壳单元中面的不同几何形状，以完成快速可靠且简便的有限元结构建模。

ii.2)建立双曲面壳单元方程

前面已经从螺旋锥齿轮几何层面分析其基本的结构特性，它在齿高、齿宽和齿厚方向的三维尺寸都是变化的，这也符合壳单元模型中的双曲面壳单元模型的几何特性。因此，考虑齿面弯曲特性的螺旋锥齿轮双曲面壳模型的几何特点，图5给出了其基本参数化描述。在坐标系(α,β,z)中，U、V和W分别表示了壳单元对应的三个主方向的变形位移量。在全局坐标系O(i₁,i₂,i₃)中，3D壳模型的一个任意点P_s可通过其位移矢量R(α,β,z)定义为

R(α,β,z):＝r(α,β)+zi_u(α,β)∈R³ (2)

式中，r(α,β)表示壳单元中面投影点P_s*的位矢，i_u(α,β)为向外的单位法向量，Z表示为点P_s到投影点P_s*的距离。由于目前的壳单元大小是有限的，被限制在以下的尺寸范围内

式中，h_G表示双曲面壳的全齿厚。

基于不同的几何表达，参考齿面的位移矢量可通过其沿全局坐标系的分量表示

r(α,β)＝r₁(α,β)i₁+r₂(α,β)i₂+r₃(α,β)i₃ (4)

式中，i₁，i₂，i₃分别是全局坐标系的单元矢量。此处，控制方程所需的拉曼参数通过参考曲面r(α,β)的第一类基本形式可以表示为

曲面的度量被定义为含有物理含义的数学项。第一类基本系数E、F和G可构成一个矩阵C_P

可以进行倒置，由于其行列式不为0。通过公式(5)和(5-1)可得，

由此，易于计算向外单位法矢为

同理，包含第二类基本系数L、M和N矩阵CC_P为

联立(3)和(5)，有一个诸如下三阶处理Θ_P

包含齿面的曲率，其确定允许定义高斯曲率Ω_G和平均曲率Ω_M。故有以下关系

通过求解，在α和β方向的主曲率k_α(α,β)和k_β(α,β)为

式中，R_α(α,β)，R_β(α,β)为齿面的主曲率半径。通过在局部正交主坐标系(α,β,z)中曲率的参数线α和β的表达，目前存在关系F＝M＝0，则公式(5-7)可简化为

表明了双曲面壳单元模型的意义，并说明了该模型可以足够精确表达螺旋锥齿面的弯曲特性。

在双曲面壳的局部坐标系中，三维应变分量表示为

式中，D_z和

分别表示微分处理，包含和z相关项及与主局部坐标相关曲率项

另外，便于论述，引入如下表达

ε:＝ε(α,β,z,t)＝[ε_αα ε_ββ ε_zz γ_αβ γ_αz γ_βz]^T (6-3)

U＝[U(α,β,t)V(α,β,t)W(α,β,t)]^T (6-4)

式中，ε包含法向应变分量ε_ααε_ββ,ε_zz和剪切应变分量γ_αβ,γ_αz,γ_βz，U包含位于壳中面(z＝0)的点的一般位移分量，t表示时间变量，则可进行二阶处理

联立以上(6-1)-(6-5)，法向应变和剪切应变关系表示为

在该双曲面壳方程中，具有恒定曲率半径的参数化系数可写为

关于ii.3)确定双曲面壳单元模型的基本几何定义和边界约束条件，具体如下：

图6给出了建立的双曲面壳模型的基本几何定义和边界约束条件。其中，β＝θ_G表示任意齿面接触点从刀盘中心到小端面之间的角度，α＝x则表示在特定

接触点的齿高。参考图3，主曲率半径R_β表示刀盘切削轨迹上的一段圆弧，

R_β:＝R_β(α,β)＝r_c (10)

在该模型中，齿高h_G和齿厚t_G分别假定线性变化且通过线性拟合可以表示其完整形状。则运用关于大端h₁和小端h₂的一个线性拟合来表示齿高为，

描述了沿齿宽方向的高度变化，其中τ_G可通过公式(1)确定。

而另一个表示曲率半径R_α，表示从齿根圆R_α1到齿根圆的R_α2的不断变化的曲率半径。根据几何关系，任一点的锥角是恒定的，可以表示为

式中，X:＝G和X:＝P分别表示大轮和小轮，Σ表示两个齿轮轴的夹角。由此，R_α的线性拟合表示为

从R_α1到R_α2逐渐减小，和齿厚拟合所呈现的线性关系一致。

根据图6的几何定义，在大端齿顶和齿根的齿厚参数t₁和t₂，在小端齿顶和齿根的齿厚参数t₃和t₄，可用来进行齿厚的线性拟合，有

或者，也可以表示为

这是由于在齿顶和齿根处具有同样的弯曲行为变化趋势所致。

对于厚度变化的双曲面壳单元模型，图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的弹性应变能Q_SE为

σ_x,、σ_θG和τ_xθG与均匀各向同质材料属性相关，可表示为

式中E为杨氏模量，υ为泊松比。根据公式(19-3)-(19-12)，存在

由于导致可变性的外力F_p所做的功Q_WF为

势能Q_PE为弹性应变能Q_SE和虚功Q_WF之差。假设Q_PE＝Q_SE-Q_WF，则假定势能的一阶差分为0，有

可采用成含未知系数的代数多项式的有限线性组合形式来表示横向偏移量(W)和剪切旋转量β_x及β_θG，必须满足预设的边界条件。

iii.2)关于啮合刚度的确定

将边界约束条件(16)-(19)代入位移假设(19-6)-(19-8)中，通过评法向和剪切应力来计算应变能(56)。而势能Q_PE和虚功Q_WF则需要利用高斯-正交法来进行计算，获得线性方程组

式中，子矩阵K^mn(m,n∈[1,3])表示齿轮啮合刚度单元，可以根据材料属性相关的体积积分和预设的调和函数进行确定。而在螺旋锥齿轮接触传动过程中，整个啮合刚度K_TOTAL除了由(59)中所确定的大小轮的齿面啮合刚度K_P和K_G外，还包括Hertz接触刚度K_H，可得

K_TOTAL＝K_G+K_P+K_H (22)

而K_H根据Hetrz接触理论进行计算

式中κ₁和κ₂分别表示对应小轮和大轮的材料参数，D₁和D₂分别表示在啮合接触点位置两个齿面的曲率半径，而μ为泊松比，E为弹性模量。

在螺旋锥齿轮的载荷接触过程中，单对齿的啮合刚度如图7所示。在接触坐标系O_G(X,Y,Z)，K₁和K₂分别表示有载荷也没有载荷的接触点，其位移矢量和接触力分别为μ₁、μ₂和F₁、F₂，则可得到下列关系

当多齿同时发生接触时，由于耦合作用的复杂啮合刚度可表示为

式中，p表示接触过程中发生接触的齿轮对的数目。

iii.2)关于齿轮接触柔性的确定

在公式(21)中，F^m(m∈[1,3])表示外力矢量，来自于无论是点接触还是线接触状态下的载荷。求解公式(21)中的系数矩阵χ中的A_mn,B_mn和C_mn，即这些有接触点变形的相互影响系数，就可以确定齿轮接触柔性。

由于每一对接触点位置都要进行计算，关于这些变形函数的未知系数的求解是十分耗时的。假定整个接触线分成N_c部分，每一部分可计算其分配的载荷，则大轮或小轮的接触柔性矩阵可写为

式中，C_is,js(i_s,j_s∈[1,N_c])为在某部分i_s的变形量，由于在某部分j_s施加载荷所致。

关于iiii)NLTCA的精确计算

基于精确有限元建模和啮合刚度与接触柔性矩阵的确定，就可以进行NLTCA的数值精确求解。其过程主要包括以下几个关键步骤：1)确定精确的齿面接触对；2)建立考虑载荷变形情况下的变形协调关系；3)建立考虑接触约束条件的受力平衡方程；4)基于改进单纯形法的NLTCA的线性规划。改进单纯形法是1953年美国数学家G.B.丹齐克提出了改进单纯形法，其基本步骤和单纯形法大致相同，主要区别是在逐次迭代中不再以高斯消去法为基础，而是由旧基阵的逆去直接计算新基阵的逆，再由此确定检验数。这样做可以减少迭代中的累积误差，提高计算精度，同时也减少了在计算机上的存储量。

步骤iiii)包括以下分步：

iiii.1)齿面接触点的精确匹配策略

本申请人已经就含误差齿面接触分析(eTCA)的求解申请了专利，其中公开了相应的TCA流程和接触点确定方法。本发明只重点论述NLTCA中的齿面接触点的精确匹配问题。在求解齿面接触点过程中，为了保证足够的齿面点的匹配精度，考虑了分别考虑各自的齿形误差ε_i(i＝1,2)。首先，大小齿轮齿面点M₂和M₁各自先接触切平面投影，则在可行域内可以按照一定的匹配规则进行目标点的匹配。此时，选择大轮齿面某离散点M₂ ^(ij)(i∈[1,n],j∈[1,m])作为主动点，来匹配另一个小轮齿面的目标点。在匹配过程中，对于任一点M₂ ^(ij)，先根据精度要求设置一定的匹配可行域范围半径(FRR)r*，然后来搜寻可行域内另一齿面上的投影点。如图8所示，其目标投影点与主动投影点的距离作为匹配半径，可表示为

r＝|r_f2(φ₂,θ₂)-ε₂n_f2(φ₂,θ₂)-r_f1(φ₁,θ₁)-ε₁n_f1(φ₁,θ₁)| (28)

式中，投影平面即齿轮接触切平面坐标系为O_f(i_f，j_f，k_f)，齿面点的参数表示其点矢r_fi和法矢n_fi(i＝1表示小轮，i＝2表示大轮)。

因此，如图9所示，齿面参数化离散点的精确匹配策略为：

i)对于任一点的MR即r_k>r*，如果不匹配，则调整FRR再次搜寻目标点直至满足条件；

ii)对于在可行域内的任意两个候选点，比较其MR，当出现r_q<r_p≤r*情况时，选取点M₁ ^(q)为匹配目标点。最终，点M₂ ^(ij)和点M₁ ^(q)匹配成功，成为合格的匹配接触对M₂ ^(ij)-M₁ ^(q)。

很明显，对于两个候选点，存在完全重合关系n_f1 ^(ij)＝n_f2 ^(ij)即r*＝0这种情况是最优最有效解，但是这只是特例。另外，如果当前匹配精度不符合要求或者需要更高精度的匹配，本文提供了更为优化的匹配策略。当在可行域内不存在一个匹配点时，FRR可以增加进行再次匹配。只需要在上次未成功匹配点附近区域以较小的MR进行搜索。但是，为了提高匹配精度，搜索半径可以逐次增大。因此，可以再上一次目标区域进行细分，搜寻满足条件的目标点。基本的细分规则为，以上次候选点为中心进行一个简单的5×5网格等分；如果失败，则再次进行(4n+1)×(4n+1)等分的细分(n为细分次数)，直至满足条件。

最后，所有的目标匹配点被确定，可进行精确编号排序，如图10所示。首先选取大小轮齿面点作为输入，并从点M_i ^(1-1)到M_i ^(n-m)进行精确编号，尤其是某一行的末尾点和下一行的起始点要注意衔接，诸如M_i ^(1-m)和M_i ^(2-1)；然后，进行有效的匹配；每一对目标匹配点按照原来的序列号进行保存。通过匹配和排序后，下列的优化关系可用来求解对应的齿面接触点

mind＝|r_f1(φ₁,θ₁)-r_f2(φ₂,θ₂)| (28-1)

式中的参数d表示实际接触过程中目标匹配点的距离。

iiii.2)变形协调条件

通过齿面接触点的精确匹配策略得到最终的齿面接触对，可标记为M₁ ^K-M₂ ^K，而它们之间的初始距离假定为D_K，它们之间的弹性变形量ω_K(1)和ω_K(2)；它们在载荷变形情况下的相对距离为Θ。则在接触点位置，有如下的变形协调关系存在

ω_K(1)+ω_K(1)+D_K-Θ≥0(K＝1,2,…,N) (29)

式中，“＝”表示已经接触；而“>”表示尚未接触，其中ω_K(1)和ω_K(2)和前文所求的接触单元柔性C_Kj有关，可得

式中F_K表示目标匹配接触对M₁ ^K-M₂ ^K的接触力。联立(68)-(69)可以得到：

[C]{F}+{D}-Θ{e}≥{0} (30-2)

式中，有

[C]＝[C_Kj]＝[C_K(1)j+C_K(2)j]，{F}＝{F₁,F₂,…,F_K,…,F_N}^T

{D}＝{D₁,D₂,…,D_K,…,D_N}^T，{e}＝{1,1,…,1,…,1}^T，{0}＝{0,0,…,0,…,0}^T

iiii.3)受力平衡条件

作用于接触区所有接触点的力之和F_K(K＝1,2,…,N)需要等于所有外部载荷，则存在以下受力平衡关系

{e}^T{F}＝P (32)

方程(32)和(30-2)是用来判定齿面是否接触的约束条件。两个齿面的弹性接触问题可以看成如何实现接触力F_K(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)，F_K≥0，Θ≥0，Y_K≥0，并要求满足接触条件(29)，即存在如下条件

Y_K＝0,F_K≥0；Y_K＝0,F_K≥0 (34)

本发明利用改进单纯形法来求解满足接触条件(29)的方程(32)和(33)，则该齿面弹性接触问题的目标函数表示为一般线性规划：

式中，对于N+1个标量方程需要人为引入一系列非负的变量参数Z＝[Z₁，Z₂，…，Z_N，Z_N+1]，使得等式约束两边的常数项都为非负值。

在本发明所提出的螺旋锥齿轮NLTCA求解过程中，与一般的线性规划不同，还有考虑接触条件(29)，采用改进的单纯形法，假定目标函数(34)中Z₀＝-Z_N+2，有

将上式作为一个约束条件，Z_N+2作为固定的基变量，用分离系数法将目标函数(34)写为单纯形表形式，采用改进单纯形法进行求解。

iiii.5)NLTCA评价项的确定

图11给出了确定NLTCA评价项的基本流程，其关键点为：1)有限元建模；2)啮合刚度、解除柔性求解；3)NLTCA方程设置及其求解。

首先，通过基于双曲面壳模型和Rayleigh-Ritz法可以确定齿面接触柔性系数C_P,G；

然后利用改进单纯形法的线性规划方法，就可以求解齿面载荷分配和载荷传动误差及其最大载荷传动误差LTE_MAX；

其次，根据载荷分配设置载荷边界条件，就可以求解最大齿面接触应力CP_MAX；

最后，根据啮合周期内的接触应力发生的时间历程就可以计算出齿面重合度ε_γ；

很显然，该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)为大轮M₂ ^(8,10)和小轮M₁ ^(7,11)，而FRR为0.25mm，计算的MR为0.173654。在表7更高精度匹配中，图12给出了具体的5×5细分示意图。设定一个新的FRR，以上一次的目标匹配点M₁ ^(7,11)为中心进行5×5细分，搜寻到最后的目标匹配点为M₁ ^(7,11) _-1 ^(4,4)。通过以上两种匹配，比较了其匹配之后求解最小距离，一次来验证其正确性。其中，目标匹配点M₂ ^(8,10)-M₁ ^(7,11)的为0.1683314mm，而M₂ ^(8,10)-M₁ ^(7,11) _-1 ^(4,4)的为0.0267345mm。因此，根据给出的齿面接触点匹配策略，最终会选取接触对M₂ ^(8,10)-M₁ ^{(7 ,11)} _-1 ^(4,4)作为结果输出来求解精确初始点。

关于双曲面壳单元和Rayleigh-Ritz法，具体如下：

一副端面铣削加工的准双曲面齿轮齿面基本参数如表8和表9所示。其中，准双曲面齿轮的HGT调整卡中，大轮采用展成法加工，小轮采用刀倾法加工，选择机械式的NO.116机床采用较为普遍的五刀法加工。基于提供的设计参数就可以齿轮副建模，用于进行有限元单元设定和双曲面壳方程求解。这里，为了来验证本文提出的基于双曲面壳理论的Rayleigh-Ritz方法的可行性，选取了NLTCA中的重要的齿面接触刚度值作为基本指标，以成熟的商用有限元分析(FEA)软件例如

的LTCA计算结果为参考，进行比较评价。

表8齿坯设计参数

表9 HGT加工调整卡参数

在NLTCA设置中，大轮凸面和凹面被设定为接触齿面，基本的材料属性为杨氏模量E＝2.09×10⁵Mpa，泊松比υ＝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×10⁶N/mm，最小值为6.869×10⁵N/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)确定变形协调条件

通过齿面接触点精确匹配策略得到的最终齿面接触对，可标记为M₁ ^K-M₂ ^K，而它们之间的初始距离假定为D_K，它们之间的弹性变形量ω_K(1)和ω_K(2)；它们在载荷变形情况下的相对距离为Θ；则在接触点位置，有如下的变形协调关系存在

ω_K(1)+ω_K(1)+D_K-Θ≥0 (K＝1,2,…,N) (29)

式中，“＝”表示已经接触；而“＞”表示尚未接触，其中ω_K(1)和ω_K(2)和前文所求的接触单元柔性C_Kj有关，可得

式中F_K表示目标匹配接触对M₁ ^K-M₂ ^K的接触；

iiii.3)确定受力平衡条件

作用于接触区所有接触点的力之和F_K(K＝1,2,…,N)需要等于所有外部载荷，则存在以下受力平衡关系

{e}^T{F}＝P (32)

方程(32)和(29)是用来判定齿面是否接触的约束条件； 两个齿面的弹性接触问题可以看成如何实现接触力F_K(K＝1,2,…,N)在一致接触柔性、外部载荷和变性相对距离的情况下，同时满足变形协调条件和受力平衡条件；

iiii.4)确定求解数值载荷齿面接触分析问题的线性规划方法；

由于方程(29)是一个大于或等于0的不等式约束，引入松弛变量Y，则可以变成等式约束，为

[C]{F}+{D}-Θ{e}-[I]{Y}＝{0} (33)

式中，

{Y}＝{Y₁,Y₂,…,Y_K,…,Y_N}^T

联立方程(32)和(33)，共得到N+1个标量方程；要求解未知参数(F，D，Y)，F_K≥0，Θ≥0，Y_K≥0，并要求满足接触条件(29)，即存在如下条件

Y_K＝0,F_K≥0；Y_K＝0,F_K≥0 (34)

利用改进单纯形法来求解满足接触条件(29)的方程(32)和(33)，则该齿面弹性接触问题的目标函数表示为一般线性规划：

式中，对于N+1个标量方程需要人为引入一系列非负的变量参数Z＝[Z₁，Z₂，…，Z_N，Z_N+1]，使得等式约束两边的常数项都为非负值；

考虑接触条件(29)，采用改进单纯形法，假定目标函数(34)中Z₀＝-Z_N+2，有

将上式作为一个约束条件，Z_N+2作为固定的基变量，用分离系数法将目标函数(34)写为单纯形表；在横线下面的每一行代表一个约束条件，共有N+2个，竖线右边为原等约束条件(34)的右边常数项；将第N+2行减去前N+1行的和，并表示在第N+2行中，得到变换后的标准形式；

iiii.5)确定数值载荷齿面接触分析评价项：

首先，通过基于双曲面壳模型和Rayleigh-Ritz法确定齿面接触柔性系数C_P,G；

然后，利用改进单纯形法的线性规划方法，求解齿面载荷分配和载荷传动误差及其最大载荷传动误差LTE_MAX；

其次，根据载荷分配设置载荷边界条件，求解最大齿面接触应力CP_MAX；

最后，根据啮合周期内的接触应力发生的时间历程计算出齿面重合度ε_γ。

2.如权利要求1所述的基于双曲面壳单元模型的数值载荷齿面接触分析方法，其特征在于：针对步骤i),从几何层面分析其基本的结构特性，分别从齿高h_G、齿厚t_G及齿宽方向w_G对齿面进行三维尺寸定义；建立螺旋锥齿轮关于节锥母线上任一点P*的齿线方程：

式中，τ_G表示齿线的弧长，S_r为径向刀位，r_c为刀盘半径，R₁为外端锥距，R₂为内端锥距。

3.如权利要求2所述的基于双曲面壳单元模型的数值载荷齿面接触分析方法，其特征在于：针对步骤ii.1)，将壳单元定位为拥有不同厚度和高度的厚圆锥壳的一部分，其齿宽方向、齿高方向和齿厚方向近似于齿轮的三维尺寸方向；由于较小的半径与厚度及长度与厚度之比，螺旋锥齿轮的齿面弯曲行为利用剪切变形理论进行建模；建立局部主系统(α,β,z)，α和β分别表示中面的正交曲线主坐标，z表示法向厚度根据壳单元中面的不同几何形状，完成双曲面壳单元模型建模，双曲面壳单元模型同时兼容了圆柱壳与圆锥壳模型的特点。

4.如权利要求3所述的基于双曲面壳单元模型的数值载荷齿面接触分析方法，其特征在于，针对步骤ii.2)，具体如下：

在坐标系(α,β,z)中，U、V和W分别表示壳单元对应的三个主方向的变形位移量；在全局坐标系O(i₁,i₂,i₃)中，3D壳模型的一个任意点P_s可通过其位移矢量R(α,β,z)定义为

R(α,β,z):＝r(α,β)+zi_u(α,β)∈R³ (2)

式中，r(α,β)表示中面投影点P_s*的位移矢量，i_u(α,β)为向外的单位法向量，Z表示为点P_s到投影点P_s*的距离；由于目前的壳单元大小是有限的，壳单元被限制在以下的尺寸范围内

式中，h_G用来表示双曲面壳的全齿厚；

基于不同的几何表达，参考齿面的位移矢量可通过其沿全局坐标系的分量表示

r(α,β)＝r₁(α,β)i₁+r₂(α,β)i₂+r₃(α,β)i₃ (4)

式中，i₁，i₂，i₃分别是全局坐标系的单元矢量；此处，控制方程所需的拉曼参数通过参考曲面r(α,β)的第一类基本形式可以表示为

曲面的度量被定义为含有物理含义的数学项；第一类基本系数E、F和G，第二类基本系数L、M和N，通过在局部正交主坐标系(α,β,z)中曲率的参数线α和β的表达，目前存在关系F＝M＝0，则公式可简化为

式中，R_α(α,β)，R_β(α,β)为齿面的主曲率半径；表明了双曲面壳单元模型的意义，并说明了该模型可以足够精确表达齿面弯曲特性；

在双曲面壳单元的局部坐标系中，法向应变和剪切应变关系表示为

在该双曲面壳单元方程中，具有恒定曲率半径的参数化系数可写为

5.如权利要求1所述的基于双曲面壳单元模型的数值载荷齿面接触分析方法，其特征在于，针对步骤ii.3)，具体如下：

主曲率半径R_β表示刀盘切削轨迹上的一段圆弧

R_β:＝R_β(α,β)＝r_c (10)

齿高h_G和齿厚t_G分别假定线性变化且通过线性拟合可以表示其完整形状，则运用关于大端h₁和小端h₂的一个线性拟合来表示齿高

描述了沿齿宽方向的高度变化，其中τ_G可通过公式(1)确定；

曲率半径R_α表示从齿根圆R_α1到齿根圆的R_α2的不断变化的曲率半径，根据几何关系，任一点的锥角是恒定的，可以表示为

式中，X:＝G和X:＝P分别表示大轮和小轮，Σ表示两个齿轮轴的夹角；由此，R_α的线性拟合表示为

从R_α1到R_α2逐渐减小，和齿厚拟合所呈现的线性关系一致；

在大端齿顶和齿根的齿厚参数t₁和t₂，在小端齿顶和齿根的齿厚参数t₃和t₄，可用来进行齿厚的线性拟合，有

或者，也可以表示为

这是由于在齿顶和齿根处具有同样的弯曲行为变化趋势所致；

对于厚度变化的双曲面壳单元模型，双曲面壳底部沿边缘固定而齿顶面、凸面、凹面、大端面和小断面都自由，可以表现螺旋锥齿轮齿面的弯曲行为特性；对于剪切旋转β_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的弹性应变能Q_SE，σ_x,、σ_θG和τ_xθG与均匀各向同质材料属性相关，其中E为杨氏模量，υ为泊松比，存在势能Q_PE为弹性应变能Q_SE和虚功Q_WF之差； 假设Q_PE＝Q_SE-Q_WF，则假定势能的一阶差分为0，有

可采用成含未知系数的代数多项式的有限线性组合形式来表示横向偏移量W和剪切旋转量β_x及β_θG，必须满足预设的边界条件；

针对步骤iii)，关于确定螺旋锥齿轮啮合刚度，具体如下：

将边界约束条件代入位移假设(16)-(19)中，通过评法向和剪切应力来计算应变能；而势能Q_PE和虚功Q_WF则需要利用高斯-正交法来进行计算，获得线性方程组

式中，子矩阵K^mn(m,n∈[1,3])表示齿轮啮合刚度单元，可以根据材料属性相关的体积积分和预设的调和函数进行确定； 而在螺旋锥齿轮接触传动过程中，整个啮合刚度K_TOTAL除了由(21)中所确定的大小轮的齿面啮合刚度K_P和K_G外，还包括Hertz接触刚度K_H，可得

K_TOTAL＝K_G+K_P+K_H (22)

而K_H根据Hetrz接触理论进行计算

式中κ₁和κ₂分别表示对应小轮和大轮的材料参数，D₁和D₂分别表示在啮合接触点位置两个齿面的曲率半径，而μ为泊松比，E为弹性模量；

在螺旋锥齿轮的载荷接触过程中，单对齿轮的啮合刚度确定如下：在接触坐标系O_G(X,Y,Z)，K₁和K₂分别表示有载荷也没有载荷的接触点，其位移矢量和接触力分别为μ₁、μ₂和F₁、F₂，则可得到下列关系

当多齿同时发生接触时，由于耦合作用的复杂啮合刚度可表示为

式中，p表示接触过程中发生接触的齿轮对的数目；

针对步骤iii)，关于螺旋锥齿轮接触柔性的确定，具体如下：

在公式(21)中，F^m(m∈[1，3])表示外力矢量，来自于无论是点接触还是线接触状态下的载荷；求解公式(21)中的系数矩阵χ中的A_mn,B_mn和C_mn，即这些有接触点变形的相互影响系数，就可以确定齿轮接触柔性；

假定整个接触线分成N_c部分，每一部分可计算其分配的载荷，则大轮或小轮的解除柔性矩阵可写为

式中，C_is,js(i_s,j_s∈[1,N_c])为在某部分i_s的变形量，由于在某部分j_s施加载荷所致。

7.如权利要求6所述的基于双曲面壳单元模型的数值载荷齿面接触分析方法，其特征在于，针对步骤iiii.1)，在求解齿面接触点过程中，为了保证足够的齿面点的匹配精度，分别考虑各自的齿形误差ε_i(i＝1,2)；首先，大小齿轮齿面点M₂和M₁各自先接触切平面投影，则在可行域内可以按照匹配规则进行目标点的匹配；此时，选择大轮齿面某离散点M₂ ^(ij)(i∈[1,n],j∈[1,m])作为主动点，来匹配另一个小轮齿面的目标点；在匹配过程中，对于任一点M₂ ^(ij)，先根据精度要求设置一定的匹配可行域范围半径FRRr*，然后来搜寻可行域内另一齿面上的投影点，其目标投影点与主动投影点的距离作为匹配半径，可表示为

r＝|r_f2(φ₂,θ₂)-ε₂n_f2(φ₂,θ₂)-r_f1(φ₁,θ₁)-ε₁n_f1(φ₁,θ₁)| (28)

式中，投影平面即齿轮接触切平面坐标系为O_f(i_f，j_f，k_f)，齿面点的参数表示其点矢r_fi和法矢n_fi，其中i＝1时表示小轮，i＝2时表示大轮；

齿面参数化离散点的精确匹配策略为：

1)对于任一点的匹配半径MR r_k＞r*，如果不匹配，则调整FRR再次搜寻目标点直至满足条件；

2)对于在可行域内的任意两个候选点，比较其MR，当出现r_q<r_p≤r*情况时，选取点M₁ ^(q)为匹配目标点；最终，点M₂ ^(ij)和点M₁ ^(q)匹配成功，成为合格的匹配接触对M₂ ^(ij)-M₁ ^(q)；式中的参数d表示实际接触过程中目标匹配点的距离。

Images