CN104091022A - 一种速度滑移工况下静压主轴系统流固耦合仿真分析方法 - Google Patents

Abstract

本发明涉及一种速度滑移工况下静压主轴系统流固耦合仿真分析方法，该方法包括：选择Navier速度滑移模型模拟轴承流固界面速度滑移现象，建立考虑速度滑移情况下静压轴承油膜压力分布的有限差分模型，求解油膜的压力分布；借助ansys有限元软件建立转子和轴瓦及轴承座的有限元模型，应用APDL参数化有限元分析技术对转子和轴瓦及轴承座有限元模型求解；提取转子和轴瓦及轴承座有限元模型相关节点的变形值，并将数据导入到油膜厚度的计算中，实现对油膜压力的重新求解；重复上述过程直至转子和轴瓦及轴承座的变形达到稳定。本发明考虑了固流界面速度滑移的影响，实现了静压轴承-转子系统的流体-结构之间的双向耦合。

Description

一种速度滑移工况下静压主轴系统流固耦合仿真分析方法

技术领域

本发明流体动力学计算领域，涉及一种适用于静压轴承-转子系统流体-结构之间的双向耦合仿真方法。

背景技术

静压轴承以支撑精度高等优点，广泛应用于机械制造和仪器制造等领域。在静压轴承中，润滑介质在轴承间隙中的流动已经进入到了微尺度流动的研究领域。然而流体在微尺度下流动时将会表现出异于宏观尺度的流动的性质，流体在微尺度下会表现出来大表面积体积比、梯度参数效应及界面效应等问题。将微尺度的性质引入到静压轴承的研究中，从而使油膜性能的研究更加准确。

静压轴承性能分析在很大程度上依赖于仿真方法，仿真方法的优劣直接决定静压轴承分析的结果的准确性。但传统的仿真方法是对静压轴承工作环境做出在一定程度上合理的假设，传统的方法忽略了在静压轴承流固界面的速度滑移现象，而且没有考虑轴承和转子之间的相互作用，没有对轴承和转子之间进行流固耦合分析，速度滑移和流固界面的相互耦合作用在一定程度上影响着油膜行为，油膜行为对决定静压轴承性能起着关键作用。在以往的一些研究中单纯研究速度滑移对轴承承载力的影响，静压轴承受滑移的影响承载力会在一定程度上得到提高，单纯研究流固界面流体和固体的相互耦合作用对静压轴承的影响，轴承的承载力在一定程度上会有所降低，而在静压轴承的实际工况作用下，两种作用都会对静压轴承产生一定的影响，而两种影响作用不能简单的相互抵消，最终的影响结果和很多因素有关，例如轴瓦及转子材料、油膜温度等，这必然会使传统的仿真方法的结果和实际工作结果之间产生一定的误差。

发明内容

针对现有技术中存在的上述问题，本发明提供一种考虑速度滑移影响下静压轴承-转子系统的流体-结构耦合仿真方法，考虑了固-流界面速度滑移的影响，具有计算简洁、通用性强、可人工控制耦合迭代过程等优点，实现了真正意义上的静压轴承-转子系统的流体-结构之间的双向耦合。

一种速度滑移工况下静压主轴系统流固耦合仿真分析方法，包括以下步骤：

(1)根据静压轴承的工况选择Navier速度滑移模型，模拟轴承-转子系统中流固界面上出现的速度滑移现象。

Navier模型假设流体的滑移速度v_s与流体在表面的剪切力成正比，即：

$v_s=b\frac{{\∂u}_1}{\∂z}{|}_{\mathrm{wall}}---\left(1\right)$

式中，u₁是流体沿固体表面的流动速度，Z是沿界面法向方向坐标，wall表示流体和固体相互作用的界面位置。b称之为滑移系数，是虚构固体表面与实际界面的距离，如图1所示。

(2)建立考虑速度滑移情况下静压轴承油膜压力分布的有限差分模型，根据油膜结构划分合理网格。

考虑速度滑移条件下的油膜压力分布符合修正后的N-S方程为：

$\frac{\∂}{\∂x}\left\{\right[-\frac{h^3}{12\mathrm{\η}}-\frac{{\mathrm{hb}}_{\mathrm{bx}}}{2\mathrm{\η}}(h-h^{*})-\frac{{\mathrm{hh}}^{*}{b}_{\mathrm{ax}}}{2\mathrm{\η}}]\frac{\∂p}{\∂x}+\frac{\mathrm{hu}}{2}\}+\frac{\∂}{\∂z}\left\{(-\frac{h^3}{12\mathrm{\η}}-\frac{h^2{b}_{\mathrm{az}}}{2\mathrm{\η}})\frac{\∂p}{\∂z}\right\}=0---\left(2\right)$

式中，u为旋转主轴表面的线速度；b_ax、b_az分别为轴瓦表面周向、轴向滑移系数；b_bx为轴颈表面周向滑移系数；h^*为油膜周向方向剪切力为零时的油膜厚度；p为油膜压力；η为润滑油的动力粘度。

采用有限差分方法将方程(2)线性化处理，得到如下线性化方程：

$P_{i,j}=\frac{1}{{2A}_{i,j}+2C_{i,j}}[(A_{i,j}+B_{i,j})P_{i+1,j}+(A_{i,j}-B_{i,j})P_{i-1,j}+C_{i,j}{P}_{i,j+1}+C_{i,j}{P}_{i,j-1}-D_{i,j}{P}_{i,j}]---\left(3\right)$

其中，

$C_{i,j}=\frac{{\mathrm{\β}}^2{H_{i,j}}^3(1+{6b}_0)}{{\mathrm{\Δ\λ}}^2}$

$\mathrm{\β}=\frac{2R}{L}$

式中，b₀为滑移系数，R为轴承半径，L为轴承长度，Δλ为轴承轴向方向网格长度，为轴承周向方向网格长度，(i，j)为油膜位置坐标，H_i,j表示(i，j)处无量纲油膜厚度。

应用超松弛迭代法将线性化方程(3)在matlab软件平台上进行数值求解，即得考虑速度滑移影响下油膜压力分布。

(3)用matlab语言编写程序将计算所得的压力分布数据导出，将数据文件存储成txt格式，便于将此计算数据在步骤(5)中作为压力边界条件加载到转子和轴瓦及轴承座有限元模型中。

(4)借助ansys有限元软件建立转子和轴瓦及轴承座的有限元模型，控制有限元模型的网格划分形式与油膜的网格划分一致，以保证数据交换过程的顺利完成。

在转子模型建立过程中，采用自下而上的建模方法，先创建关键点，依次构成线和面，用lsize控制转子有限元模型的网格划分，将面绕轴线旋转出周向规则的网格，然后在转子模型的表面建立表面效应单元，用于加载与表面相互垂直的压力载荷，定义转子有限元模型的材料属性，转子的有限元模型建立完成。

在轴瓦及轴承座模型的建立中，静压腔的深度与轴承的尺寸相比很小，所以采用有限元方法对轴瓦建模时忽略轴瓦表面的静压腔，有限元模型的建立过程与转子模型的建立过程相同，但需根据实际情况分别对轴瓦及轴承座定义有限元模型的材料属性。

(5)应用APDL参数化有限元分析技术编写相应程序实现matlab与ansys之间的数据交换，编写APDL程序将matlab数据读入ansys，将油膜压力分布的matlab计算结果作为边界条件加载到转子和轴瓦及轴承座有限元模型中，实现对转子和轴瓦及轴承座有限元模型的求解。

(6)应用APDL参数化有限元分析技术提取转子和轴瓦及轴承座有限元模型相关节点的变形值。有限元模型在某节点的变形可以分解为X,Y,Z三个方向，节点在X,Y方向上的变形将影响到油膜厚度的变化。

(7)将ansys导出的节点变形数据导入到matlab油膜厚度的计算中，实现对油膜压力分布的重新求解。

(8)重复步骤(3)～(7)直至转子和轴瓦及轴承座的变形达到稳定，实现真正意义上静压轴承-转子系统的流体-结构之间的双向耦合。

与现有技术相比，本发明具有以下优点：

本发明考虑了传统设计中忽略的固液界面速度滑移现象，将油膜分析尺度引入到微尺度的研究领域；并将流体和固体分析过程进行耦合，使得分析结果更加准确，更加有利于指导静压轴承应用于实践中。应用本发明仿真方法和传统仿真方法在不忽略速度滑移和不进行流固耦合分析的情况下计算某静压轴承的承载力，并与实验方法得到的承载力进行比较，本发明仿真方法计算得到承载力误差为2.97％，而传统仿真方法得到的误差分别为10.89％和38.6％。实验表明，本发明仿真方法的精度明显高于传统仿真方法。

附图说明

图1为Navier滑移模型示意图；

图2为本发明所涉及方法的流程图。

具体实施方式

本发明所述方法由安装在计算机上的软件程序实现。所述计算机上安装由美国TheMathWorks公司出品的商业数学软件MATLAB和由美国ANSYS公司开发Mechanical APDLProduct Lanucher软件。

本发明所述方法的的流程图如图2所示，具体包括以下步骤：

步骤1，根据静压轴承的工况选择Navier速度滑移模型，模拟轴承-转子系统中流固界面上出现的速度滑移现象。

步骤2，建立考虑速度滑移情况下静压轴承油膜压力分布的有限差分模型，根据油膜结构划分合理网格。

步骤3，用matlab语言编写程序将计算所得的压力分布数据导出，将数据文件存储成txt格式，便于将此计算数据作为压力边界条件加载到转子和轴瓦及轴承座有限元模型中。

步骤4，借助ansys有限元软件建立转子和轴瓦及轴承座的有限元模型，控制有限元模型的网格划分形式与油膜的网格划分一致，以保证数据交换过程的顺利完成。

步骤5，应用APDL参数化有限元分析技术编写相应程序实现matlab与ansys之间的数据交换，编写APDL程序将matlab数据读入ansys，将油膜压力分布的matlab计算结果作为边界条件加载到转子和轴瓦及轴承座有限元模型中，实现对转子和轴瓦及轴承座有限元模型的求解。

步骤6，应用APDL参数化有限元分析技术提取转子和轴瓦及轴承座有限元模型相关节点的变形值。有限元模型在某节点的变形可以分解为X,Y,Z三个方向，节点在X,Y方向上的变形将影响到油膜厚度的变化。

步骤7，将ansys导出的节点变形数据导入到matlab油膜厚度的计算中，实现对油膜压力分布的重新求解。

步骤8，重复步骤3～7，直至转子和轴瓦及轴承座的变形达到稳定，实现了真正意义上静压轴承-转子系统的流体-结构之间的双向耦合。

表1给出了三种不同条件下通过仿真计算得到的一定结构的静压轴承承载力，并列出通过一定的实验方法得到的静压轴承承载力数据。三种条件分别为：无滑移并不考虑流固界面的相互耦合；考虑速度滑移但不进行流固界面的相互耦合分析；考虑速度滑移并进行流固界面的相互耦合分析。表1同时还给出了实验方法得出的静压轴承的承载力，即采用测微仪检测主轴偏心量，并通过砝码对轴承施加静力外载荷，当主轴的偏心率为0.1时静压轴承的承载力。传统的仿真方法一般忽略静压轴承流固界面的速度滑移，同时也不会进行流固耦合分析，本发明仿真方法将两者都作为影响静压轴承性能的因素并考虑在仿真过程中。从表1中可以看出实验测得的静压轴承的承载力为1010N，通过传统仿真方法在不忽略速度滑移和不进行流固耦合分析的情况下轴承的承载力为1120N，考虑速度滑移但不对流固界面的相互耦合作用进行分析时，静压轴承的承载力为1400N，而采用本发明仿真方法得到的静压轴承的承载力为980N。通过计算可得三种不同条件下得到的承载力相对实验结果的误差分别为10.89％、38.6％和2.97％，因此，本发明仿真方法的准确度明显高于传统方法。

表1 不同条件下静压轴承的承载力

Claims (3)

1.一种速度滑移工况下静压主轴系统流固耦合仿真分析方法，其特征在于包括以下步骤：

步骤1，根据静压轴承的工况选择Navier速度滑移模型，模拟轴承-转子系统中流固界面上出现的速度滑移现象；

Navier模型假设流体的滑移速度vs与流体在表面的剪切力成正比，即：

$v_s=b\frac{{\∂u}_1}{\∂z}{|}_{\mathrm{wall}}---\left(1\right)$

式中，u₁是流体沿固体表面的流动速度，Z是沿界面法向方向坐标，wall表示流体和固体相互作用的界面位置；b称之为滑移系数，是虚构固体表面与实际界面的距离；

步骤2，建立考虑速度滑移情况下静压轴承油膜压力分布的有限差分模型，根据油膜结构划分合理网格；

考虑速度滑移条件下的油膜压力分布符合修正后的N-S方程为：

$\frac{\∂}{\∂x}\left\{\right[-\frac{h^3}{12\mathrm{\η}}-\frac{{\mathrm{hb}}_{\mathrm{bx}}}{2\mathrm{\η}}(h-h^{*})-\frac{{\mathrm{hh}}^{*}{b}_{\mathrm{ax}}}{2\mathrm{\η}}]\frac{\∂p}{\∂x}+\frac{\mathrm{hu}}{2}\}+\frac{\∂}{\∂z}\left\{(-\frac{h^3}{12\mathrm{\η}}-\frac{h^2{b}_{\mathrm{az}}}{2\mathrm{\η}})\frac{\∂p}{\∂z}\right\}=0---\left(2\right)$

式中，u为旋转主轴表面的线速度；b_ax、b_az分别为轴瓦表面周向、轴向滑移系数；b_bx为轴颈表面周向滑移系数；h^*为油膜周向方向剪切力为零时的油膜厚度；p为油膜压力；η为润滑油的动力粘度；

采用有限差分方法将方程(2)线性化处理，得到如下线性化方程：

$P_{i,j}=\frac{1}{{2A}_{i,j}+2C_{i,j}}[(A_{i,j}+B_{i,j})P_{i+1,j}+(A_{i,j}-B_{i,j})P_{i-1,j}+C_{i,j}{P}_{i,j+1}+C_{i,j}{P}_{i,j-1}-D_{i,j}{P}_{i,j}]---\left(3\right)$

其中，

$C_{i,j}=\frac{{\mathrm{\β}}^2{H_{i,j}}^3(1+{6b}_0)}{{\mathrm{\Δ\λ}}^2}$

$\mathrm{\β}=\frac{2R}{L}$

式中，b₀为滑移系数，R为轴承半径，L为轴承长度，Δλ为轴承轴向方向网格长度，为轴承周向方向网格长度，(i，j)为油膜位置坐标，H_i,j表示(i，j)处无量纲油膜厚度；应用超松弛迭代法将线性化方程(3)在matlab软件平台上进行数值求解，即得考虑速度滑移影响下油膜压力分布；

步骤3，用matlab语言编写程序将计算所得的压力分布数据导出，将数据文件存储成txt格式，便于将此计算数据在步骤5中作为压力边界条件加载到转子和轴瓦及轴承座有限元模型中；

步骤4，借助ansys有限元软件建立转子和轴瓦及轴承座的有限元模型，控制有限元模型的网格划分形式与油膜的网格划分一致，以保证数据交换过程的顺利完成；

步骤5，应用APDL参数化有限元分析技术编写相应程序实现matlab与ansys之间的数据交换，编写APDL程序将matlab数据读入ansys，将油膜压力分布的matlab计算结果作为边界条件加载到转子和轴瓦及轴承座有限元模型中，实现对转子和轴瓦及轴承座有限元模型的求解；

步骤6，应用APDL参数化有限元分析技术提取转子和轴瓦及轴承座有限元模型相关节点的变形值；有限元模型在某节点的变形可以分解为X、Y、Z三个方向，节点在X、Y方向上的变形将影响到油膜厚度的变化；

步骤7，将ansys导出的节点变形数据导入到matlab油膜厚度的计算中，实现对油膜压力分布的重新求解；

步骤8，重复步骤3～7，直至转子和轴瓦及轴承座的变形达到稳定，实现真正意义上静压轴承-转子系统的流体-结构之间的双向耦合。

2.根据权利要求1所述的一种速度滑移工况下静压主轴系统流固耦合仿真分析方法，其特征在于，在步骤4所述转子模型建立过程中，采用自下而上的建模方法，先创建关键点，依次构成线和面，用lsize控制转子有限元模型的网格划分，将面绕轴线旋转出周向规则的网格，然后在转子模型的表面建立表面效应单元，用于加载与表面相互垂直的压力载荷，定义转子有限元模型的材料属性，转子的有限元模型建立完成。

3.根据权利要求1所述的一种速度滑移工况下静压主轴系统流固耦合仿真分析方法，其特征在于，在步骤4所述轴瓦及轴承座模型的建立中，静压腔的深度与轴承的尺寸相比很小，所以采用有限元方法对轴瓦建模时忽略轴瓦表面的静压腔，有限元模型的建立过程与转子模型的建立过程相同，但需根据实际情况分别对轴瓦及轴承座定义有限元模型的材料属性。