周如能，顾永鹏，任革学，于溯源
(清华大学 a.能源与动力工程系；
b.核能与新能源技术研究院；
c.航天航空学院，北京 100084)
箔片动压气体轴承具有无油、高温和高速性能好等优点，被广泛应用于飞机的空气循环机和燃料电池的压缩机等领域[1-3]，其中波箔式箔片轴承应用最为广泛[4-5]。由于波箔结构的特殊性，波箔动压气体轴承系统一部分非线性来源于箔片结构的接触摩擦行为[6]，系统动力学建模中箔片结构的接触摩擦模型是准确预测波箔动压气体轴承静、动态特性的关键[7]。
文献[8]最早提出了波箔气体轴承的理论模型——简单弹性基础模型，后来被拓展到动力学计算中，并引入了黏性阻尼近似模拟结构耗散特性[9-10]，但结构耗散特性与轴承-转子系统的运动状态关系较大，根据某一工况估算的等效黏性阻尼不能准确表征实际的能量耗散机制。这一局限性使很多学者提出了考虑波箔结构及波箔之间相互作用的模型耦合摩擦模型，进而对波箔动压气体轴承系统进行静、动力学分析[11-13]。
摩擦模型较多，其选择受系统求解方法的限制。传统的显式分别迭代传参求解方法可以采用最经典的库仑摩擦模型，在每个时间步内迭代判断每个摩擦节点的状态并得到摩擦力[12]，但各子系统变量之间存在时间滞后，需极小的时间步长来保证收敛性，计算效率低，故文献[14-15]提出了全耦合求解方法，系统所有变量可同时求解，大幅提高了计算效率。但由于库仑摩擦模型在零速度时的不连续性而无法直接使用，需要改进摩擦模型。按是否采用微分方程库仑摩擦模型的改进可分为无需微分方程的代数摩擦模型和采用微分方程的动态摩擦模型[16]。最典型的代数摩擦模型是在经典库仑摩擦的间断处引入近似函数进行平滑处理，文献[17]采用代数摩擦模型实现了波箔动压气体轴承的全耦合建模与求解，但平滑模型无法准确描述静摩擦，即无法准确预测箔片节点的黏着状态，进而影响轴承性能的预测。波箔动压气体轴承分析中使用的动态摩擦模型主要包括Petrov和Ewins提出的动态摩擦模型[18]和LuGre动态摩擦模型[19-20]，前者同样是基于符号函数的近似衍生得到，理论上无法准确描述静摩擦，后者可以很好地描述斯特里贝克效应以及黏着-滑移现象。文献[19]最早将LuGre动态摩擦模型引入波箔动压气体轴承，并与link-spring模型进行联合仿真，但其分析局限于箔片内部的力学行为，未考虑转子和气膜的行为，且求解方法仍为分别迭代传参求解。
本文在上述研究的基础上，将LuGre动态摩擦模型引入全耦合的波箔动压气体轴承模型，提出了考虑转子、动压气膜、箔片结构、波箔之间的相互作用以及箔片节点之间黏着-滑移接触摩擦行为的综合动力学模型[20]。基于该模型，进一步对比了LuGre动态摩擦模型和平滑近似代数摩擦模型在程序实现、稳态特性以及轴承性能预测方面的差异。
波箔动压气体轴承结构如图1所示，以轴承中心O为原点建立坐标系Oxyz，θ为位置角(与x轴负方向的夹角)。工作时转子高速旋转，在转子与顶箔之间形成一层流体动压气膜，流体动压力作用于转子上使转子悬浮的同时也作用于顶箔上，而顶箔和波箔的受力变形会导致波箔与顶箔及轴承座之间出现复杂的接触摩擦行为。
1—气膜；
2—焊接端；
3—顶箔；
4—波箔；
5—轴承座。
为了求解这一包含可压缩流体润滑、接触摩擦等复杂非线性行为的流固耦合问题，气膜压力和箔片变形均采用标准的有限单元法进行空间离散，如图2所示。气膜压力沿气膜厚度方向的变化可忽略，不同位置的气膜压力可表示为p(θ,z)。对于第一代波箔动压气体轴承，一般假设箔片的曲率效应及变形沿轴承轴向的变化可忽略，为分析变形，以轴承周向为xf轴，径向为yf轴建立平面坐标系。为提高计算效率，箔片采用二维的梁单元建模。在标准的有限元框架下，气膜压力p(θ,z)和箔片变形u可以分别用单元节点向量插值得到，即
(1)
式中：Np，Nq分别为气膜压力和箔片变形的形函数矩阵；
pe和qe分别为单元节点的压力向量和变形向量。气膜压力和箔片变形的广义坐标向量分别标记为p和q。
1.1 转子的控制方程
因2套轴承对称支承一个刚性转子时，2套轴承位置转子响应相同，在此仅建立一套轴承模型，转子等效为质点，记xr为转子的位移矢量， 则转子
图2 有限单元法离散示意图
的控制方程为
(2)
式中：Mr为转子质量矩阵；
Wr为轴承静载荷，本文指转子重力；
Fr为气膜压力对转子的合力；
Fub为转子不平衡力。
Fr可通过气膜压力p(θ,z)积分得到，即
z∈[-L/2,L/2]，
(3)
Aθ=[cosθ,sinθ]T，
式中：R为轴承半径；
L为轴承长度；
pa为环境压力；
Ωe为单元面积。
不平衡力矢量Fub为
(4)
式中：u为不平衡量；
ω为转子角速度；
θ0为不平衡力的初始相位。
1.2 气膜的控制方程
气膜压力p(θ,z)的控制方程为等温可压缩理想流体的雷诺方程，即
(5)
U=[ωR/2,0]T，
式中：h为气膜厚度；
μ为空气动力黏度。
忽略转子和轴承的热膨胀以及离心膨胀，假设轴承间隙C为常数，气膜厚度h可表示为
h=C(1+excosθ+eysinθ)+hf，
(6)
[ex,ey]T=xr/C，
hf=-Nq,yfqe，
式中：ex，ey分别为在x，y方向量纲一的转子位移；
hf为箔片变形导致气膜厚度的变化量；
Nq,yf为顶箔单元对应于yf方向变形的形函数矩阵。
由于气膜与外界环境相通，气膜控制方程的边界条件为在所有边界上p=pa。
气膜压力采用双线性四节点单元进行划分，利用布勃诺夫-伽辽金法得到单元雷诺方程的弱形式为
(7)
Bp=∇Np，
气膜的整体控制方程可以通过组装所有单元方程得到
(8)
1.3 箔片的控制方程
采用二维欧拉梁单元建立箔片模型，单元的控制方程为
(9)
Ce=βKe，
式中：Me，Ce，Ke分别为单元的质量矩阵、阻尼矩阵(此处采用比例阻尼)和刚度矩阵；
Fe为单元的外部力矢量；
β为比例阻尼系数。此外，箔片材料弹性模量需根据平面应力假设修正。
由于波箔大部分节点不受力，且波箔动力学特性可以忽略，模型引入Guyan缩减法缩减波箔的自由度[20]，在保证计算精度的同时可大幅提高计算效率，且缩减后的波箔子结构可以视为超级单元。
组装所有的顶箔单元方程及缩减后的波箔子结构方程得到箔片整体的控制方程为
(10)
对于平滑近似代数摩擦模型，(10)式可简化为
(11)
(10)，(11)式的边界条件为轴承尾端的波箔和顶箔均焊接固定于轴承座。
1.4 接触摩擦的控制方程
根据波箔缩减策略，采用点对二维梁单元的接触模型来表征箔片结构的接触摩擦行为。接触力包括法向接触力fn和切向摩擦力ft，基于赫兹接触模型得到法向接触力fn为
(12)
式中：k为法向接触刚度系数；
δ为接触变形；
c为法向接触阻尼系数。
切向摩擦力ft分别采用LuGre动态摩擦模型和平滑近似代数摩擦模型表征。
1.4.1 LuGre动态摩擦模型
LuGre动态摩擦模型可以认为是Dahl模型的改进模型，可以表征斯特里贝克效应以及相应的黏着-滑移运动[21]。对于单个节点摩擦，LuGre动态摩擦模型的控制方程为
(13)
(14)
式中：σ0，σ1，σ2分别为摩擦刚度系数、微观阻尼系数和黏性阻尼系数；
zt为接触点的内部摩擦状态变量，描述了接触表面上的微观变形；
vt为两接触体的相对速度；
vc为转换速度，应尽可能小，以保证摩擦模型具有较好的微观阻尼[22]；
g(vt)为表征库仑摩擦和斯特里贝克效应的函数。
g(vt)一般可表示为
g(vt)=fsl+(fst-fsl)e-(vt/vs)2，
(15)
fsl=μffn，
式中：fsl，fst分别为滑动摩擦力和静摩擦力；
vs为斯特里贝克速度；
μf为摩擦因数。
定义黏着系数γ=fst/fsl表征斯特里贝克效应。当滑动速度为常数时，稳态切向摩擦力ft,s可以表示为
ft,s=g(vt)sgn(vt)+σ2vt，
(16)
这与大多数代数摩擦模型一致。
将所有摩擦节点的状态变量定义为矢量z，得到LuGre动态摩擦的控制方程为
(17)
1.4.2 平滑近似代数摩擦模型
采用最常用的双曲正切函数近似模型[17]，切向摩擦力为
ft=μffntanh(rvt)，
(18)
式中：r为平滑系数，取值越大越接近经典库仑摩擦模型，但也将导致严重的数值困难，平滑系数的取值直接决定了模型的计算精度和效率。
1.5 系统耦合方程的建立及求解
当采用LuGre动态摩擦模型时，联立转子、气膜、箔片和接触摩擦的控制方程，得到系统的控制方程为
(19)
定义系统状态变量为
(20)
将(19)式整理为一般形式为
(21)
当采用平滑近似代数摩擦模型时，系统的控制方程为
(22)
定义系统状态变量为
(23)
同样可将(22)式整理为(21)式的一般形式。
基于MATLAB的变阶变步长的隐式积分器ode15s[23]对(21)式的所有状态变量同时积分求解，隐式积分器需要获得系统方程的雅克比矩阵，仿真程序中采用了解析的雅克比矩阵∂F/∂y提高计算效率。此外，本模型还忽略了顶箔与波箔之间的分离行为，并引入Gümbel条件修正，即在(3)式中忽略p 在气体轴承几何发散区会存在低压区，该区域顶箔与波箔之间的法向接触力趋近于0。对于LuGre动态摩擦模型，法向接触力过小会使(14)式中的g(vt)趋近于0，进而导致数值困难。需定义一个极小的正实数作为接触容差，δ<的节点对应的摩擦状态自由度和方程需要被删除掉，并在δ≥时添加回系统方程。在仿真中还需要执行一个检测函数，当检测到接触状态变化时终止仿真程序，并在修改系统方程自由度后从断点继续仿真。此外，LuGre控制方程对系统雅克比矩阵的贡献项也需要进行较为繁琐的推导与编程。故LuGre动态摩擦模型会增加系统自由度，且程序实现复杂。 而平滑近似代数摩擦模型既不会增加系统自由度，也无需在仿真中调整系统方程，程序实现简单。 2种摩擦模型的切向摩擦力随两接触体的相对速度变化如图3所示，LuGre动态摩擦模型未考虑黏性阻尼项，且摩擦力ft基于滑动摩擦力μffn进行了归一化。由图3可知：LuGre动态摩擦模型能完全表征静摩擦(节点黏着状态)，平滑近似代数摩擦模型理论上静态时的摩擦力为0，仅能依靠较低的滑移速度来近似表征静摩擦。 图3 2种摩擦模型的切向摩擦力随两接触体的相对速度的变化趋势 以某波箔动压气体轴承为研究对象，其主要结构参数见表1。箔片材料参数为：密度7 800 kg/m3，弹性模量207 GPa，泊松比0.29。运行条件：pa=101 325 Pa，μ=1.95×10-5Pa· s。LuGre动态摩擦模型基本参数见表2。气膜周向和轴向的单元数分别为50和20，顶箔被划分为100个梁单元。此外，k=1×106N/mm，c=1×102N· s/mm，=1×10-12mm，β=1×10-4s。ode15s函数的相对容差设为1×10-6。为方便对比，黏性阻尼和黏着系数的取值使LuGre动态摩擦模型的稳态摩擦力和库仑摩擦完全一致。代数摩擦模型中的平滑系数r取10 s/mm[17]，这一取值已经导致了严重的数值困难，平均计算时间约为LuGre摩擦模型的5倍，且摩擦因数越大计算效率越低。系统建模及求解方法的准确性已在文献[20]中得到验证。 表1 波箔动压气体轴承主要结构参数 表2 LuGre动态摩擦模型基本参数 2.3.1 转子静态性能 本节所研究的静态性能为在不同静载荷下转子的静态偏心率，转子转速为30 000 r/min，摩擦因数为0.1，主要任务是求解(21)式的静平衡方程F(y0)=0得到其静态解y0。由于箔片的摩擦节点较多，导致系统在同一工况下存在多个静态平衡点，且无法直接将方程在其平衡点附近进行线化迭代求解； 2种摩擦模型转子的静态偏心率随静载荷的变化曲线如图4所示(为方便对比，图中还给出了无摩擦的计算结果， 即平滑近似代数摩擦模型的理论结果)：1)随静载荷增加，2种摩擦模型的转子静态偏心率差异逐渐增大； 图4 2种摩擦模型的转子静态偏心率随静载荷的变化曲线 2.3.2 转子稳定性 所研究的稳定性包括不同摩擦因数下系统的局部稳定性(系统在平衡位置处受小扰动后回到静态平衡位置稳定运行的能力)和全局稳定性(系统受任意非破坏性的冲击后回到静态平衡位置稳定运行的能力)[24]。转子质量为10.2 kg，静载荷为100 N，不平衡量为0。 局部稳定性和全局稳定性的关系取决于Hopf分岔行为(系统平衡点局部失稳时，根据系统参数的不同所形成的空间分岔形式包括超临界Hopf分岔和亚临界Hopf分岔)，如图5所示：1)对于超临界Hopf分岔，在稳定极限环存在的转速区域内，任意初始条件的转子轨迹都将收敛于该稳定极限环， 系统的全局稳定性和局部稳定性一致； (a)超临界Hopf分岔 局部稳定性的计算需要确定转子的平衡位置并施加小扰动。平衡位置通过时域仿真在固定转速下从坐标原点准静态逐步加载得到。随后，给予转子1‰的扰动并仿真足够长的时间来判断系统的局部稳定性。若转子轨迹收敛，系统在该转速下局部稳定； 全局稳定性的计算主要关心的是受各种非破坏性冲击后的响应，冲击效果可以等效为转子位置主导的初始条件。对于一般的非破坏性冲击，转子初始位置基本都在轴承间隙圆内。转子初始时在固定转速下被约束到间隙圆内的不同位置，并根据此时轴承力、静载荷和约束力达到静平衡的条件，计算得到初始的气膜压力、箔片变形以及摩擦状态。随后，移除约束力释放转子并仿真足够长的时间来判断系统的全局稳定性。若从所有初始位置释放的转子都回到了静态平衡位置，系统在该转速下全局稳定； 2种摩擦模型局部失稳转速随摩擦因数的变化曲线如图6所示：1)在低摩擦因数区，2种模型局部失稳转速一致，这是因为箔片摩擦节点在低摩擦因数下主要处于滑移状态，2种摩擦模型均能很好地表征动力学行为； 图6 2种摩擦模型局部失稳转速随摩擦因数的变化曲线 2种摩擦模型全局失稳转速随摩擦因数的变化曲线如图7所示：1)在低摩擦因数区，2种模型全局失稳转速一致，这是由于几乎所有箔片摩擦节点均处于滑移状态； 图7 2种摩擦模型的全局失稳转速随摩擦因数的变化曲线 对比分析了波箔动压气体轴承建模中LuGre动态摩擦模型和平滑近似代数摩擦模型在程序实现、稳态特性以及性能预测方面的差异，得到以下结论： 1)LuGre动态摩擦模型可以准确表征黏着状态，在任何仿真工况下预测结果均比较准确，且计算效率较高； 2)平滑近似代数摩擦模型无需增加系统自由度，程序实现简单，且在低速准静态或低摩擦因数、大扰动等以滑移为主的工况下预测结果比较准确；2.1 程序实现
2.2 稳态特性
2.3 性能预测
此外，平滑近似代数摩擦模型直接求静态解无意义，因为在理论上该模型的静态摩擦力为0：故计算的平衡点为通过时域仿真在固定转速下从坐标原点准静态逐步加载得到。
2)缓慢加载情况下的平滑近似代数摩擦模型处于LuGre动态摩擦模型和无摩擦模型之间，即在非零的低速滑移情况下平滑近似代数摩擦模型可以在一定程度上表现出接近黏着特性的行为。
2)对于亚临界Hopf分岔， 在不稳定极限环存在的转速区域内， 扰动较小的转子轨迹将收敛于稳定的平衡点，而扰动较大的转子轨迹将向外发散或收敛于外侧可能存在的稳定极限环，系统局部稳定，全局不稳定。
反之，则为局部不稳定。随转子转速升高，系统在某一转速下彻底失去局部稳定性，这一阈值转速即为局部失稳转速。
反之，则为全局不稳定。类似地，随转子转速升高，系统在某一转速下彻底失去全局稳定性，这一阈值转速即为全局失稳转速。
2)随摩擦因数增大，平滑近似代数摩擦模型预测的局部失稳转速先大于LuGre动态摩擦模型， 然后小于LuGre动态摩擦模型，这是由于随摩擦因数增大， 处于黏着状态的摩擦节点增多，平滑近似代数摩擦模型无法准确描述黏着状态。上述分析说明平滑近似代数摩擦模型在低摩擦因数区可准确预估局部失稳转速，在中摩擦因数区会高估局部失稳转速，在高摩擦因数区会低估局部失稳转速。
2)在中摩擦因数区，2种模型全局失稳转速出现差异，但没有局部失稳转速大，这是因为与局部的小扰动相比，全局稳定性计算中转子的扰动幅度较大，箔片摩擦节点大多处于滑移状态；
3)在高摩擦因数区，即使是大幅扰动下很多节点也主要处于黏着状态，平滑近似代数摩擦模型的全局失稳转速明显低于LuGre动态摩擦模型。
但该模型会增加系统自由度，且需要对零接触力进行特殊处理并执行检测函数，导致其程序实现较复杂。
但该模型无法准确表征黏着状态，在高摩擦因数、小扰动等以黏着为主的工况下预测结果偏差较大，且计算效率较低。