包含完整变量扰动的流体润滑频变动力学特性计算方法与流程

1.本发明属于高速无油轴承技术领域，特别涉及一种包含完整变量扰动的流体润滑频变动力学特性计算方法。


背景技术：

2.在能源动力、国防军工、航空航天、车辆舰船等领域中，高速旋转机械得到了广泛的应用，随着对其性能和效率的要求提高，高速轴承和密封的性能分析和设计，以及稳定性的控制已成为关键技术。工质流体自润滑轴承(如混合气体、液氧、超临界二氧化碳等)在彻底避免工质被油污染的同时，还能降低密封技术要求从而缩短转子跨距以提高其刚性，成为国内外旋转机械领域的一个重要方向。
3.然而特种参数下的流体通常具有复杂的热物性，并由此带来润滑膜流动的多种附加效应，比如真实气体效应、湍流与惯性流效应、稀薄效应等。具体到数学模型则体现为润滑膜的控制方程中不只含有膜厚和压力两个显函数变量，还存在与压力、温度关系复杂而无法代入表示的热物性变量，以及由热物性和膜厚的非线性函数表示的附加效应修正。
4.尽管国内外已有很多研发人员对具有附加效应的润滑问题进行研究，但绝大多数只限于轴承或干气密封的稳态特性。然而要对转子系统的临界特性和稳定性进行分析，离不开动态刚度和阻尼系数表示的动力学特性。而且无论是可压缩流体润滑，还是润滑膜工作面存在动态变化(比如可倾瓦轴承的瓦块摆动、箔片轴承的变形)，动力学刚度和阻尼系数都具有频率效应。现有的动特性求解方法要么对控制方程进行准静态处理，只能得到扰动频率无限接近于零的静刚度，要么在频率扰动处理过程中只考虑了部分附加效应的动态变化。
5.因此，需提供一种包含完整变量扰动的偏导数法来解决上述问题。


技术实现要素：

6.为了解决上述问题，本发明目的在于提供一种用包含完整变量扰动的偏导数法，主要目的在于计算任意扰动频率下流体润滑轴承的动力学特性，所述方法用于求取任意扰动频率下流体润滑轴承动力学系数，包括如下步骤：
7.s1：结合润滑介质的薄膜流动特性建立包含附加效应的瞬态雷诺方程，并对其进行频率扰动分析，推导瞬态雷诺方程中出现的完整变量扰动展开表达式，进而得到稳态雷诺方程和包含完整变量动态变化的扰动雷诺方程；其中，附加效应包括热物性变化、湍流效应、惯性流效应、稀薄效应等，所述完整变量扰动展开表达式不只是膜厚和压力扰动，还包括由膜厚和压力动态映射关系构造的热物性扰动和扰动形式的附加效应函数；
8.s2：针对具体轴承结构形式建立动力学模型，并对动力学模型进行结构扰动，确定s1中的扰动膜厚与轴承力扰动的关系，获得隐含在转子扰动和结构扰动中的扰动变量所对应的复数压力分布控制方程；
9.s3：数值求解s1中的稳态雷诺方程和s2中的扰动变量对应的复数压力分布控制方
程，得到扰动变量对应的复数压力分布；
10.s4：利用s3中获得的扰动变量对应的复数压力分布求取动力学刚度和阻尼系数。
11.本发明所提供的包含完整变量扰动的流体润滑频变动力学特性计算方法，还具有这种的特征，包括如下步骤，所述s1包括：
12.s1.1：建立润滑介质薄膜流动的瞬态雷诺方程，包含由热物性变量和膜厚的非线性函数表示的附加效应；
13.s1.2：结合润滑薄膜流动特性选取特征参数对瞬态雷诺方程进行无量纲化；
14.s1.3：转子系统临界特性与稳定性分析对应的状态为，转子在旋转的同时，轴颈中心以任意扰动频率绕着平衡位置在小邻域内正向涡动，处于扰动状态下的润滑膜通过含扰动频率与旋转频率比值的复指数表示压力和膜厚的动态变化，得到频率扰动展开形式的无量纲瞬态膜厚和压力；
15.具体为：转子以频率ω旋转，轴颈中心绕着平衡位置以频率υ在小邻域内涡动，处于扰动状态下的润滑膜通过含扰动频率与旋转频率比值的复指数表示压力和膜厚的动态变化，即：
[0016][0017]
其中为频率扰动展开形式的无量纲瞬态膜厚和压力，为无量纲稳态膜厚和压力，扰动膜厚和扰动压力分别为和复指数中的涡动比ω＝υ/ω，无量纲时间
[0018]
s1.4：将控制方程中的热物性变量统一进行一阶频率扰动展开，并通过热物性对压力的无量纲偏导数构造热物性变量扰动与扰动压力之间的动态映射，得到无量纲方程中含油的密度和黏度扰动展开形式；
[0019]
具体为：将控制方程中的热物性变量统一表示为并进行一阶频率扰动展开并通过如下表达式构造热物性变量扰动与扰动压力之间的动态映射。
[0020][0021]
其中，为热物性对压力的无量纲偏导数。
[0022]
s1.5：将扰动膜厚和热物性变量扰动代入附加效应的非线性函数表达式，借助等价无穷小理论，推导由稳态附加效应函数与扰动膜厚和热物性变量扰动线性组合叠加表示的附加效应动态变化；
[0023]
s1.6：将s1.3至s1.5得到的所有变量的频率扰动展开形式代入s1.2得到的无量纲瞬态雷诺方程，利用方程给定的等式关系分离零阶和一阶项，分别得到的稳态雷诺方程和的扰动雷诺方程。
[0024]
本发明所提供的包含完整变量扰动的流体润滑频变动力学特性计算方法，还具有这种的特征，包括如下步骤，所述s2包括：
[0025]
s2.1：针对具体轴承结构形式建立动力学模型，将所述动力学模型无量纲化，获得无量纲轴承动力学模型；
[0026]
s2.2：对s2.1获得的无量纲轴承动力学模型进行结构扰动，确定扰动膜厚与轴承力扰动的关系；
[0027]
s2.3：将s1中获得的扰动雷诺方程对隐含在转子扰动和结构扰动中的扰动变量求偏导数，并借助s2.2中扰动膜厚与轴承力扰动的关系，获得各扰动变量所对应的复数压力分布控制方程。
[0028]
本发明所提供的包含完整变量扰动的流体润滑频变动力学特性计算方法，还具有这种的特征，包括如下步骤，所述s3包括：
[0029]
s3.1：将s1中获得的稳态雷诺方程离散，采用超松弛迭代法求解稳态压力分布然后基于稳态压力分布利用润滑介质的热物性模型或数据库更新稳态热物性分布进而基于无量纲稳态膜厚和稳态热物性分布更新稳态附加效应函数；重复上述过程直至稳态压力分布收敛；
[0030]
s3.2：采用数值微分方法对求解域内的网格节点在稳态压力分布处计算热物性对压力的无量纲偏导数
[0031]
s3.3：将s3.1和s3.2得到的所有稳态变量分布作为输入条件，用数值方法求解s2所获得的各扰动变量对应的复数压力分布控制方程。
[0032]
本发明所提供的包含完整变量扰动的流体润滑频变动力学特性计算方法，还具有这种的特征，包括如下步骤，所述s4包括：
[0033]
基于所选取的坐标系，对步骤s3获得的各扰动变量对应的复数压力分布进行数值积分，再利用旋转矩阵变换得到该坐标系下的动力学刚度和阻尼系数。
[0034]
有益效果：
[0035]
本发明所提供的包含完整变量扰动的流体润滑频变动力学特性计算方法，从具有附加效应的薄膜润滑雷诺方程出发，通过建立动态映射给出一种不局限于特定润滑流体或热物性模型的静动特性分析通用方法—完整变量扰动的偏导数法。该方法能够统一处理结构扰动和润滑膜附加效应的动态变化，适用于求解轴承和干气密封在不同扰动频率下的动力学刚度和阻尼系数。
附图说明
[0036]
为了更清楚地说明本发明的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0037]
图1为本发明实施例所提供的包含完整变量扰动的流体润滑频变动力学特性计算方法的流程图；
[0038]
图2为本发明实施例所提供的包含完整变量扰动的流体润滑频变动力学特性计算方法中扰动状态下超临界二氧化碳润滑整周式径向箔片轴承的示意图；
[0039]
图3为本发明实施例所提供的包含完整变量扰动的流体润滑频变动力学特性计算方法中热物性变量扰动与扰动压力之间的动态映射示意图。
[0040]
其中：1：顶层平箔；2：支撑拱箔；3：润滑膜；4：扰动状态下轴颈瞬时位置；5：轴颈平衡位置。
具体实施方式
[0041]
为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，以下实施例结合附图对本发明所提供的包含完整变量扰动的流体润滑频变动力学特性计算方法作具体阐述。
[0042]
超临界二氧化碳(s
‑
co2)作为已有应用的特种参数润滑介质中热物性和润滑薄膜流动特性最复杂的之一，以s
‑
co2工质自润滑箔片轴承为例，在轴颈扰动状态下整周式径向箔片轴承的示意图如图2所示。
[0043]
计算时需要事先给定轴承的结构和运行参数，包括轴承间隙c、轴颈直径d、轴承长l、箔片无量纲柔度环境压力p
a
、环境温度t
a
、偏心率ε0、转速，涡动比ω
[0044]
如图1所示，为包含完整变量扰动的流体润滑频变动力学特性计算方法的流程图，具体实施步骤如下：
[0045]
1)结合润滑介质的薄膜流动特性建立包含附加效应的瞬态雷诺方程，并对其进行频率扰动分析，推导稳态雷诺方程和包含完整变量动态变化的扰动雷诺方程。
[0046]
1.1)建立考虑附加效应的超临界二氧化碳薄膜流动的瞬态雷诺方程，通用形式的可压缩湍流润滑雷诺方程如下：
[0047][0048]
其中，r为轴颈半径，环向角度坐标θ如图2所示，p为润滑膜压力，ρ和μ分别为密度和黏度，k
θ
和k
z
描述湍流效应，由热物性变量ρ、μ和膜厚h的非线性函数表示。
[0049][0050][0051]
上式中局部雷诺数re＝ρωrh/μ，常数α1，β1，α2，β2的值由湍流润滑模型确定；上式以ng
‑
pan模型为例，当re<1000时上述所有常数均等于零。
[0052]
1.2)结合润滑薄膜流动特性选取适当的特征参数对瞬态雷诺方程进行无量纲化。则引入以下变换，无量纲轴向坐标λ＝z/r；无量纲压力、密度和黏度下标a表示环境参数；无量纲膜厚以及由旋转频率ω定义的无量纲时间代入步骤1.1)的控制方程，可得无量纲瞬态雷诺方程如下：
[0053][0054]
其中λ＝μ
a
ω/2p
a
ψ2是轴承数，ψ＝c/r为间隙比。
[0055]
1.3)转子系统临界特性与稳定性分析对应的状态为，转子以频率ω旋转，轴颈中心绕着平衡位置以频率υ在小邻域内涡动，处于扰动状态下的润滑膜通过含扰动频率的复指数表示压力和膜厚的动态变化，即：
[0056][0057]
上式中分别为频率扰动展开形式的无量纲瞬态膜厚和压力，为无量纲稳态膜厚和压力，扰动膜厚和扰动压力分别为和复指数中的涡动比ω＝υ/ω，无量纲时间如步骤1.2)所述。
[0058]
1.4)将控制方程中的热物性变量统一表示为并进行一阶频率扰动展开并通过如下表达式构造热物性变量扰动与扰动压力之间的动态映射，如附图3所示。
[0059][0060]
其中，为热物性对压力的无量纲偏导数。
[0061]
具体到步骤1.2)的无量纲方程中含有的密度和黏度扰动展开如下：
[0062][0063][0064]
1.5)将扰动膜厚和热物性变量扰动代入附加效应的非线性函数表达式，借助等价无穷小理论，推导由稳态附加效应函数与线性组合叠加表示的附加效应动态变化。
[0065]
具体到步骤1.2)的无量纲方程，应将和代入湍流系数的表达式。
[0066]
步骤1.1)湍流系数k
θ
、k
z
的公式中局部雷诺数re的动态变化反映了湍流对动力学特性的影响，引入由轴承间隙、环境密度和环境黏度定义的平均雷诺数re
a
＝ρ
a
ωrc/μ
a
，可得由无量纲变量表示的动态局部雷诺数。
[0067][0068]
由幂函数的泰勒公式，结合步骤1.4)的黏度扰动展开式，有
[0069][0070]
将上式以及步骤1.3)的膜厚扰动和步骤1.4)的密度扰动代入动态局部雷诺数，略去二阶及以上的扰动项，得到动态局部雷诺数re的扰动形式。
[0071][0072]
其中是对应于轴颈平衡位置的稳态雷诺数，由以下线性组合构成的扰动雷诺数是扰动密度、扰动黏度和扰动膜厚的联合响应。
[0073]
借助幂函数的泰勒公式不难推出re
β
。
[0074][0075]
雷诺方程中的湍流系数存在于分母，使得方程表现出更明显的非线性。将两个方向的湍流系数统一写为k＝12 αre
β
，将re
β
代入湍流系数k的表达式同时在其分母中加入一个可作为小量而存在的项并定义与轴颈平衡位置对应的稳态湍流系数k0＝12 αre
0β
，可得湍流系数的复扰动展开：
[0076][0077]
轴颈扰动下湍流系数的动态变化通过扰动雷诺数re
d
起作用。可以看出，从和到扰动湍流系数k
θ
，k
z
是复合动态映射。
[0078]
1.6)将步骤1.3)至步骤1.5)中得到的所有变量的频率扰动展开形式代入步骤1.2)的无量纲瞬态雷诺方程，利用方程给定的等式关系分离零阶和一阶项，分别得到的稳态雷诺方程和的扰动雷诺方程。
[0079][0080][0081]
稳态雷诺方程和扰动雷诺方程的边界条件分别为：
[0082]
[0083][0084]
对于箔片轴承，由于顶层平箔不能承受小于环境压力pa的润滑膜压力，还应补充雷诺边界条件如下：
[0085][0086]
2)针对具体轴承结构形式建立动力学模型，并对其进行结构扰动，确定步骤1)中的扰动膜厚与轴承力扰动的关系。获得隐含在转子扰动和结构扰动中的各扰动变量所对应的复数压力分布控制方程。
[0087]
2.1)对于整周式径向箔片轴承，顶层平箔与轴颈表面形成的间隙为气膜厚度，动态变形的底层拱箔可建模为阻尼柔性支承。考虑柔度和等效阻尼的综合箔片模型如下：
[0088][0089]
其中w
t
是平箔位移，以轴承间隙c为特征量对其无量纲化。将无量纲压力和无量纲时间代入上式，可得无量纲箔片动力学模型如下：
[0090][0091]
其中无量纲柔度和无量纲等效阻尼与对应的有量纲形式有如下转换。
[0092][0093][0094]
在计算中给定箔片的无量纲柔度即可得到轴承性能，因而本发明不局限于确定柔度α的具体箔片模型。以临界阻尼比作为箔片结构的等效阻尼。
[0095]
2.2)对步骤2.1)的无量纲轴承动力学模型进行结构扰动，确定扰动膜厚与轴承力扰动的关系。平箔位移的动态变化展开为频率扰动形式如下：
[0096][0097]
将上式与步骤1.3)的压力扰动展开代入步骤2.1)的无量纲箔片动力学模型对其进行结构扰动，利用等式关系分离零阶和一阶项，得到平箔稳态位移和扰动位移如下：
[0098][0099][0100]
进而可得步骤1.3)中平衡位置稳态膜厚和扰动膜厚的具体表达式。
[0101][0102]
2.3)将步骤1.6)中得到的扰动雷诺方程对隐含在转子扰动和结构扰动中的各扰动变量求偏导数，并借助步骤2.2)中扰动膜厚与轴承力扰动的关系，获得各扰动变量所对应的复数压力分布控制方程。
[0103]
步骤2.2)中扰动膜厚含有和两个扰动变量，将其分别对和求偏导数。
[0104][0105]
其中，对应的扰动压力偏导数为：
[0106][0107]
结合以上两式，对步骤1.6)的扰动雷诺方程求偏导数，得到和两个扰动变量对应的复数压力分布p
ε
和的控制方程，如下：
[0108][0109]
[0110]
以上两个方程的边界条件与步骤1.6)的扰动雷诺方程一致，将分别替换为p
ε
和即可。
[0111]
3)数值求解步骤1)中的稳态雷诺方程和步骤2)中各扰动变量对应的复数压力分布控制方程。
[0112]
3.1)以五点差分格式离散步骤1.6)得到的稳态雷诺方程，采用超松弛迭代法求解稳态压力分布然后基于更新超临界二氧化碳润滑膜的稳态密度分布和黏度分布进而基于和更新稳态湍流系数k
θ0
和k
z0
分布；多次重复以上过程直到稳态压力分布收敛。
[0113]
3.2)采用数值微分方法对求解域内的网格节点在处计算密度和黏度对压力的无量纲偏导数和
[0114]
3.3)将步骤3.1)和步骤3.2)得到的所有稳态变量分布作为输入条件，用有限差分法求解步骤2.3)中复数压力分布p
ε
和的控制方程。
[0115]
4)基于图2所示的x
‑
y坐标系，通过以下公式对步骤3)求得的复数压力分布p
ε
和进行数值积分。
[0116][0117]
再利用旋转矩阵[a]变换得到该坐标系下的无量纲动力学刚度和阻尼系数。
[0118][0119][0120]
动力学刚度和阻尼系数与对应无量纲形式的转换如下：
[0121][0122]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。比如从包含完整变量动态变化的扰动雷诺方程中去掉对结果影响不明显的扰动项。以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变型，这些改进和变型也应视为本发明的
保护范围。