CN106844818A - 基于粗糙表面的直齿轮三维接触刚度计算方法 - Google Patents

Abstract

基于粗糙表面的直齿轮三维接触刚度计算方法，首先采用六面体网格划分，要求齿面上每个正方形网格面积相等；基于有限元法计算光滑齿面接触压力分布，提取接触区节点压强；基于分形理论计算单个正方形网格法向接触刚度，切向接触刚度；计算齿面接触刚度。本发明进一步精确了齿轮接触刚度计算模型，解决了原先计算过程中按照赫兹理论光滑接触的弊端，揭示了粗糙齿面接触过程中粗糙度参数对于齿轮刚度特性的影响规律和计算方法。同时基于有限元的三维接触刚度计算模型，比原先二维计算模型或者计算公式更精确，并且考虑齿轮实际装配和制造过程中误差的影响，从而使齿轮接触刚度更加准确，为齿轮动力学分析奠定了基础。

Description

基于粗糙表面的直齿轮三维接触刚度计算方法

技术领域

本发明属于齿轮力学分析与动力学研究领域，具体涉及一种基于粗糙表面的直齿轮三维接触刚度计算方法。

背景技术

齿轮副内部激励是在齿轮副啮合过程中在系统内部产生的，也是齿轮系统动力学研究中的重点和难点。齿轮系统的内部激励主要分为三种：时变刚度激励、时变误差激励和啮合冲击激励。对于时变刚度激励主要由三部分构成：轮齿的弯曲刚度、剪切刚度和接触刚度，目前应用较为广泛的是IOS6336啮合刚度计算公式以及将齿轮的刚度激励等效为傅里叶展式。接触刚度的计算采用的是基于光滑表面的赫兹接触理论进行推导的，因此不能精确模拟实际轮齿粗糙表面接触状态的变形。所以为了建立更有效的齿轮动力学模型，有必要深入研究基于粗糙表面的直齿轮三维接触刚度的计算方法。

发明内容

本发明的目的是提供一种基于粗糙表面的直齿轮三维接触刚度的计算方法，该方法既能有效的精确齿轮接触刚度模型的准确性，同时又充分考虑到由于装配误差、制造误差等因素下的轮齿啮合载荷分布不均的问题。

为了实现上述目的，本发明采用的技术方案是：

S1采用六面体网格划分，要求齿面上每个正方形网格面积A₀相等；

S2基于有限元法计算光滑齿面接触压力分布，提取接触区节点压强P_i；

S3基于分形理论计算单个正方形网格法向接触刚度K_N，切向接触刚度K_T；

S4计算齿面接触刚度K。

所述S1具体包括以下步骤：

齿轮对采用六面体网格划分，约束六面体网格的边长为lmm，因此得到齿轮 对的接触面上均是标准的正方形网格，网格面积A₀＝l²mm²。

所述步骤2)具体包括以下步骤：

将网格划分好的齿轮对导入有限元仿真软件中，设置静力学分析环境，通过定义每个齿轮的材料属性、接触对、固定被动齿轮的中心孔，然后在主动齿轮轴线上施加扭矩T，以载荷步的方式得到当前齿轮特定接触面的应力和应变云图，此时提取特定接触面上接触区节点压强。

所述步骤3)具体包括以下步骤：

在有限元分析中使用的网格已经足够小，将提取到的接触区节点压强P_i看作是齿面上正方形网格内的接触压强，因此齿面上不同的正方形网格的接触压力表示为：

F_i＝P_iA₀ (1)

由于机械加工的齿面具有分形特征，存在局部与整体自相似的特性。粗糙齿面的轮廓由W-M函数表示为：

其中，n为频率指数，是最低及最高截止频率对应的序列，表示粗糙表面轮廓高度，D表示分形维数，表示谱密度的尺寸参数，G表示粗糙度参数，x表示采样长度坐标。采用相同的材料和加工方法得到齿轮样品方块，在三维形貌测量仪上获得齿面形貌点坐标，再采用功率谱密度方法拟合W-M函数得到D和G两个参数。

在齿轮对的接触模型中两个粗糙表面的接触被简化为一个刚性表面与粗糙表面的接触。虚线是中性面表示理想光滑齿轮对的接触位置，将中性面的上方的粗糙表面凸出的部分定义为微凸体。粗糙表面是由无数多个大小不一的圆形微凸体组成，假设每个微凸体之间彼此分立，相互作用忽略不计。当粗糙表面与刚性表面相互接触时，不同微凸体在压力的作用下发生弹性或塑性形变，则微凸体横截面积a′的分布规律满足：

式中，a′为微凸体变性前与刚性表面相交的横截面积，a表示微凸体变形后与刚性面的接触面积定义为真实接触面积，当微凸体发生弹性变形时，a′＝2a,当发生塑性形变时，a′＝a，D表示分形维数，a′_l表示最大横街面积，表示域拓展因子，可由超越方程(4)求得：

[ψ^(2-D)/2-(1+ψ^-D/2)^-(2-D)/D]/[(2-D)/D]＝1 (4)

根据微凸体横截面积的不同，将微凸体的变形分为弹性和塑性变形，则决定微凸体发生弹性或塑性变形的临界横截面积是由齿轮对的材料属性决定的：

H表示较软材料的硬度，，Y表示屈服强度值；表示等效弹性模量,,E₁、E₂、v₁、v₂分别表示两接触齿面的弹性模量和泊松比。如果微凸体的横截面积a′>a′_c，这发生弹性变形，若横截面积a′≤a′_c，这发生塑性变形。粗糙表面单个微凸体法向变形量由W-M函数中峰-谷间幅值来表示，

δ＝2^3-DG^(D-1)(lnγ)^1/2π^(D-2)/2(a′)^(2-D)/2 (6)

单个微凸体曲率R为，

对于单个微凸体的弹性或塑性形变，其法向载荷f与横截面积a′满足如下关系，

f_p＝Ha′ (9)

式中，下角标e和p分别代表弹性和塑性变形，表示等效弹性模量，H表示较软材料的硬度，，R和δ分别表示单个微凸体的曲率和法向变形量。

单个正方体网格总的弹性变形力F_E和总塑性变形力F_P表示为：

式中符号的意义与上文相同，由此得到单个正方形网格内总的接触压力,

F_i＝F_E+F_P (12)

所述步骤3)还包括以下步骤：

根据式(1)和(12)可以得到单个正方形网格当中最大横截面积a′_l。

对于单个微凸体，其变形分为完全弹性变形及完全塑性变形两个个阶段，其中塑性变形阶段刚度为0。根据刚度定义，发生完全弹性变形的单个微凸体法向接触刚度可表示为，

因此通过在完全弹性区域进行积分，得到单个正方形网格内法向接触刚度

研究表明单个微凸体法切向变形量δ_t表示为

式中μ表示静摩擦系数，G′表示结合面的等效剪切模量且满足 1/G′＝(2-ν₁)/G₁+(2-ν₂)/G₂，其中ν₁、ν₂、G₁、G₂分别表示两材料的泊松比和剪切模量，t和f分别表示单个微凸体的切向和法向载荷，r表示单个微凸体的真实接触区域半径且a表示真实接触面积。

则单个微凸体的切向刚度可以表示为，

对于单个微凸体切向载荷t与法向载荷f满足如下关系，

T＝τ_bA_r，τ_b表示较软材料的剪切强度，A_r表示正方形网格的实际接触面积

则单个正方形网格的切向接触刚度K_T可以表示为，

所述步骤4)具体包括以下步骤：

设特定齿面上有m个正方形网格，则接触齿面的法向接触刚度K₁表示为，

接触齿面的切向接触刚度K₂为

则齿面的接触刚度可以表示为

K＝K₁cosα+K₂sinα (22)

式中α表示齿轮的压力角。

相对于现有技术，本发明的有益效果为：

本发明提出的基于粗糙表面的直齿轮三维接触刚度计算方法，进一步精确了齿轮接触刚度计算模型，解决了原先计算过程中按照赫兹理论光滑接触的弊端，揭示了粗糙齿面接触过程中粗糙度参数对于齿轮刚度特性的影响规律和计算方法。同时本发明是基于有限元的三维接触刚度计算模型，比原先的二维计算模型或者计算公式更精确，并且考虑到齿轮实际装配和制造过程中误差的影响，从而使齿轮接触刚度更加准确，为齿轮动力学分析奠定了基础。

建立的直齿轮接触刚度模型中首次包含了粗糙度有关参数，间接阐述了加工方式与齿轮接触刚度的关系。

通过特殊有限元网格划分方法，首次将分形理论分别应用于接触区不同正方向网格中，通过对不同接触区内网格的累加，得到了齿轮三维接触刚度，避免了因齿轮对因装配或制造误差引起的偏载对啮合刚度的影响。

采用的技术方案中S1、S2、S3对除了直齿轮以外的其他类型齿轮同样适用。

附图说明

图1是齿轮对简化接触模型图。

图2是单个微凸体接触区域图。

图3是直齿轮接触刚度转化示意图。

具体实施方式

下面结合附图1-3和实施例对本发明作进一步说明。

一种基于粗糙表面的直齿轮三维接触刚度计算方法，包括以下步骤：

1)齿轮对特定网格划分，要求各网格面积A₀相等；

2)基于有限元法计算光滑齿面接触压力分布，提取接触区节点压强P_i；

3)基于分形理论计算单个正方形网格法向接触刚度K_N，切向接触刚度K_T；

4)计算齿面接触刚度K。

所述步骤1)具体包括以下步骤：

齿轮对采用六面体网格划分，约束六面体网格的边长为2mm，因此可以得 到齿轮对的接触面上均是标准的正方形网格，网格面积A₀＝4mm²。

所述步骤2)具体包括以下步骤：

将网格划分好的齿轮对导入有限元仿真软件中，设置静力学分析环境，通过定义每个齿轮的材料属性，接触对，固定被动齿轮的中心孔，然后在主动齿轮轴线上施加扭矩T，以载荷步的方式可以得到当前齿轮特定接触面的应力和应变云图，此时可以提取特定接触面上接触区节点压强。

所述步骤3)具体包括以下步骤：

在有限元分析中使用的网格已经足够小，可以将提取到的接触区节点压强P_i看作是齿面上正方形网格内的接触压强，因此齿面上不同的正方形网格的接触压力可以表示为：

F_i＝P_iA₀ (1)

由于机械加工的齿面具有分形特征，存在局部与整体自相似的特性。粗糙齿面的轮廓可以由W-M函数表示为：

其中n为频率指数，是最低及最高截止频率对应的序列，表示粗糙表面轮廓高度，D表示分形维数，表示谱密度的尺寸参数，G表示粗糙度参数，x表示采样长度坐标。采用相同的材料和加工方法得到齿轮样品方块，在三维形貌测量仪上获得齿面形貌点坐标，再采用功率谱密度方法拟合W-M函数得到D和G两个参数。

在齿轮对的接触模型中两个粗糙表面的接触被简化为一个刚性表面与粗糙表面的接触，如图1所示。虚线是中性面表示理想光滑齿轮对的接触位置，将中性面的上方的粗糙表面凸出的部分定义为微凸体，如图2所示。粗糙表面是由无数多个大小不一的圆形微凸体组成，假设每个微凸体之间彼此分立，相互作用忽略不计。当粗糙表面与刚性表面相互接触时，不同微凸体在压力的作用下发生弹性或塑性形变，则微凸体横截面积a′的分布规律满足：

式中，a′为微凸体变性前与刚性表面相交的横截面积，a表示微凸体变形后与刚性面的接触面积定义为真实接触面积，当微凸体发生弹性变形时，a′＝2a,当发生塑性形变时，a′＝a，D表示分形维数，a′_l表示最大横街面积，表示域拓展因子，可由超越方程(4)求得：

[ψ^(2-D)/2-(1+ψ^-D/2)^-(2-D)/D]/[(2-D)/D]＝1 (4)

根据微凸体横截面积的不同，将微凸体的变形分为弹性和塑性变形，则决定微凸体发生弹性或塑性变形的临界横截面积是由齿轮对的材料属性决定的：

H表示较软材料的硬度，，Y表示屈服强度值；表示等效弹性模量,,E₁E₂v₁v₂分别表示两接触齿面的弹性模量和泊松比。如果微凸体的横截面积a′>a′_c，这发生弹性变形，若横截面积a′≤a′_c，这发生塑性变形。粗糙表面单个微凸体法向变形量可由W-M函数中峰-谷间幅值来表示，

δ＝2^3-DG^(D-1)(lnγ)^1/2π^(D-2)/2(a′)^(2-D)/2 (6)

单个微凸体曲率R为，

对于单个微凸体的弹性或塑性形变，其法向载荷f与横截面积a′满足如下关系，

f_p＝Ha′ (9)

式中，下角标e和p分别代表弹性和塑性变形，表示等效弹性模量，E₁E₂v₁v₂分别表示两接触齿面的弹性模量和泊松比，H表示较软材料的硬度，，R和δ分别表示单个微凸体的曲率和法向变形量。

单个正方体网格总的弹性变形力F_E和总塑性变形力F_P可以表示为：

式中符号的意义与上文相同，由此得到单个正方形网格内总的接触压力,

F_i＝F_E+F_P (12)

所述步骤3)具体包括以下步骤：

根据式(1)和(12)可以得到单个正方形网格当中最大横截面积a′_l。

对于单个微凸体，其变形分为完全弹性变形及完全塑性变形两个个阶段，其中塑性变形阶段刚度为0。根据刚度定义，发生完全弹性变形的单个微凸体法向接触刚度可表示为，

因此通过在完全弹性区域进行积分，得到单个正方形网格内法向接触刚度

研究表明单个微凸体法切向变形量δ_t可表示为

式中μ表示静摩擦系数，G′表示结合面的等效剪切模量且满足1/G′＝(2-ν₁)/G₁+(2-ν₂)/G₂，其中ν₁ ν₂ G₁ G₂分别表示两材料的泊松比和剪切模量，t和f分别表示单个微凸体的切向和法向载荷，r表示单个微凸体的真实接触区域半径且式中a表示真实接触面积。

则单个微凸体的切向刚度可以表示为，

对于单个微凸体切向载荷t与法向载荷f满足如下关系，

其中T＝τ_bA_r，τ_b表示较软材料的剪切强度，A_r表示正方形网格的实际接触面积

则单个正方形网格的切向接触刚度K_T可以表示为，

所述步骤4)具体包括以下步骤：

设特定齿面上有m个正方形网格，则接触齿面的法向接触刚度K₁表示为，

接触齿面的切向接触刚度K₂为

根据图3可得齿面的接触刚度可以表示为

K＝K₁ cosα+K₂ sinα (22)

式中α表示齿轮的压力角。

Claims (6)

1.基于粗糙表面的直齿轮三维接触刚度计算方法，其特征在于：本方法包括如下步骤，

S1采用六面体网格划分，要求齿面上每个正方形网格面积A₀相等；

S2基于有限元法计算光滑齿面接触压力分布，提取接触区节点压强P_i；

S3基于分形理论计算单个正方形网格法向接触刚度K_N，切向接触刚度K_T；

S4计算齿面接触刚度K。

2.根据权利要求1所述的基于粗糙表面的直齿轮三维接触刚度计算方法，其特征在于：

所述S1具体包括以下步骤：

齿轮对采用六面体网格划分，约束六面体网格的边长为lmm，因此得到齿轮对的接触面上均是标准的正方形网格，网格面积A₀＝l²mm²。

3.根据权利要求1所述的基于粗糙表面的直齿轮三维接触刚度计算方法，其特征在于：所述步骤2)具体包括以下步骤：

将网格划分好的齿轮对导入有限元仿真软件中，设置静力学分析环境，通过定义每个齿轮的材料属性、接触对、固定被动齿轮的中心孔，然后在主动齿轮轴线上施加扭矩T，以载荷步的方式得到当前齿轮特定接触面的应力和应变云图，此时提取特定接触面上接触区节点压强P_i。

4.根据权利要求1所述的基于粗糙表面的直齿轮三维接触刚度计算方法，其特征在于：所述步骤3)具体包括以下步骤：

在有限元分析中使用的网格已经足够小，将提取到的接触区节点压强P_i看作是齿面上正方形网格内的接触压强，因此齿面上不同的正方形网格的接触压力表示为：

F_i＝P_iA₀ (1)

由于机械加工的齿面具有分形特征，存在局部与整体自相似的特性；粗糙齿面的轮廓由W-M函数表示为：

$\begin{array}{cc}z\left(x\right)=G^{(D-1)}\underset{n=0}{\overset{n_{max}}{\&Sigma;}}\left(\frac{cos2{\mathrm{\&pi;\&gamma;}}^{n}x}{{\mathrm{\&gamma;}}^{(2-D)n}}\right)& 1<D<2,\mathrm{\&gamma;}>1\end{array}---\left(2\right)$

其中，n为频率指数，n＝0和n_max是最低及最高截止频率对应的序列，z表示粗糙表面轮廓高度，D表示分形维数，γ表示谱密度的尺寸参数，G表示粗糙度参数，x表示采样长度坐标；采用相同的材料和加工方法得到齿轮样品方块，在三维形貌测量仪上获得齿面形貌点坐标，再采用功率谱密度方法拟合W-M函数得到D和G两个参数；

在齿轮对的接触模型中两个粗糙表面的接触被简化为一个刚性表面与粗糙表面的接触；虚线是中性面表示理想光滑齿轮对的接触位置，将中性面的上方的粗糙表面凸出的部分定义为微凸体；粗糙表面是由无数多个大小不一的圆形微凸体组成，假设每个微凸体之间彼此分立，相互作用忽略不计；当粗糙表面与刚性表面相互接触时，不同微凸体在压力的作用下发生弹性或塑性形变，则微凸体横截面积a′的分布规律满足：

$n\left(a^{\&prime;}\right)=\frac{1}{2}{\mathrm{D\&psi;}}^{(2-D)/2}{a}_L^{\&prime;D/2}{a}^{\&prime;-(2+D)/2}0<a^{\&prime;}\&le;a_L^{\&prime;}---\left(3\right)$

式中，a′为微凸体变性前与刚性表面相交的横截面积，a表示微凸体变形后与刚性面的接触面积定义为真实接触面积，当微凸体发生弹性变形时，a′＝2a,当发生塑性形变时，a′＝a，D表示分形维数，a′_l表示最大横街面积，表示域拓展因子，可由超越方程(4)求得：

[ψ^(2-D)/2-(1+ψ^-D/2)^-(2-D)/D]/[(2-D)/D]＝1 (4)

根据微凸体横截面积的不同，将微凸体的变形分为弹性和塑性变形，则决定微凸体发生弹性或塑性变形的临界横截面积是由齿轮对的材料属性决定的：

$a_c^{\&prime;}=2a_c=2{\left(\frac{2E^{*}}{H}\right)}^{2/(D-1)}{G}^2---\left(5\right)$

H表示较软材料的硬度，H＝2.8Y，Y表示屈服强度值；E^*表示等效弹性模量,E₁、E₂、v₁、v₂分别表示两接触齿面的弹性模量和泊松比；如果微凸体的横截面积a′>a′_c，这发生弹性变形，若横截面积a′≤a′_c，这发生塑性变形；粗糙表面单个微凸体法向变形量由W-M函数中峰-谷间幅值来表示，

δ＝2^3-DG^(D-1)(lnγ)^1/2π^(D-2)/2(a′)^(2-D)/2 (6)

单个微凸体曲率R为，

$R=\frac{a^{\&prime;D/2}}{2^{4-D}{\mathrm{\&pi;}}^{D/2}{G}^{D-1}{\left(ln\mathrm{\&gamma;}\right)}^{1/2}}---\left(7\right)$

对于单个微凸体的弹性或塑性形变，其法向载荷f与横截面积a′满足如下关系，

$f_e=\frac{4}{3}{E}^{*}{R}^{\frac{1}{2}}{\mathrm{\&delta;}}^{\frac{3}{2}}---\left(8\right)$

f_p＝Ha′ (9)

式中，下角标e和p分别代表弹性和塑性变形，表示等效弹性模量，H表示较软材料的硬度，H＝2.8Y，R和δ分别表示单个微凸体的曲率和法向变形量；

单个正方体网格总的弹性变形力F_E和总塑性变形力F_P表示为：

$F_E={\&Integral;}_{a_c^{\&prime;}}^{a_L^{\&prime;}}{f}_{e}n\left(a^{\&prime;}\right){\mathrm{da}}^{\&prime;}\left\{\begin{array}{cc}\frac{2^{(9-2D)/2}D}{3(3-2D){\mathrm{\&pi;}}^{(3-D)/2}}{E}^{*}{\left(\ln \mathrm{\&gamma;}\right)}^{1/2}{G}^{(D-1)}\&times;{\mathrm{\&psi;}}^{1-\frac{D}{2}}{a}_L^{\&prime;\frac{D}{2}}(a_L^{\&prime;\frac{3}{2}-D}-a_c^{\&prime;\frac{3}{2}-D})& (D\&NotEqual;1.5)\\ 2E^{*}{\mathrm{\&psi;}}^{\frac{1}{4}}{\mathrm{\&pi;}}^{-\frac{3}{4}}{\left(\ln \mathrm{\&gamma;}\right)}^{\frac{1}{2}}{G}^{\frac{1}{2}}{a}_L^{\&prime;\frac{3}{4}}\ln \frac{a_L^{\&prime;}}{a_c^{\&prime;}}& (D=1.5)\end{array}\right.---\left(10\right)$

$F_P={\&Integral;}_0^{a_c^{\&prime;}}{f}_{p}n\left(a^{\&prime;}\right){\mathrm{da}}^{\&prime;}=\frac{{\mathrm{HD\&psi;}}^{(2-D)/2}}{2-D}{a_L^{\&prime;}}^{D/2}{a}_c^{\&prime;(2-D)/2}---\left(11\right)$

由此得到单个正方形网格内总的接触压力,

F_i＝F_E+F_P (12)。

5.根据权利要求4所述的基于粗糙表面的直齿轮三维接触刚度计算方法，其特征在于：

所述步骤3)还包括以下步骤：

根据式(1)和(12)可以得到单个正方形网格当中最大横截面积a′_l；

对于单个微凸体，其变形分为完全弹性变形及完全塑性变形两个个阶段，其中塑性变形阶段刚度为0；根据刚度定义，发生完全弹性变形的单个微凸体法向接触刚度可表示为，

$k_n=\frac{{\mathrm{df}}_e}{d\mathrm{\&delta;}}=\frac{{\mathrm{df}}_e/{\mathrm{da}}^{\&prime;}}{d\mathrm{\&delta;}/{\mathrm{da}}^{\&prime;}}=\frac{4E^{*}(3-D)}{3\sqrt{2\mathrm{\&pi;}}(2-D)}{a}^{\&prime;\frac{1}{2}}---\left(13\right)$

因此通过在完全弹性区域进行积分，得到单个正方形网格内法向接触刚度

$K_N={\&Integral;}_{a_c^{\&prime;}}^{a_L^{\&prime;}}{k}_{n}n\left(a^{\&prime;}\right){\mathrm{da}}^{\&prime;}=\frac{4E^{*}(3-D)D}{3\sqrt{2\mathrm{\&pi;}}(1-D)(2-D)}\&times;{\mathrm{\&psi;}}^{1-D/2}{a}_l^{\&prime;D/2}(a_l^{\&prime;(1-D)/2}-a_c^{\&prime;(1-D)/2})---\left(14\right)$

研究表明单个微凸体法切向变形量δ_t表示为

${\mathrm{\&delta;}}_t=\frac{3\mathrm{\&mu;}f}{16G^{\&prime;}r}\&lsqb;1-{(1-\frac{t}{\mathrm{\&mu;}f})}^{2/3}\&rsqb;---\left(15\right)$

式中μ表示静摩擦系数，G′表示结合面的等效剪切模量且满足1/G′＝(2-ν₁)/G₁+(2-ν₂)/G₂，其中ν₁、ν₂、G₁、G₂分别表示两材料的泊松比和剪切模量，t和f分别表示单个微凸体的切向和法向载荷，r表示单个微凸体的真实接触区域半径且a表示真实接触面积；

则单个微凸体的切向刚度可以表示为，

$k_t=\frac{dt}{{\mathrm{d\&delta;}}_t}---\left(16\right)$

对于单个微凸体切向载荷t与法向载荷f满足如下关系，

$\frac{t}{f}=\frac{T}{F_i}---\left(17\right)$

T＝τ_bA_r，τ_b表示较软材料的剪切强度，A_r表示正方形网格的实际接触面积

$A_r={\&Integral;}_0^{a_c^{\&prime;}}{a}^{\&prime;}n\left(a^{\&prime;}\right){\mathrm{da}}^{\&prime;}+\frac{1}{2}{\&Integral;}_{a_c^{\&prime;}}^{a_L^{\&prime;}}{a}^{\&prime;}n\left(a^{\&prime;}\right){\mathrm{da}}^{\&prime;}=\frac{D}{4-2D}{\mathrm{\&psi;}}^{1-\frac{D}{2}}{a_L^{\&prime;}}^{D/2}(a_L^{\&prime;1-D/2}+a_c^{\&prime;1-D/2})---\left(18\right)$

则单个正方形网格的切向接触刚度K_T表示为，

$K_T={\&Integral;}_{{a_c}^{\&prime;}}^{{a_L}^{\&prime;}}{k}_{t}n\left(a^{\&prime;}\right){\mathrm{da}}^{\&prime;}=\frac{8G^{\&prime;}{\mathrm{D\&psi;}}^{(2-D)/2}}{\sqrt{2\mathrm{\&pi;}}(1-D)}{(1-\frac{T}{\mathrm{\&mu;}P})}^{1/3}\&times;{a_L}^{\&prime;D/2}({a_L}^{\&prime;(1-D)/2}-{a_c}^{\&prime;(1-D)/2})---\left(19\right).$

6.根据权利要求1所述的基于粗糙表面的直齿轮三维接触刚度计算方法，其特征在于：

所述步骤4)具体包括以下步骤：

设特定齿面上有m个正方形网格，则接触齿面的法向接触刚度K₁表示为，

$K_1=\underset{i=1}{\overset{m}{\&Sigma;}}{K}_N---\left(20\right)$

接触齿面的切向接触刚度K₂为

$K_2=\underset{i=1}{\overset{m}{\&Sigma;}}{K}_T---\left(21\right)$

则齿面的接触刚度可以表示为

K＝K₁cosα+K₂sinα (22)

式中α表示齿轮的压力角。