一种基于流固耦合的直齿轮弹流润滑仿真分析方法

摘要

本发明涉及一种基于流固耦合的直齿轮弹流润滑仿真分析方法，属于计算机辅助设计技术领域，采用流固耦合的方法针对直齿轮齿面不同啮合位置精准建模并进行弹流润滑分析，所述方法包括以下步骤：步骤1：针对直齿轮啮合位置进行模型等效；步骤2：针对直齿轮齿面不同啮合位置精准建模；步骤3：绘制模型有限元网格；步骤4：将真实齿轮齿面粗糙度写入模型；步骤5：根据齿轮真实工况设置仿真模型边界条件；步骤6：模型后处理；步骤7：验证与归纳。本发明创造性地采用流固耦合的方法针对直齿轮齿面不同啮合位置进行弹流润滑分析，更加真实地反映齿轮齿面不同啮合位置润滑的实际情况，为提供高质量的喷油润滑设计方案提供帮助。

说明书

技术领域

本发明涉及计算机辅助设计技术领域，具体涉及一种基于流固耦合的直齿轮弹流润滑仿真分析方法。

背景技术

准确掌握弹流润滑情况，可以显著提升齿轮、轴承等零件的使用寿命。齿轮弹流润滑的分析途径主要有两种：有限元方法和理论计算方法。有限元方法大多忽略结构模型的变形，仅采用流体力学方法对流体区域进行计算，不能得到经典的弹流润滑形态，如压力二次高峰、油膜颈缩等现象。部分有限元方法即使考虑了耦合系统对弹流的影响，也仅限于理想弹流润滑分析，忽略了实际弹流过程的热效应、流变效应和粗糙形貌效应等因素对弹流润滑的影响，所以得到的油膜厚度解及摩擦系数解不够准确，且研究方法均没有考虑齿轮表面油液量的影响，但在实际工作中，啮合轮齿之间油液量大小会直接影响齿面的摩擦系数，进而影响齿轮的啮合生热量。不同箱体之间由于形状、内部轴系布置及润滑系统的不同，箱体内各个位置的油液分布差异巨大。因此，在考虑齿轮表面油液量对于不同啮合区摩擦系数影响的基础上，对齿轮箱内部齿轮啮合面喷射参数优化十分重要。

而理论计算方法是基于Reynolds方程的数值解法来求解，仅根据一类弹流接触问题进行建模，如点接触或线接触，而且对计算区域进行了理想化假设，脱离了工程实际，对于其它复杂的弹流接触问题也难得到准确的解。

因此，开发一种不仅能考虑齿轮服役工况下不同因素影响，还能运用于复杂的接触问题来获取不同啮合区摩擦系数，以此优化齿轮箱内部齿轮啮合面喷射参数的弹流润滑仿真方法显得格外重要。

发明内容

有鉴于此，本发明的目的在于提供一种基于热流固耦合的直齿轮弹流润滑仿真分析方法，通过对MATLAB粗糙度数据的写入以及考虑真实情况下直齿轮的载荷、滑动速度等边界条件，考虑了粗糙度效应和滑滚比，从而有效地模拟了齿轮服役工况下不同啮合位置的弹流润滑情况，并通过线接触弹流润滑油膜厚度和摩擦系数探索齿轮箱内部滑油喷射参数的最佳位置和喷射量，解决了现有技术无法得到准确解、适应性差、考虑因素少的问题。

本发明的目的通过以下技术方案实现：

一种基于流固耦合的直齿轮弹流润滑仿真分析方法，采用流固耦合的方法对齿轮齿面不同啮合位置进行解析分析获取边界条件，通过弹流润滑分析获取整个齿面润滑数据，所述方法包括以下步骤：

步骤1：等效

将直齿轮啮合等效为两当量弹性圆柱体通过油膜接触；

步骤2：建模

根据等效尺寸建立不同啮合位置的微米级弹流模型，所述微米级弹流模型包括等效结构模型和齿轮啮合时的滑油流体模型；

步骤3：绘制模型网格；

步骤4：写入粗糙度；

通过实验得到齿面粗糙度数据，将粗糙度数据以坐标形式一一对应写入流固耦合面的节点上；

步骤5：设置边界条件

采用瞬态结构有限元技术，对包含轴承箱体的齿轮啮合进行动力学仿真分析，提取齿轮服役工况下其啮合界面的压力分布有限元仿真分析结果，将有限元仿真结果写入弹流模型；

步骤6：后处理

对弹流模型进行流固耦合仿真计算，用MATLAB读取计算结果，绘制三维油膜压力和三维油膜厚度图；

步骤7：验证与归纳

分析各工况的油膜厚度和不同啮合位置的油膜、油压、摩擦系数，与理论计算结果进行对比分析。

进一步，步骤3中，两个当量弹性圆柱面绘制三层密网格，其余区域为一层网格。

进一步，步骤2中，根据赫兹接触的接触半宽建立有限元模型。

进一步，步骤5中，所述有限元仿真结果包括润滑油油温、约束、加载力和壁面速度。

进一步，步骤5中，还包括通过解析法计算齿轮服役工况下不同啮合位置的卷吸速度和相对滑动速度，并将卷吸速度和相对滑动速度作为弹流模型的其他边界条件。

进一步，步骤7中，所述理论计算结果包括采用FORTRAN语言编写弹流润滑理论的计算结果和采用道森-希金森膜厚经验公式的计算结果。

进一步，步骤7中，所述各工况包括从啮入到啮出的八个啮合位置的弹流润滑仿真模型。

本发明的有益效果在于：

本发明通过MATLAB将真实粗糙度数据写入，考虑了粗糙度效应，从而更加真实地模拟了齿轮接触界面三维形貌弹流润滑，并通过油膜厚度和摩擦系数实验验证了仿真方法的正确性，实现了闭环工作，使本发明的可靠性得到保障，为其他复杂弹流模型的分析方法提供了方法依据。

总之，本发明创造性地采用流固耦合的方法针对齿轮齿面不同啮合位置精准建模并进行弹流润滑分析，更加真实地反映齿轮齿面不同啮合位置润滑的实际情况，提供有效可行的的弹流润滑设计方案，为提供高质量的喷油润滑设计方案提供帮助。

附图说明

图1为本发明的流程图；

图2为齿轮啮合某一时刻弹流润滑的等效模型；

图3为弹流模型编号示意图；

图4为渐开线的形成示意图；

图5为齿轮不同啮合位置综合曲率半径示意图；

图6为图2的局部有限元模型；

图7为图6结构模型的边界条件设置图；

图8为实验测量真实情况下齿面粗糙度形貌图；

图9为MATLAB处理后的齿面粗糙度形貌图；

图10为齿轮结构动力学仿真齿面载荷形貌图；

图11为某一啮合位置速度分解图；

图12为图6中流体模型的边界条件设置图；

图13为三维油膜压力图；

图14为三维油膜厚度图；

图15为中心线上的油膜压力和油膜厚度图；

图16为油膜厚度等高线图；

图17为弹流润滑理论计算油膜压力和油膜厚度曲线图；

图18为齿面不同啮合位置仿真最小膜厚与经验最小膜厚对比图；

图19为带粗糙度的油膜厚度等高图；

图20为齿轮齿面不同啮合位置光滑弹流模型摩擦系数折线图；

图21为齿轮齿面啮入和啮出侧粗糙弹流模型最小膜厚折线图；

图22为齿轮齿面啮入和啮出侧粗糙弹流模型摩擦系数折线图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

在本发明的上述描述中，需要说明的是，术语“一侧”、“另一侧”等指示的方位或位置关系为基于附图所示的方位或位置关系，或者是该发明产品使用时惯常摆放的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性。

此外，术语“相同”等术语并不表示要求部件绝对相同，而是可以存在微小的差异。术语“垂直”仅仅是指部件之间的位置关系相对“平行”而言更加垂直，并不是表示该结构一定要完全垂直，而是可以稍微倾斜。

本实施例提供一种基于热流固耦合的直齿轮弹流润滑仿真分析方法，如图1所示，具体包括以下步骤：

1、针对齿轮啮合位置进行模型等效

以直齿轮啮合过程为基础，将某一时刻的直齿轮接触等效成当量弹性圆柱体I与当量弹性圆柱体II通过油膜接触，如图2所示，当量弹性圆柱体I简称为接触体I，当量弹性圆柱体II简称为接触体II,R1为接触体I的曲率半径，R2为接触体II的曲率半径。

2、针对齿轮齿面不同啮合位置精准建模

根据线接触等效后的尺寸建立从啮入到啮出不同啮合位置的多弹流模型弹流模型，如图3，模型所需要的综合曲率半径通过解析法计算，计算方法如下:

如图4所示，发生线所示位置I-I，运动到位置Ⅱ-Ⅱ，发生线上点A的运动轨迹AK即为渐开线，K点所对应的曲率半径的大小即为

在求解精确建立热弹流仿真模型时，需要用到曲率半径，计算如下：

ρ＝r

其中：r

因此，啮合点K处的综合曲率半径为：

将齿轮参数代入上式可计算出从啮入到啮出的任意啮合位置处两齿轮接触点曲率半径及对应的综合曲率半径。用MATLAB绘制齿面任意啮合位置处的综合曲率半径分布情况如图5示。

根据结构动力学仿真获得接触区域的接触载荷，如下表

根据赫兹接触接触半宽公式为：

其中

E1，E2为接触模型I，和接触模型Ⅱ的弹性模量

w为线载荷

u1和u2为接触模型I，和接触模型Ⅱ的泊松比

ρ1和ρ2为接触模型I，和接触模型Ⅱ的接触曲率半径

接触模型的尺寸一般为赫兹接触接触宽度的三到五倍。

接触模型上述所述弹流模型包括结构模型和流体模型；

接触体I 1500um×1500um×1500um；

接触体II 1500um×1500um×1500um；

流体1500um×1500um；

初始最小油膜厚度为1.5um，如图6所示。

3、ADINA绘制模型网格

根据模型整体尺寸大小合理选择网格尺寸大小，此处网格尺寸大小为15um，靠近流固耦合面网格加密。

由于圆柱体为弹性材料，故圆柱体面绘制三层密网格，其余区域用一层网格代替，网格细节部分如图6所示。

流体模型的网格为四面体，方便在ADINA中设置自适应动网格。

4、将真实齿轮齿面粗糙度写入模型

采用高集成多功能摩擦磨损试验台，收集样件表面的三维微观形貌数据，测量齿轮齿面真实粗糙度数据，如图8所示，通过matlab将粗糙度数据处理成符合模型施加条件的数据点，如图9所示。

使用MATLAB读入模型文件，并把粗糙度数据以坐标形式写入流固耦合面的节点上，保证流固耦合面上的节点坐标一致，保存为新的网格模型文件。

5、根据齿轮真实工况设置ADINA边界条件

对服役工况下的齿轮进行结构动力学仿真，获取齿面载荷分布，如图10。再通过解析法获取齿轮啮合过程中任意啮合位置的相对滑动速度，绝对滑动速度，具体如下：

如图11所示，啮合点K处的速度大小为v

v

v

N

v

其中，ω

v

由下式可求得啮合点K处的相对切向速度大小为：

v

综上，结合啮合特性，计算出模型在服役工况下，任意啮合位置的齿面相对滑动速度，如下表所示，其中R1为主动轮接触曲率半径，R2为从动轮接触曲率半径，R为接触综合曲率半径，V1为主动轮绝对滑动速度，V2为从动轮绝对滑动速度，Vt为相对滑动速度，SSR为滑滚比。

不同转速工况齿面啮合点边界条件参数表

在ADINA中进行上述提供边界条件的设置，如图12所示。在流体模型中选择设置自适应动网格。分别设置求解控制参数，采用direct耦合方式进行流固耦合数值计算。

7、模型后处理

用MATLAB读模型计算结果，并绘制三维油膜压力图、三维油膜厚度图及中心线上的油膜压力和油膜厚度图,分别如图13、图14及图15所示。

同时，油膜厚度的等高线图，如图16所示，在三维线接触弹流润滑仿真计算结果中，油膜厚度和油膜压力在接触线方向上均匀分布；在啮合方向的中心区域处，油膜厚度分布平坦、数值较小，油膜内部出现中心压力峰，油膜压力数值较大；在啮合方向的近出口区域处，油膜厚度发生局部“颈缩”现象，“颈缩”程度沿接触线方向一致。

由上述结果可知，在出口位置出现了压力二次高峰和油膜颈缩现象，这是经典的弹流润滑特征，说明该仿真方案与弹流润滑十分吻合。

8、验证与归纳

取环境温度90℃，首先对模型进行基本假设，然后使用FORTRAN语言编写弹流润滑理论计算公式，进行齿轮单个啮合点的弹流润滑理论计算。

1)基本假设

为减小计算难度，把单个啮合点的弹流润滑理论计算模型简化为二维的线接触弹流润滑模型，进行以下基本假设：

润滑油为牛顿流体，流动特性遵守牛顿粘性定律；

润滑油的惯性力及体积力忽略不计；

压力沿膜厚方向不变，润滑油沿膜厚方向速度不计；

油膜厚度方向上的油液粘度和密度为常数；

剪切力和速度梯度只在油膜厚度方向上比较显著；

固体表面光滑，忽略表面粗糙度的影响；

润滑油的流动处于层流状态；

与油膜相邻的两侧固体表面曲率半径远大于油膜厚度，固忽略油膜曲率的影响；

流体区域和两侧固体区域是等温的。

2)基本方程

基本方程包括：Reynolds方程、膜厚方程、变形方程、粘压方程、密压方程、载荷方程，主要形式如下所示。

Reynolds方程：

式中，x为啮合方向，ρ为润滑油密度，h为油膜厚度，η为润滑油粘度，p为油膜压力，u

膜厚方程：

式中，R为综合曲率半径，υ为表面弹性变形，由Boussinesq积分给出，如变形方程。变形方程:

式中，E为综合弹性模量，

粘压方程:

式中，η

密压方程:

式中，ρ

载荷方程:

3)理论结算结果对比

线接触弹流润滑问题的计算方法中，粘压方程、密压方程和变形方程都与压力有关，因此一般也是先给定一个初始的压力分布，如Hertz接触压力。计算得润滑油粘度、密度和弹性变形后，通过膜厚方程获得膜厚值。然后把膜厚代入Reynolds方程求解出新的压力分布，再不断对前一次的压力分布进行迭代修正、重新计算弹性变形和膜厚值，反复迭代，直至两次迭代得到的压力差十分接近，迭代结束。最终获得油膜压力和油膜厚度分布。

FORTRAN理论计算结果如图17所示，FORTRAN理论计算所得到的中心油膜压力为813.458Mpa，中心油膜厚度796.1nm，最小油膜厚度为786.5nm。

由图15和图17可知，从油膜厚度和压力分布形态来看，仿真计算结果同理论计算结果相似。对比FORTRAN理论计算结果和仿真计算结果，从具体数值来看，仿真计算结果同理论计算结果有一定差异，具体如下表所示。

啮入点处理论与仿真计算油膜压力、油膜厚度的对比

由表可知，在油膜压力计算方面，仿真计算的中心油膜压力为807.049Mpa，与理论计算的中心油膜压力813.458Mpa相差不大；在油膜厚度计算方面，仿真计算的中心油膜厚度为676.6nm，与理论计算的中心油膜厚度796.1nm相差119.5nm，仿真计算的最小油膜厚度为575.2nm，与理论计算的最小油膜厚度726.5nm相差151.3nm。由此看来，仿真计算的油膜压力近似于理论计算结果，但仿真计算的油膜厚度比理论计算结果小100多纳米。

分析原因如下：数值理论计算方法把线接触弹流润滑问题简化为二维的润滑问题，采用基于多个假设的Reynolds方程进行弹流润滑计算，完全忽略了润滑油油膜厚度方向上流体属性(速度梯度、粘度梯度和压力梯度)的变化，而仿真方法不仅可以完整模拟线接触弹流润滑问题的三维接触特性，而且可以使用完整的N-S动量守恒方程进行流体运动分析，考虑了润滑油的惯性力、体积力，还考虑了膜厚方向上润滑油粘度和密度的差异，使计算结果更加准确。

而本专利研究对象为高速重载直齿轮，使用带假设的二维数值理论方法计算的油膜厚度必然也有偏大的趋势，由此看出，本文提出的三维齿轮弹流润滑仿真方法的计算准确性高于二维数值理论方法。

4)道森-希金森膜厚公式

道森-希金森公式适用于重载荷、高弹性模数及高粘压系数的情况。

其中

R为综合曲率半径

w为线载荷

E′为综合弹性模量

η

图18为齿面啮合过程中不同啮合位置的仿真最小油膜厚度和道森-希金森经验公式计算的最小油膜厚度对比图，由图可知，仿真最小膜厚和经验最小膜厚整体上趋势一致，但仿真经验最小膜厚要略大于仿真最小膜厚，这是由于在重载荷和高滑动速度的情况下，经验公式中影响因素的系数可能会出现偏差，尤其是考虑环境温度的情况下，实际膜厚是要偏小的。

带粗糙度效应的仿真结果：图19带粗糙度的油膜厚度等高高图。由结果可知，在粗糙峰处，明显出现了油膜扭曲与断裂分离的现象，另一方面，粗糙度在一定程度上使得接触线区域出现峰顶碰撞的现象，符合真实情况下混合润滑下的摩擦情况。

5)分析各个工况不同啮合位置仿真结果

下面分析各个工况不同啮合位置的油膜、油压、摩擦系数。

通过计算从啮入到啮出的八个啮合位置的光滑表面弹流润滑仿真模型，提取每个啮合位置仿真模型的摩擦系数和对应的滑滚比SSR，并绘制成曲线图，如图20所示，可以看出靠近啮入和啮出位置的摩擦系数相对更大。

因此，单独提取啮入和啮出侧的啮合位置进行进一步仿真分析探究，通过matlab给两个位置施加粗糙度，保存为新弹流润滑仿真模型。考虑到真实情况下齿轮喷油润滑过程受齿轮转动风阻、喷油压力以及喷油量等多个因素影响，使得齿轮齿面啮合过程中并非处于充足供油的状态，因此要考虑不充分供油情况下的润滑及摩擦情况。油气比例的计算是滑油和空气混合后会有不同的粘度和密度，用按比例变化流体材料属性的方法来模拟不同油气比例的情况，此处认定油气比例为1时，供油充足无空气。

这里设定当油气比例为1时，供油充足，齿面完全覆盖上油膜。因此下一步进行不同油气比例下(即不同供油情况下)齿面弹流润滑仿真。

如图21，研究对象仍为靠近啮入和啮出的位置，提取不同油气比例下带粗糙度的弹流润滑仿真结果流体部分的油膜厚度，并与光滑表面的弹流润滑仿真结果做对比，结果可知在油气比例高的情况下，光滑面的最小膜厚相对粗糙面要大，反之，在油气比例低的情况下，光滑面的最小膜厚相对粗糙面要小。

如图22，提取不同油气比例下带粗糙度的弹流润滑仿真结果摩擦系数，并与光滑表面的弹流润滑仿真结果做对比，结果可知在油气比例高的情况下，带粗糙度的弹流润滑仿真模型，摩擦系数要远大于光滑表面的摩擦系数，尤其是当油气比例降低的情况下，摩擦系数升高的速度明显提升。

因此根据以上分析，可以为齿轮啮合界面喷油参数如喷油量、喷射角度、喷射位置，以及最佳喷油效果提出量化的指标，为齿轮啮合高效润滑冷却喷油润滑设计方案提供帮助。

从上述结果可以看出，本弹流模型仿真方法可以较为准确模拟齿轮不同啮合位置线接触的油膜厚度和摩擦系数，并可以通过MATLAB粗糙度数据写入来模拟实际情况下的线接触弹流润滑情况，且针对不同啮合位置不同供油情况进行分析，为齿轮啮合高效润滑冷却喷油润滑设计方案提供帮助。

本实施例通过对ADINA的材料模型进行二次开发，考虑了热效应，流变效应，通过MATLAB写入粗糙度数据，考虑了真实齿轮齿面粗糙度效应，从而有效地模拟了实际情况下的弹流润滑，并通过FORTRAN语言编写弹流润滑理论计算公式和道森-希金森膜厚公式验证了仿真方法的可行性和正确性，为其他复杂弹流模型的分析方法提供了方法依据。同时，针对不同啮合位置不同供油情况进行分析，为齿轮啮合高效润滑冷却喷油润滑设计方案提供帮助。

相较于理论计算方法，本实施例不但计算效率高，展示的结果丰富，而且适用性广，可广泛用于复杂齿轮弹流线接触模型中。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。