一种基于Fluent的径向油膜轴承的动网格更新方法与流程

本发明属于油膜轴承及其转子系统性能计算技术领域，具体涉及一种基于fluent的径向油膜轴承的动网格更新方法。

背景技术：

进行油膜轴承动态性能计算时，转子的轴颈中心在涡动，轴颈表面各点的绝对速度和加速度不一样，计算流场的形状和运动学参数时刻在改变，利用静网格模型不能进行计算。动网格模型可以用来模拟由于边界运动而引起的流域随时间的变化，在单相流和多相流下均适用。边界运动可以是规定动作（例如，用户可以定义刚体中心的速度和角速度），也可以是需要待求解的动作（例如，六自由度模型下以求解结果确定下一时刻的速度、加速度、角速度、角加速度）。对于动网格模型，用户需要确定初始网格和运动区域。常用的流体性能计算软件fluent提供了三种动网格更新算法：smoothing、layering、remeshing。根据运动边界的变化，fluent软件对指定的网格自动更新。smoothing方法即为光顺法，特点是对网格节点间的拓扑关系不改变，只进行网格形状的改变，当网格出现畸变时容易出现发散。layering方法是随着动边界的移动，在边界处发生网格的增加或合并。remeshing方法是将控制区内的所有网格重新划分，一般适用于非结构化网格。

考虑到油膜径向尺寸很小的特点，从理论上来说，smoothing动网格更新方式更适合用于轴承性能计算中，但会出现网格畸变。曾有人提出了基于弹性变形的动网格更新方法。该方法与smoothing方法中的laplace方法类似。网格节点位移按以下公式计算：

（1）

（2）

式中，为膜厚方向网格的层数，为网格节点所在的层数，和是轴颈中心的位移。

基于弹性变形的动网格调整法存在两个问题：首先，一旦偏心率过大，垂直于轴颈的网格线就会出现严重倾斜，这种网格更新算法容易出现发散、停止、闪退等现象，导致计算结果不准确；第二，这样的计算方式容易造成累积误差，长时间计算会有网格畸变甚至负体积的情况。

纵观流场计算的动网格算法用于油膜轴承性能计算时，由于油膜轴承径向尺寸较小，网格发生变形时，极容易出现网格畸变和负体积网格，导致轴承性能无法计算或者累积误差较大。

技术实现要素：

本发明为了解决现有技术中存在轴承油膜网格长宽比过大容易在动网格中发生网格畸变的问题，提供一种基于fluent的径向油膜轴承的动网格更新方法。

为解决上述技术问题，本发明采用如下技术方案：一种基于fluent的径向油膜轴承的动网格更新方法，包括以下步骤：

（1）创建轴承油膜网格文件；

（2）fluent软件中udf控制宏程序的编写：包括三个宏：定义轴颈表面各点速度的宏define_profile、自定义动网格更新方式的宏define_grid_motion和用于创建网格节点数据文本文件的自定义宏define_grid_motion；

（3）启动流体计算软件fluent，导入网格文件，挂载用于创建网格节点数据文本的宏define_grid_motion，得到网格节点数据文件；

（4）利用excel工具对步骤（3）中得到的数据文件进行处理；

（5）启动流体计算软件fluent，导入网格文件，挂载udf程序并设置udf控制条件；

（6）利用udf动网格程序控制流场计算；

（7）判断轴颈中心轨迹是否收敛。

所述的步骤（1）具体为：根据研究对象选取有限单元的类型和网格划分方法，借助于网格划分软件gambit将计算流场划分为混合网格，对于轴承油膜的进油口和进油槽部分采用非结构化网格，对油膜部分采用结构化网格。

所述的步骤（2）中，由于轴颈中心存在变位速度，其坐标在时刻改变，所以轴颈表面各点的速度根据相对运动的公式进行计算，计算出各点的速度数值不一样，通过编写宏define_profile程序控制轴颈表面各点的速度。

所述的步骤（2）中，利用宏define_grid_motion进行动网格程序的更新编写，控制指向轴颈中心的网格线方向始终不改变，从而避免网格倾斜，减少计算累计误差。

所述的步骤（2）中，在fluent软件中，利用自定义动网格技术计算轴承油膜特性时，需要建立4个文本文件以记录每个时间步的动力学数据，以便于用于下一时间步的计算；文件1记录轴颈受的油膜力；文件2记录每一时间步的轴颈中心的坐标；文件3记录每一时间步的轴颈中心的变位速度；文件4记录轴瓦的坐标。

所述的步骤（3）中，在fluent软件中，要精确控制油膜中的所有网格节点，需要获取网格节点的全局编号、所在层数和相对于轴颈中心的单位向量，并把这些数据保存在文件5中，以方便查找和控制；以上五个文本文件需要利用宏define_grid_motion程序进行编写，在动网格预览中执行已编译的程序，得到所需的上述网格节点文本文件。

所述的步骤（5）中，控制条件的设置分别为：

problemsetup设置：选择mixture模型，定义润滑油物性参数，液相与气相之间的质量传递选择singhal气化方式，设置进油口压力、出油口压力以及操作压力，转速选择udf中已定义的宏define_profile，动网格选择udf中已定义的宏define_grid_motion；

solution设置：压力-速度耦合项选择piso，监视器中的连续性残差、速度残差、气相残差均设置，监视进油口和出油口的质量流量，其余均保持默认设置；初始化流场置零；

时间步长设置为0.0001~0.0005，时间步数为10000；其余默认，点击calculate计算。

使网格不发生扭曲变形，就需要指向轴颈中心的网格线方向始终不改变，始终垂直于轴颈表面，同时要保证在移动前后径向方向的网格线沿圆周方向上等均分布，使径向网格线的夹角相等，即移动前的网格线与移动后的网格线共线或平行；所述的步骤（6）中的计算方法具体为：

假设移动前和移动后轴承的中心点分别为和，移动前和移动后的轴颈的中心点为和，油膜网格在径向方向上分成层网格，假设为时刻第层网格上任一待更新节点，为节点更新后对应的第层网格节点；根据此刻的轴颈受的油膜力、轴颈坐标、轴颈中心变位速度、轴颈质量可以得出下一时刻的轴颈中心坐标，根据轴颈中心坐标和该节点相对轴颈中心的单位向量，做射线分别交轴颈圆和轴承圆于点、，根据移动前后网格线共线或平行，可知与平行或共线；又知待求点所在的层数为，即已知，为轴承半径，为轴颈半径；坐标求解算法步骤如下：

（1）、求出；

（2）、已知的单位法向量，由于与共线，所以得出的单位法向量即为；

（3）、求解与的夹角余弦

（4）、利用三角形关系式求解

（5）、求出

（6）、，已知坐标，即可求出坐标。

所述的步骤（5）具体为：首先设定一个初始速度和初始位移，通过udf函数define_grid_motion计算出油膜力和，再根据牛顿第二定律计算出该时间步的轴颈中心坐标、速度；如果需要，再计算出轴瓦的中心坐标以及速度；将计算出来的数据分别导入对应的文本文件：文件1、文件2、文件3和文件4；同时利用文件5精确控制网格，进行网格更新，在该时间步内迭代完后进入下一个时间步；下一个时间步用到文件2、文件3和文件4最后一次保存的数值；以此循环，直至轴颈中心趋近于一定点或者稳定的轨迹；当轴颈中心轨迹趋于一定点，或者轴颈中心轨迹稳定为某一封闭曲线时，判断为收敛，否则继续执行步骤（4）直到稳定为止。

采用上述技术方案，本发明解决了轴承油膜网格长宽比过大容易在动网格运动中发生网格畸变的问题，保证了移动前后径向方向的网格线沿圆周方向上等均分布。利用该方法计算动网格节点坐标的过程中，使网格线始终垂直于轴颈表面，网格不发生扭曲变形，避免了网格计算的累计误差，加快了轴承性能的计算速度。

附图说明

图1是本发明的更新流程图；

图2是本发明中流场网格的更新算法原理示意图；

图3是轴瓦受到正弦激振器激振示意图；

图4是轴承油膜网格图；

图5是轴瓦受到正弦激振的轴颈中心轨迹图。

具体实施方式

结合实施例对本发明的技术方案做进一步说明。

轴瓦受激振器激振控制结构简图如图3所示。利用位移激振器来控制轴瓦的运动，来考察激振器的振幅和相位对转子轴心轨迹的影响。算例的参数见表1。

表1轴承、转子和润滑油物性参数

如图1-4所示，fluent软件动网格更新流程中的网格节点更新是根据轴颈、轴瓦的位移确定的。网格更新方式的核心是求解节点坐标的算法，其网格节点求解原理图如图2所示。基于本发明的一种基于fluent的径向油膜轴承的动网格更新方法，计算油膜轴承所支撑的转子的轴心轨迹，包括以下步骤：

（1）根据研究对象选取有限单元的类型和网格划分方法，借助于网格划分软件gambit将计算流场划分为混合网格，对于轴承油膜的进油口和进油槽部分采用非结构化网格，对油膜部分采用结构化网格。

（2）fluent软件中udf控制宏程序的编写：由于轴颈中心存在变位速度，其坐标在时刻改变，所以轴颈表面各点的速度根据相对运动的公式进行计算，计算出各点的速度数值不一样，通过编写宏define_profile程序控制轴颈表面各点的速度；利用宏define_grid_motion进行动网格程序的更新编写，控制指向轴颈中心的网格线方向始终不改变，从而避免网格倾斜，减少计算累计误差；利用宏define_grid_motion获取网格节点的全局编号、初始坐标以及所在层数。定义轴颈表面各点速度的宏define_profile（命名为rot）、自定义动网格更新方式的宏define_grid_motion（命名为dynamic）和用于创建网格节点数据文本文件的自定义宏define_grid_motion（命名为node）。

（3）设置数据传递文本文件：文件1为force.txt读入每一时间步轴颈受的油膜力；文件2为weiyi.txt，读写每一时间步的轴颈中心的位移；文件3为vel.txt，读写每一时间步的轴颈中心的变位速度；文件4为bearing.txt，读写轴瓦的位移。

（4）打开fluent14.0软件，导入网格模型，选择尺寸为mm，挂载udf.c文件，在动网格中选择所有节点区域，挂载宏node，执行一次网格预览，生成网格节点数据文件node.txt（数据包含节点相对原点的矢量和所在的层数，称为文件5）。

（5）对生成的数据node.txt在excel表格中重新排列（按全局编号从小到大的顺序）。

（6）选择计算模型为mixture，定义物性参数，设置进油压力，定义轴颈边界旋转速度时挂载宏rot，动网格部分的宏选择dynamic，初始化流场，时间步设置为0.0001，点击calculate进行计算。其中动网格的技术核心是提出的动网格更新方法，具体计算方法如下：

参照图2，假设移动前和移动后轴承的中心点分别为和，移动前和移动后的轴颈的中心点为和，油膜网格在径向方向上分成层网格，假设为时刻第层网格上任一待更新节点，为节点。根据此刻的轴颈受的油膜力、轴颈坐标、轴颈中心变位速度、轴颈质量可以得出下一时刻的轴颈中心坐标，根据轴颈中心坐标和该节点相对轴颈中心的单位向量，做射线分别交轴颈圆和轴承圆于点、，根据移动前后网格线共线（平行），可知与平行（共线）。又知待求点所在的层数为，即已知，为轴承半径，为轴颈半径。坐标求解算法步骤如下：

（1）、求出；

（2）、已知的单位法向量，由于与共线，所以得出的单位法向量即为；

（3）、求解与的夹角余弦

（4）、利用三角形关系式求解

（5）、求出

（6）、，已知坐标，即可求出坐标。

（7）计算结束后，可以根据weiyi.txt在matlab中生成轴颈中心轨迹图线，以此判断收敛或稳定性。图5（a）给出了激振器振幅为1微米、不同相位差下的转子轴心轨迹图。图5（b）给出了相位差为120⁰、不同激振幅值下的转子轴心轨迹图。

本实施例并非对本发明的形状、材料、结构等作任何形式上的限制，凡是依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均属于本发明技术方案的保护范围。