CN113268901A - 基于格子Boltzmann动压气体轴承间隙微流动仿真方法 - Google Patents

Abstract

本发明涉及一种基于格子Boltzmann动压气体轴承间隙微流动仿真方法，选用多松弛时间格子Boltzmann模型，然后将动压气体轴承简化为刚性轴承，在其转静子间隙流场内生成贴体网格，并采用广义插值格子玻尔兹曼方法进行求解；此外为了实现微流动的模拟，将松弛时间与Kn数进行关联，同时基于非平衡外推格式，将Langmuir滑移模型引入格子Boltzmann方法(LBM)中，可以满足间隙微流动在滑移区和过渡区的模拟；最后可以得到流场内的速度分布情况以及气膜压力，从流场特征角度分析轴承机理，为更好的设计轴承提供一定的理论基础。克服了转静子间隙气体流动处于过渡区流态时，基于连续性假设的常规宏观方法不再适用的问题。

Description

基于格子Boltzmann动压气体轴承间隙微流动仿真方法

技术领域

本发明涉及动压气体轴承仿真研究技术领域，尤其是一种基于格子Boltzmann动压气体轴承间隙微流动仿真方法。

背景技术

气体轴承是指借助于轴承滑动副表面之间形成的压力气膜将负荷支撑起来的轴承，工作时滑动副表面之间完全由气膜分开。由于两个相对运动部件之间有一层气膜，其中一个处于“悬浮”状态，因此也被称为气浮轴承。与滚动轴承和油滑动轴承相比，气体轴承利用气体作为润滑剂，摩擦力矩和摩擦系数很小，对使用环境和使用部位也没有任何污染。同时，具有速度高、精度高、功耗低、寿命长等优点，是一种理想的支撑元件。

目前针对动压气体轴承的模拟研究基本都是采用基于连续性假设得出的宏观雷诺方程以及其修正形式，但是由于转静子间隙通常很小，间隙内流动状态如果一旦属于过渡区，那么将不再满足连续性假设，这时雷诺方程也不再适用。因此需要采用更加准确合适的数值方法，格子Boltzmann方法(LBM)可以模拟宏观尺度到介观尺度，已被证明可以模拟较大Kn数的微流动，有望成为一种替代方法。格子Boltzmann方法不但可以看作由离散的格子气自动机方法的发展而来，而且也可以认为是连续玻尔兹曼方程的一种特殊的有限差分格式。具有程序易于实施，并行性能好，边界条件处理简单等优点。

发明内容

针对现有技术的缺陷，本发明提供了一种基于格子Boltzmann动压气体轴承间隙微流动仿真方法，解决了转静子间隙内流态模拟不准确的技术问题。

本发明采用的技术方案如下：

一种基于格子Boltzmann动压气体轴承间隙微流动仿真方法，选取合适的多松弛时间格子Boltzmann模型，然后将动压气体轴承简化为刚性轴承，在其转静子间隙流场内生成贴体网格，并采用广义插值格子玻尔兹曼方法进行求解。为了实现微流动的模拟，将松弛时间与Kn数进行关联，同时基于非平衡外推格式，将Langmuir滑移模型引入格子Boltzmann方法中，以满足间隙微流动在滑移区和过渡区的模拟；最后可以得到流场内的速度分布情况以及气膜压力，从流场特征角度分析轴承机理，为更好的设计轴承提供一定的理论基础。

具体包括以下步骤：

步骤一：建立轴承转子与静子之间间隙流场的物理模型，根据轴颈半径和以及和轴承套之间的间隙尺寸，基于求解椭圆型偏微分方程生成自适用的贴体网格；

步骤二：采用多松弛时间格子Boltzmann模型，根据维度和离散速度数量，选择对应的演化方程以及秦勒展开形式的麦克斯韦平衡态分布函数；

步骤三：根据贴体网格物理平面和欧拉网格计算平面的转换关系，求解计算平面内各格点的离散速度，得出插值系数；

步骤四：设置和计算初始值，根据大气温度，压力，并进行无量纲转换，得到流体初始速度和密度，以及克努森数Kn、雷诺数Re、马赫数Ma、松弛因子τ_s等，然后初始化速度分布函数；

步骤五：根据演化方程得到贴体网格各格点的碰撞后分布函数f_i ^*(ξ，t)，

t、ξ分别代表时间和计算平面中的位置，结合步骤三得出的所述插值系数，计算碰撞后(ξ₁(x-c_tΔt)，ξ₂(x-c_iΔt))位置的粒子分布函数，最终得到下一时间步长的粒子分布函数f_i(ξ，t+Δt)；

步骤六：边界条件的处理：当转子与静子之间间隙内气体流型属于连续区时，采用非平衡外推方法：

f_i(x_b，t)＝f_i ^(eq)(ρ(x_b，t)，u_w)+[f_i(x_f，t)-f_i ^(eq)(x_f，t)]

式中：x_b是边界格点，x_f是与边界法向上最近的流体格点，ρ(x_b，t)是边界处的密度，u_w是边界的速度，f_i ^eq是平衡态分布函数；

当所述气体流型演变成滑移流或者过渡流时，采用滑移边界条件，当Kn≤0.1，即处于滑移区时，采用一阶滑移边界条件：

f_i(x_b，t)＝α(f_i ^eq(x_b，t)-f_i ^eq(x_f，t))+f_i(x_f，t)

当Kn≥0.1，即处于过渡区时，采用二阶滑移边界条件：

上两式中，α、β为滑移系数，定义如下：

其中，A₁和A₂取决于气、壁相互作用的参数，N＝1/δ_x为格子数，x_ff为边界法向方向上第二近的流体格点；

步骤七：计算流体的宏观速度和密度，并根据间隙流场内各处的密度修正Kn数；

步骤八：重复步骤五至七，直至满足收敛条件。

其进一步技术方案为：

所述步骤一中，所述椭圆型偏微分方程如下：

式中，ξ和η为计算平面的变量，x和y为物理平面的变量。

所述步骤二中，在时间t和位置x处，具有离散速度c_i的速度分布函数f_i(x，t)的演化表示为：

f_i(x+c_iδ_t，t+δ_t)-f_i(x，t)＝Ω_i(f)+δ_tF_i

其中，Ω_i(f)为离散碰撞算子，F_i是离散外力项，δ_t是时间步长；采用多重松弛时间的碰撞算子，其定义如下：

其中，M为速度空间到矩空间的正交变换矩阵，S为矩空间的对角松弛矩阵，

为平衡态分布函数：

式中，ω_j为权系数，c_s为声速，ρ为格子密度，u为格子速度。

所述步骤三中，求解计算平面内各格点的离散速度

然后对

在时间步长Δt上采用二步龙格库塔法计算积分，根据积分结果得出所述插值系数。

所述步骤三中，对于流场内部点采用的是二阶迎风插值法，而在与边界相邻的格点上采用二阶中心插值法。

所述步骤五中，计算碰撞后的粒子分布函数f_i ^*(ξ，t)的方程式为：

所述步骤八中，所述收敛条件为：

式中，u为格子速度。

所述步骤七中，随着间隙通道内各处压力的变化，Kn数也会随着发生改变，进而影响到松弛因子τ_s的数值，采用下式进行修正：

式中，ρ为格子密度。

本发明的有益效果如下：

本发明提出了采用介观的格子Boltzmann方法，模拟研究动压气体轴承转静子间隙微通道内的气体流动特性，更好的体现滑移区和过渡区的微尺度效应，可以捕捉轴颈和轴承套两侧的速度滑移，提高了准确度。

附图说明

图1是本发明的介观模型的求解流程图。

图2是旋转Couette流的贴体网格划分示意图。

图3是本发明方法求解的不同Kn数情况下理论解和数值解的切向速度分布。

图4是本发明方法求解的不同偏心率下间隙内圆周方向的无量纲压力分布。

图5是本发明方法求解的轴承轴颈侧切向速度沿圆周方向的分布。

图6是本发明方法求解的不同转速下间隙内圆周方向的无量纲压力分布。

具体实施方式

以下结合附图说明本发明的具体实施方式。

本实施例的一种基于格子Boltzmann动压气体轴承间隙微流动仿真方法，如图1所示，包括以下步骤：

步骤一：

建立轴承转子与静子之间间隙流场的物理模型，传统的弹性箔片动压气体轴承间隙内会存在可变性的波箔，为了方便模拟，此处忽略波箔的影响，将动压气体轴承简化为刚性轴承，根据轴颈和轴承套的尺寸参数，基于椭圆型偏微分方程生成贴体网格；

具体地，已知的轴颈半径为10.95mm、轴承套的半径为11mm、偏心率为0.2～0.8，然后基于椭圆型偏微分方程生成转静子间隙流场的贴体网格，最终网格数为27520×20，可采用最简单的椭圆型方程Laplace方程：

其中，ξ和η为计算平面的变量，x和y为物理平面的变量。

步骤二：采用多松弛时间格子Boltzmann模型，根据维度和离散速度数量，选择对应的演化方程以及秦勒展开形式的麦克斯韦平衡态分布函数；

在不失一般性的情况下，在时间t和位置x处，具有离散速度c_i(i代表离散速度方向)的速度分布函数f_i(x，t)的演化可以表示为：

f_i(x+c_iδ_t，t+δ_t)-f_i(x，t)＝Ω_i(f)+δ_tF_i

其中，Ω_i(f)为离散碰撞算子，F_i是离散外力项，δ_t是时间步长。为了更好的模拟微流动，这里采用多重松弛时间(MRT)的碰撞算子，其定义如下：

对于二维情况，采用标准的D2Q9模型，M为速度空间到矩空间的正交变换矩阵，具体取值如下：

S为矩空间的对角松弛矩阵：

式中，ω_j为权系数，取值为ω₀＝4/9，ω₁＝ω₂＝ω₃＝ω₄＝1/9，ω₅＝ω₆＝ω₇＝ω₈＝1/36，

为声速，ρ为格子密度，u为格子速度，离散速度ci为：

步骤三：根据贴体网格物理平面和欧拉网格计算平面的转换关系，求解计算平面内各格点的离散速度，得出插值系数；

计算平面内各格点的离散速度

定义如下：

对于上式，可利用计算平面和物理平面的如下关系计算：

其中，J为雅克比行列式：

时间步长Δt必须满足CFL数：

式中，λ为CFL数(0＜λ≤1)，这里取λ＝1。

然后对

在时间步长Δt上进行积分，为了满足LBM的二阶数值精度，采用二步龙格库塔法计算积分，其积分的具体步骤如下：

：第一步

：第二步

根据上式计算插值系数。

步骤四：

设置和计算初始值，根据大气温度，压力，并进行无量纲转换，得到流体初始速度和密度(压力)，以及克努森数Kn、雷诺数Re、马赫数Ma、松弛因子τ_s等，然后初始化速度分布函数；

克努森数Kn的变化范围是0.018～0.25、转轴转速取30000r/min～70000r/min，流体初始速度为0，密度(压力)为1、松弛因子τ_s根据Kn数计算，然后初始化速度分布函数。

需要注意的一点是，在偏心的情况下，转静子间隙各处的尺寸不同，所以Kn数也会不同，因此在这里基用气膜厚度的公式进行计算。假设同心情况下已知Kn为Kn₀，那么在偏心率为ε的时候，从间隙最大处到最小处，不同角度θ处的Kn数有：

Kn(θ)＝Kn₀/(1+εcos(θ))

步骤五：

根据演化方程得到贴体网格各格点的碰撞后分布函数f_i ^*(ξ，t)，t、ξ分别代表时间和计算平面中的位置；

结合步骤三得出的所述插值系数，计算碰撞后(ξ₁(x-c_tΔt)，ξ₂(x-c_iΔt))位置的粒子分布函数，最终得到下一时间步长的粒子分布函数f_i(ξ，t+Δt)。

步骤六：

边界条件的处理：

当转子与静子之间间隙内气体流型属于连续区时，采用非平衡外推方法：

f_i(x_b，t)＝f_i ^(eq)(ρ(x_b，t)，u_w)+[f_i(x_f，t)-f_i ^(eq)(x_f，t)]

式中：x_b是边界格点，x_f是与边界法向上最近的流体格点，ρ(x_b，t)是边界处的密度，u_w是边界的速度，f_i ^eq是平衡态分布函数；

当气体流型演变成滑移流或者过渡流时，采用滑移边界条件，当Kn≤0.1，即处于滑移区时，采用一阶滑移边界条件：

f_i(x_b，t)＝α(f_i ^eq(x_b，t)-f_i ^eq(x_f，t))+f_i(x_f，t)

当Kn≥0.1，即处于过渡区时，采用二阶滑移边界条件：

上两式中，α、β为滑移系数，定义如下：

其中，A₁和A₂取决于气、壁相互作用的参数，N＝1/δ_x为格子数；x_ff为边界法向方向上第二近的流体格点；

步骤七：计算流体的宏观速度和密度，并根据间隙流场内各处的密度修正Kn数；

计算宏观速度和密度，计算式如下：

式中：ρ为格子密度、u_α为宏观速度。

随着间隙通道内各处压力的变化，Kn数也会随着发生改变，进而影响到松弛因子τ_s的数值，因此根据密度采用下式修正Kn数：

ρ为格子密度。

步骤八：重复步骤五至七，直至满足一定的收敛条件：

最后输出速度、压力等数据。

首先可利用两个圆柱间的旋转微Couette流验证新建立模型的准确性和有效性，贴体网格的生成情况如图2所示，网格数为188×20。在模拟中，半径为R1的内圆柱以恒定的角速度ω顺时针旋转，而半径为R2的外圆柱保持静止，其中R1∶R2＝3∶5。该流动的无量纲切向速度理论解为：

其中，

σ₁和σ₂分别为内外圆柱的切向动量调节系数，λ是平均分子自由程，r是径向距离。图3展示了Kn＝0.05和Kn＝0.1两种情况下的切向速度分布，从图中可以看出，不论在滑移区还是过渡区，数值解都和理论解有着较好的一致性，证明了本实施模型在微流动中的准确性。

上述实施例的步骤三中，对于流场内部点采用的是二阶迎风插值法，而在与边界相邻的格点上采用二阶中心插值法。插值函数为：

其中，

md＝sgn(Δξ_up，i，ξ)，nd＝sgn(Δξ_up，i，η)

分别表示插值过程中插值点的取值方向，ξ_m，η_n表示插值点的位置。

a)对于二阶迎风插值，插值系数为：

a_i，1，α＝-|Δξ_up，i，α|(|Δξ_up，i，α|-2)

b)对于二阶中心插值，插值系数为：

上述实施例的步骤五中，克努森数Kn、雷诺数Re、马赫数Ma的关系式如下：

松弛因子τ_s需要与Kn数建立联系，以模拟微流动，其关系式如下。

式中：H为特征长度，δ_x为格子步长。

在基于模型验证准确的前提下，对于上述实施例的动压气体轴承进行模拟计算，控制其余参数一致，仅改变偏心率这一变量，计算得到了不同偏心率下动压气体轴承间隙内圆周方向的无量纲压力分布情况，如图4所示，当偏心率从0.2增加到0.8时，气膜压力的极值逐渐提高，且极大值和极小值的所处的位置均慢慢向间隙最小处(180°)靠近，此外波峰和波谷的形状也变得越来越尖锐，说明压力极大值和极小值的分布更加集中。这些结论与宏观雷诺方程得到的结论一致，再一次说明了本模型的有效性。

图5是以间隙最大处为起点，沿圆周方向遍布一周时，轴颈侧表面气体无量纲速度的分布情况。可以看出，随着间隙尺寸的减小，轴颈侧表面气体速度也逐渐减小，表明发生的滑移量逐渐增大，并在最小间隙处，达到最大滑移量。随着偏心率的增加，表面各处的滑移速度都有着不同程度的增大，在最小间隙处，滑移量的增大最为明显。

在偏心率为0.4以及其余参数一致的前提下，仅改变轴承的旋转速度，可以模拟得到在不同转速下气体轴承间隙圆周方向的气膜压力分布情况，如图6所示。随着转速的提高，气膜压力的极值也均有所提高，不过提升的速率逐渐减慢。压力极大值的位置基本不发生变化，但极小值的位置逐渐向远离间隙最小处(180°)的位置移动。结合偏心率变化规律，说明在研究范围内，偏心率增大和转速的提高均有益于轴承性能的提升。

Claims (8)

1.一种基于格子Boltzmann动压气体轴承间隙微流动仿真方法，其特征在于，包括以下步骤：

步骤一：建立轴承转子与静子之间间隙流场的物理模型，根据轴颈半径和以及和轴承套之间的间隙尺寸，基于求解椭圆型偏微分方程生成自适用的贴体网格；

步骤二：采用多松弛时间格子Boltzmann模型，根据维度和离散速度数量，选择对应的演化方程以及秦勒展开形式的麦克斯韦平衡态分布函数；

步骤三：根据贴体网格物理平面和欧拉网格计算平面的转换关系，求解计算平面内各格点的离散速度，得出插值系数；

步骤四：设置和计算初始值，根据大气温度，压力，并进行无量纲转换，得到流体初始速度和密度，以及克努森数Kn、雷诺数Re、马赫数Ma、松弛因子τ_s等，然后初始化速度分布函数；

步骤五：根据演化方程得到贴体网格各格点的碰撞后分布函数f_i ^*(ξ，t)，

t、ξ分别代表时间和计算平面中的位置，结合步骤三得出的所述插值系数，计算碰撞后(ξ₁(x-c_tΔt)，ξ₂(x-c_iΔt))位置的粒子分布函数，最终得到下一时间步长的粒子分布函数f_i(ξ，t+Δt)；

步骤六：边界条件的处理：当转子与静子之间间隙内气体流型属于连续区时，采用非平衡外推方法：

f_i(x_b，t)＝f_i ^(eq)(ρ(x_b，t)，u_w)+[f_i(x_f，t)-f_i ^(eq)(x_f，t)]

式中：x_b是边界格点，x_f是与边界法向上最近的流体格点，ρ(x_b，t)是边界处的密度，u_w是边界的速度，f_i ^eq是平衡态分布函数；

当所述气体流型演变成滑移流或者过渡流时，采用滑移边界条件，当Kn≤0.1，即处于滑移区时，采用一阶滑移边界条件：

f_i(x_b，t)＝α(f_i ^eq(x_b，t)-f_i ^eq(x_f，t))+f_i(x_f，t)

当Kn≥0.1，即处于过渡区时，采用二阶滑移边界条件：

上两式中，α、β为滑移系数，定义如下：

β＝-A₂(KnN)²α

其中，A₁和A₂取决于气、壁相互作用的参数，N＝1/δ_x为格子数，x_ff为边界法向方向上第二近的流体格点；

步骤七：计算流体的宏观速度和密度，并根据间隙流场内各处的密度修正Kn数；

步骤八：重复步骤五至七，直至满足收敛条件。

2.根据权利要求1所述的基于格子Boltzmann动压气体轴承间隙微流动仿真方法，其特征在于，所述步骤一中，所述椭圆型偏微分方程如下：

式中，ξ和η为计算平面的变量，x和y为物理平面的变量。

3.根据权利要求2所述的基于格子Boltzmann动压气体轴承间隙微流动仿真方法，其特征在于，所述步骤二中，在时间t和位置x处，具有离散速度c_i的速度分布函数f_i(x，t)的演化表示为：

f_i(x+c_iδ_t，t+δ_t)-f_i(x，t)＝Ω_i(f)+δ_tF_i

其中，Ω_i(f)为离散碰撞算子，F_i是离散外力项，δ_t是时间步长；采用多重松弛时间的碰撞算子，其定义如下：

其中，M为速度空间到矩空间的正交变换矩阵，S为矩空间的对角松弛矩阵，

为平衡态分布函数：

式中，ω_j为权系数，c_s为声速，ρ为格子密度，u为格子速度。

4.根据权利要求3所述的基于格子Boltzmann动压气体轴承间隙微流动仿真方法，其特征在于，所述步骤三中，求解计算平面内各格点的离散速度

然后对

在时间步长Δt上采用二步龙格库塔法计算积分，根据积分结果得出所述插值系数。

5.根据权利要求3所述的基于格子Boltzmann动压气体轴承间隙微流动仿真方法，其特征在于，所述步骤三中，对于流场内部点采用的是二阶迎风插值法，而在与边界相邻的格点上采用二阶中心插值法。

6.根据权利要求3所述的基于格子Boltzmann动压气体轴承间隙微流动仿真方法，其特征在于，所述步骤五中，计算碰撞后的粒子分布函数f_i ^*(ξ，t)的方程式为：

。

7.根据权利要求6所述的基于格子Boltzmann动压气体轴承间隙微流动仿真方法，其特征在于，所述步骤八中，所述收敛条件为：

式中，u为格子速度。

8.根据权利要求1所述的基于格子Boltzmann动压气体轴承间隙微流动仿真方法，其特征在于，所述步骤七中，随着间隙通道内各处压力的变化，Kn数也会随着发生改变，进而影响到松弛因子τ_s的数值，采用下式进行修正：

式中，ρ为格子密度。

Images