齿轮有限元接触分析是典型的状态非线性问题，求解过程为准确捕捉齿面接触情况及应力分布状态，就必须增大有限元模型的网格密度．但这势必会导致模型的规模维度增大，计算时间延长．故传统的网格划分方式不能满足齿轮有限元接触分析的准确性与高效性要求，因此许多学者探索了齿面啮合接触区网格局部加密的有限元精确建模方法．有限元模型局部加密的关键是对网格密度不同的边界进行协调过渡处理，通常有四种方式来实现协调过渡．第一种是通过六面体单元退化为四面体单元来实现过渡[1-2]，ANSYS等通用商业软件采用这种方式，虽具备一定自适应性，但模型节点和单元规模极为庞大计算效率低．第二种是通过高低阶单元连接实现过渡，文献[3-4]根据这种方式来加密渐开线齿面网格．齿面最外层每四个8节点六面体单元与轮齿内部一个20节点六面体单元过渡连接实现齿面网格加密．但由于对整个齿面都进行了加密，在计算规模上缩减空间不大．第三种是通过建立多点约束方程方式实现过渡[5-7]，文献[6]采用子结构法将齿面加密区与非加密区装配成整体，并利用形函数插值的方式建立多点约束方程；文献[7]在Abaqus中通过分界面建立绑定接触来建立多点约束方程，这种方式提高了求解效率，但不能保证网格节点的严格连通性，在连接处会产生较大的应力梯度，无法保证求解精度．第四种是通过非规则形状的网格剖分实现过渡[8-9]，文献[8]提出一种基于分级剖分法的有限元精确齿面建模方法，对齿面接触带进行非规则六面体网格剖分加密来提高求解效率，但在齿顶和齿根关键部位仍需多点约束耦合，因此该算法还欠缺一定自适应性．目前在众多的网格划分方式中，栅格法[10-11]具有高度的自动化，易实现网格密度控制，且具有较强的自适应性，广泛用于复杂二维和三维模型的网格划分．因此为解决上述四种加密模式存在的局限性问题，本研究结合栅格法的特点，提出将栅格法-27分法加密模板扩展应用到齿面局部加密中，建立了齿面有限元局部加密算法，给出了齿面局部加密的建模流程，形成了一种基于栅格法的人字齿轮有限元接触分析的精确建模方法．1 人字齿面有限元局部加密算法1.1　确定齿面接触区域和栅格法加密模板人字齿面由两个左右旋向的渐开线螺旋面组成，第i个齿的双螺旋齿廓面表达[12]为r=rbe(θi+λ)-rbθe1(θi+λ)+pλk;θi=θ+θi0, (1)式中：θi0=π/2-sb/(2rb)+(i-1)×2π/Z，其中，rb为基圆半径，sb为基圆齿厚；θ=α+inv α，α为齿面某一点的压力角，其中inv α为渐开线函数；Z为齿轮齿数；p=Ph/(2π)，其中Ph为螺旋线导程；(θi,λ)为参变量，其中λ为螺旋面的旋向符号．齿面接触线可由齿轮啮合面与齿面的交线确定，啮合面由啮合角α'和基圆半径rb1具体[8]为-xtan α'+y-rb1/cosα'=0．(2)根据齿轮参数和齿面方程(1)进行参数化有限元建模[13]，每个瞬时的接触线由式(1)和式(2)联立求解确定．根据赫兹接触理论，将啮合面沿其法线上下平移两倍的接触半宽b，两个平行平面截取的齿面范围确定为潜在接触区．如图1所示为右旋螺旋齿廓的有限元网格模型，将潜在接触区(两条红线)映射到有限元网格中，可以看出其包含在一个阶梯状分布的六面体网格带中．10.13245/j.hust.220316.F001图1齿面潜在接触区域示意图对包含潜在接触区的阶梯状六面体网格带进行加密，一方面要增大网格密度使网格形状逼近渐开线螺旋齿面；另一方面要实现加密网格带和非加密区的协调过渡，使网格节点一一对应，不存在孤立节点．假设提取的初始数据中包含n个节点和m个单元，节点坐标的形式为{(xi,yi,zi)|i=1，2，…，n}，单元编号的形式为{kij|i=1,2,…,m;j=1,2,…,8}．为清楚表示齿面单元节点间的拓扑关系，通过层次结构来描述网格单元．将按顺序排列的8节点集合定义为单元元素Vi={ki1,ki2,ki3，ki4,ki5,ki6,ki7,ki8}对应的下标i=1，2，…，m，将m个单元元素Vi定义为单元集合E={Vi|i=1,2,⋯,m}．根据阶梯状网格带的特点，为保证加密后与相邻非加密单元的节点连通性，先对网格带中完全包含在接触区范围内单元选用27分法[11]中的面加密模板进行面加密；再对阶梯网格带的拐角边界选用相邻边加密模板加密；最后对除拐角外的其他边界单元选取边加密模板加密．如图2给出了加密单元与非加密单元间的网格衔接过渡状态．下面将展开说明基于这三种加密模板的齿面局部加密算法．10.13245/j.hust.220316.F002图2基于三种加密模板的齿面局部加密示意图1.2　单元编号索引矩阵及单元协调过渡加密算法为实现自动化齿面加密流程，提出单元协调过渡加密算法，通过单元间的节点连通关系构建编号索引矩阵，计算索引值判别过渡单元的加密模式．首先，定义齿面最外层全部ms个六面体单元为单元集合Es，阶梯网格带中完全包含在接触区范围内的m*个单元为集合E*，并令其全部采用面加密模板；将属于单元集合Es但不属于E*的六面体单元定义为过渡单元集合Eq，表示为：Es={Vii=1,2,⋯,ms};E*={Vii=1,2,⋯,m*};Eq={ViVi∈Es⋂Vi∉E*}.然后建立单元编号索引矩阵Be，假设集合Eq中包含m个单元Vi，则Be为m×8阶的矩阵，具体为Be=[biji=1,2,⋯,m;j=1,2,⋯,8]．以索引值判断过渡集合Eq中单元的加密模式，即确定集合Eq中单元节点与其相邻的集合E*节点的连通关系．定义集合E*所包含的全部节点为K*集合，若集合Eq中单元元素Vi的节点kij属于集合K*，则索引矩阵对应位置bij=1，表示节点相接；若集合Eq中Vi的节点kij不属于集合K*，则在索引矩阵对应位置bij=0，表示节点不相接．因此，将单元编号索引矩阵每行中元素为1的个数定义为单元编号索引值SVi．当SVi=0或SVi=2时，分别表示此单元与集合E*中单元不相接或仅有一条边相接，单元无须加密；SVi=4表示此单元与集合E*中单元有一个面相接，采用边加密模板；SVi=6表示此单元与集合E*中单元有两个面相接，采用相邻边加密模板．以索引值判别集合Eq中单元的加密模式，即单元协调过渡加密算法．以图3为例说明单元协调过渡加密算法．图中取网格带边界相邻的9个单元，假设单元1，2，3，6属于完全包含在潜在接触区域的单元集合E*，采用面加密模板加密，而单元4，5，7，8，9属于过渡单元集合Eq，过渡单元集合Eq的编号索引矩阵Be为Be=1100110011101110000000000100010011001100．10.13245/j.hust.220316.F003图3单元协调过渡加密的示意图根据单元编号索引值来判断单元所属的加密模式：SV5=6则单元5进行相邻边加密，SV4=SV9=4则单元4，9进行边加密，SV7=0，SV8=2则单元7，8不加密，图3给出了应用加密模板进行加密后的效果．1.3　编号轮转算法依据单元协调过渡加密算法判断单元所属的加密模式后，提出编号轮转算法来判别过渡单元的待加密边．根据每行索引矩阵中1元素的位置即可判断待加密边的位置，然后利用节点循环移位的方式重新排列节点编号，统一节点编号后即可应用加密模板计算加密节点坐标．因为索引矩阵每行的前四位与后四位一致，所以仅取每行索引矩阵的前四位定义为判别向量QVi=[bi1,bi2,bi3,bi4]．当QVi=[0,1,1,0]或QVi=[0,1,1,1]时，将单元编号转换为Vi={ki4,ki1,ki2,ki3,ki8,ki5,ki6,ki7}；当QVi=[0,0,1,1]或QVi=[1,0,1,1]时，将编号转换为Vi={ki3,ki4,ki1,ki2,ki7,ki8,ki5,ki6}；当QVi=[1,0,0,1]或QVi=[1,1,0,1]时，则将单元编号转换为Vi={ki2,ki3,ki4,ki1,ki6,ki7,ki8,ki5}．依据判别向量将单元节点重新排列的上述过程即为编号轮转算法．加密后的节点有可能偏离标准渐开线螺旋齿面，利用齿面节点映射算法[8]修正坐标将网格节点映射到齿廓上．修正后的加密节点坐标整理排序并转化为文本格式输出，通过APDL程序读取数据建立节点node，定义合理的单元类型和材料属性，用E命令连接节点node为单元element完成网格重构．2 齿面局部加密的有限元建模流程上述单元协调过渡加密算法和编号轮转算法能实现对齿面潜在接触区单元的局部加密，并且能在ANSYS中进行参数化建模，具体的齿面局部加密建模流程如图4所示，具体步骤如下．10.13245/j.hust.220316.F004图4一侧螺旋齿面局部细化的建模流程图a.根据齿廓共轭啮合定律确定啮合面，利用接触半宽的预估值确定齿面的潜在接触区域．结合齿面参数利用APDL条件语句可快速判断齿面潜在接触区域；编写提取单元编号和节点坐标APDL子程序导出初始数据，并删除待加密齿面层单元．b.利用单元编号索引值判断过渡单元的加密模式，通过编号轮转算法对单元节点重新编号；然后利用空间8节点等参元映射计算加密子节点坐标；最后通过节点映射算法修正加密节点坐标．空间8节点等参单元[14]是描述空间直棱六面体的完备单元，满足描述上述有限元齿面网格单元的完备性和连续性要求．空间8节点等参数单元从全局坐标系映射到参考坐标系ξηζ的过程如图5所示，将任意六面体单元映射变换为边长为2的立方体单元，参考坐标系ξηζ的原点位于它的形心处．10.13245/j.hust.220316.F005图5加密模板和空间8节点等参单元示意图8节点等参单元内任意一点的全局坐标系与局部坐标系的坐标变换关系表达为：x=∑i=1n=8Nixi;y=∑i=1n=8Niyi;z=∑i=1n=8Nizi．(3)n=8表示为8节点等参数单元，其形函数的具体形式为Ni=(1+ξiξ)(1+ηiη)(1+ζiζ)/8(i=1,2,⋯,8), (4)式中ξi，ηi和ζi为节点i在参考坐标系下的局部坐标，对应于8个节点其值分别为-1或1，即：ξ1,2,⋯,8=1,-1,-1,1,1,-1,-1,1;η1,2,⋯,8=-1,-1,1,1,-1,-1,1,1;ζ1,2,⋯,8=-1,-1,-1,-1,1,1,1,1.上述加密模板中每个加密子节点，在其8节点等参元参考坐标系中的相对位置(ξ,η,ζ)是确定的．因此可以将每个加密子节点局部坐标(ξ,η,ζ)代入形函数计算公式(4)求得Ni，然后利用坐标变换式(3)计算出全局坐标系下的节点坐标．c.编写APDL面加密、边加密和相邻边加密网格重构子程序，分别读取加密节点坐标的文本数据进行网格重构，其他非加密区用原单元重构．图6为表1的人字齿有限元模型局部加密后的效果图，可见本文提出的齿面局部加密算法具备较好的自适应性，能够对齿顶、齿根和齿面各部分进行有效加密，同时保证整个模型网格节点的连通性．局部细化后的网格质量通过ANSYS中的网格检查命令CHECK和MCHECK进行网格质量检查，检测所有单元均满足ANSYS软件中所设定的单元形状参数限定．单元的网格质量评估准则采取计算单元的雅可比矩阵[11]来评估有限元的网格质量，计算所有单元的最小雅可比矩阵行列式的值，其值均大于零，满足有限元分析的基本要求，说明划分单元的网格质量满足有限元精确接触分析的要求．10.13245/j.hust.220316.F006图6有限元人字齿面局部加密的效果图3 外啮合人字齿有限元接触分析以一对外啮合人字齿轮为例，应用所提出的有限元齿面局部加密方法进行精确建模，其人字齿轮的基本参数如下：主动轮齿数为43，从动轮齿数为42，两齿轮法面模数为3.5 mm，法面压力角为22.5°，主动轮螺旋角为26.969°，从动轮螺旋角为-26.969°，法面齿顶高系数为1，法面顶隙系数为0.25，两齿轮变位系数都为0，主从动轮齿宽为147 mm，退刀槽宽度为33 mm，齿轮弹性模量设置为206 GPa，泊松比设置为0.3．在接触齿面间建立接触对；在主动轮轴线齿宽中点处建立mass21单元，与主动轮轮体内表面节点耦合，保留其轴向转动自由度并施加扭矩T=51.18  N⋅m；在从动轮齿宽中点处建立mass21单元与从动轮轮体内表面节点耦合施加全自由度约束；不计摩擦，法向惩罚刚度系数取１，穿透系数取0.1，采用增广拉格朗日法进行求解．在ANSYS软件中进行接触分析，并进行后处理结果提取．图7为人字齿轮双螺旋齿面的接触应力云图，其中齿面的最大接触应力在左旋齿面上为1 501.49 MPa，而右旋齿面的最大接触应力为 1 501.32 MPa，两者相差0.01%，轮齿的左右旋齿面近似平均分配了施加在主动轮的工况载荷．10.13245/j.hust.220316.F007图7人字齿轮的有限元接触应力云图(色标单位：MPa)细化建模和未细化建模有限元求解的最大接触应力对比结果为：未细化的人字齿模型有限元求解的齿面最大接触应力为3 203.1 MPa，上述细化后模型有限元求解的最大接触应力为1 501.49 MPa．可与看出：未细化模型无法精确捕捉齿顶与齿根接触部分的应力梯度变化，求解的最大接触应力远远超出预估值．而网格局部细化的模型准确地逼近了齿廓曲面，将应力梯度变化较好地捕捉，通过对比说明了建模方法的精确性．分别将ANSYS求解的齿面节点应力和节点坐标数据进行提取，分析多对啮合齿面上的接触应力分布状态．图8为人字齿轮齿宽方向B上的最大接触应力σmax的分布曲线，可以看出：左旋齿面和右旋齿面的应力沿齿宽方向呈对称分布，并且在齿面啮入和啮出部位有明显的应力集中．图9为人字齿轮接触应力σ投影到啮合平面上(齿轮两基圆的内公切线方向即齿廓啮合线方向D和齿宽方向B)的分布状态，图中曲面忽略其应力集中部分，齿面的接触应力沿接触线呈均匀分布，并且应力曲面每一截面呈现赫兹接触理论的抛物线形状，忽略其应力集中部分后平均接触应力为725.75 MPa，与赫兹接触理论的计算值744.34 MPa相对误差为2.56%．10.13245/j.hust.220316.F008图8人字齿轮齿宽方向的接触应力分布曲线10.13245/j.hust.220316.F009图9人字齿轮啮合平面上投影的应力分布状态齿轮的啮入啮出部位出现了应力集中现象，因此在局部细化算法中通过坐标修正对主从动齿轮进行了齿顶修形，修形高度均为0.55 mm，最大修形量均为0.01 mm，修形曲线均选择Walker曲线[15]．齿顶修形有效避免了齿顶和齿根啮合过程中的边缘接触情况，从图10的修形后的齿宽方向上的接触应力曲线和图11的修形后啮合平面上投影的接触应力分布状态可以看出：最大接触应力从1 501.49 MPa下降到873.38 MPa，且每个啮入啮出部位都有明显改善．10.13245/j.hust.220316.F010图10齿廓修形后齿宽方向的接触应力分布曲线10.13245/j.hust.220316.F011图11齿廓修形后啮合平面上投影的应力分布状态人字齿轮的建模与分析实例显示了本加密方法能保证节点连通性的同时实现对齿根、齿顶和齿面各部位的自适应加密，有效协调了接触分析求解精度和效率之间的矛盾，为齿轮的加强接触分析提供有效的建模方法．