赵润超，焦映厚，曲秀全，张 赛，武祥林

(振动与噪声控制实验室(哈尔滨工业大学)，哈尔滨 150001)

燃气轮机周向拉杆转子是通过数对环状分布的拉杆，由预紧力将轮盘压紧而组合成的复杂结构. 轮盘间由预紧力引起的接触效应、机组结构的非连续性是燃机转子非线性行为的主要来源. 因此，为准确描述原模型特征，在建模过程中须考虑轮盘接触效应的影响.

Greenwood和Williamson提出了一种基于统计学分布的粗糙面弹性接触模型(GW模型)，给出了接触变形与表面形貌的关系，建立了弹性及塑性接触的判别准则，该模型可较为准确地反映粗糙表面的接触特性[1]. Pullen和Williamson提出了塑性接触处理方法，详细推导了高度分布的关系并与实验观测的进行对比，研究了粗糙表面塑性接触的特性[2]. 汪光明等提出了一种将拉杆转子部件相互连接处虚拟成铰接头和抗弯弹簧的力学模型，用动态子结构法对拉杆转子的横向振动进行理论计算[3]. 李辉光等基于弹塑性理论对具有粗糙表面的长方微元体有限元接触分析，给出了计算粗糙面接触刚度的方法[4]. 程礼等将接触面虚拟为具有非线性刚度项的抗弯弹簧，动力学结果表明转子结构的不连续和接触效应是造成转子振动双稳态特征的主要因素[5]. 石坤等针对机械结合部提出了一种观念虚拟材料的参数化模型，根据GW模型和Hertz接触理论，推导出了虚拟材料的弹性模量和泊松比[6].

目前的接触效应处理大多只考虑弹性或塑性过程，忽略了以弹塑性变形为主的接触过程. 且大多学者将轮盘接触条件处理为抗弯弹簧-铰链，在有限元模型编制时处理较为复杂. 对此，本文考虑弹性-弹塑性-塑性全部变形过程，发展了一种基于虚拟材料层的燃机转子轮盘接触模型. 通过控制虚拟材料层的材料参数(弹性模量、泊松比、密度、厚度)变化反映转子轮盘的力学接触特性.

1 轮盘表面单峰接触模型

在单微凸体的接触中，接触法向变形量ω的变化对接触面积、接触载荷、接触刚度影响很大. 随着变形量ω的增大，接触状态由弹性接触到弹塑性接触，再到塑性接触变化. 弹性接触和塑性接触可采用赫兹接触理论解释，而弹塑性接触理论涉及诸多不确定因素，很难通过现有理论解释. 然而，在接触状态发生变化时，其力学性质和变形量是连续且光滑过渡的，这为弹塑性的接触理论提供了一个解决思路，接触面积随变形量变化的示意图如图1所示.

图1 接触状态过渡示意图

1.1 单峰弹性接触阶段

当变形量ω小于弹性转为弹塑性变形的临界点ωe时，接触区域发生完全弹性变形，根据Hertz接触理论，接触面积、接触载荷、平均接触压力可表示为

ae=πR′ω,

(1)

(2)

(3)

式中:R′为模型简化后微凸球的半径，E′为模型简化后的弹性模量.

1.2 单峰塑性接触阶段

当变形量ω大于弹塑性转为塑性变形的临界点ωp时，接触区域发生完全塑性变形，材料力学性能将不能完全恢复到变形前状态[7]. 根据塑性变形接触理论，接触面积、接触载荷、平均接触压力可表示为

ap=2πR′ω,

(4)

fp=Hap,

(5)

(6)

式中，H为两接触体中表面硬度较低值.

1.3 单峰弹塑性接触理论

在弹塑性过渡临界点ωe和ωp处，接触面积和平均接触压力的过渡是连续且平滑的，不会产生突变，接触面积和平均接触压力的数值和其变化率分别等于弹性状态和塑性状态下的值. 构造关于接触面积aep的三次Hermite插值多项式aep(ω)：

aep(ω)=πRΓ(ω),

Γ(ω)=a1ω3+a2ω2+a3ω+a4.

(7)

代入边界条件，求解得Γ(ω)中各系数为

假设接触圆区域的半径rep与接触变形量ω的关系为

(8)

式中，δ是弹塑性接触半径的变系数.

根据接触力学理论[8]，弹塑性阶段平均接触压力pep和变形量ω的有如下关系：

(9)

式中，m1和m2是与接触半径rep相关的常量.

将式(8)代入式(9)整理得

pep=m3+m4lnω.

(10)

根据平均接触压力的过渡平顺且光滑的条件，可得平均接触压力的边界方程为

(11)

联立式(11)两方程，求解系数m3、m4并代入式(10)得

(12)

弹塑性阶段接触载荷的表达式为

fep=aep·pep=

(13)

2 轮盘表面实际接触分析

假设轮盘接触面的名义接触面积为An，其表面的单峰微凸体个数为N, 则在超过一定高度范围内，单峰微凸体的数量期望为nd.

N=ηAn,

(14)

(15)

式中:η为微凸体个数分布的面积密度，Φ(z)为概率密度函数，d为微凸球平均表面与另一表面的距离.

大量实验结果表明，接触面各微凸体的高度值服从高斯正态分布. 概率密度函数Φ(z)可用正态分布函数表示为

(16)

式中:μ为两轮盘接触面的平均高度，z为微凸体高度，σ为两接触面间高度分布标准差.

接触面的整体力学特性是全部单峰微凸体力学特性按照分布函数的积分，因此接触面整体力学特性可表示为弹性阶段、弹塑性阶段和塑性阶段的积分之和.

轮盘接触面的实际接触面积为

(17)

轮盘接触面的实际接触载荷为

(18)

根据法向刚度的定义，轮盘接触面的法向接触刚度为：

(19)

将式(13)代入上式，整理轮盘接触面的法向接触刚度为

(20)

(a)接触载荷不同模型对比 (b)不同塑性指数接触刚度对比 (c)接触刚度与载荷关系

图2(a)给出了本文模型、ZMC模型和赫兹弹性模型的接触刚度随无量纲变形量变化关系的对比，结果表明赫兹弹性模型未考虑弹塑性变形，结果具有较大误差. 本文模型相较于ZMC模型而言，对弹塑性及塑性变形的发生更为敏感，能够良好地反映实际接触状态.

单位面积接触刚度在不同塑性指数下变化规律如图2(b)所示，接触刚度反映了接触载荷对变形量变化的敏感程度. 在塑性指数较低时，轮盘接触面的法向接触刚度的变化率较小，随着塑性指数的不断增加，法向接触刚度的变化率逐渐变大.

在给定接触面的材料参数和表面形貌参数后，其接触的法向刚度与接触载荷为一一对应的关系，如图2(c)所示，在接触变形量较小时，即接触载荷较小时，法向接触刚度值迅速增加，随着接触载荷逐渐增大，法向接触刚度的变化趋势区域平稳. 对于高塑性的材料，法向接触刚度在接触载荷足够大时几乎相等，而对于低塑性材料，法向接触刚度在接触载荷足够大时值分离度较大.

3 虚拟材料层参数计算

目前对燃机转子接触效应的研究大多是将接触面的刚度转换为无厚度的铰链-弹簧模型，但此模型需要实验标定刚度参数且不足以全面反映轮盘接触面的力学特征. 另外，刚度不是材料的固有参数，弹性模量、剪切模量、泊松比、密度才是材料的固有弹性参数[9]. 本文引入虚拟材料层建立轮盘接触面模型，通过虚拟层的材料参数(弹性模量E、泊松比μ、密度ρ、厚度l)表征接触力学特性(图3).

图3 含虚拟材料的接触界面示意图

3.1 虚拟材料层弹性模量

图4 虚拟材料层等效变形量示意图

在实际轮盘接触层中，对模型单元施加法向载荷，总形变Δlfact由三部分组成：左轮盘材料的变形Δl1、接触层刚度单元变形Δl2和右轮盘材料变形Δl3，其中弹性层的厚度远远小于左右轮盘单元的厚度，将其近似为无厚度的弹簧. 在虚拟材料层中，由法向载荷引起的变形为Δlvirtual. 各变形量为

(21)

(22)

式中:A为两轮盘名义接触面积，K为轮盘实际接触面接触刚度(K=k·A，k为接触层单位面积接触刚度).

(23)

3.2 虚拟材料层泊松比

在实际轮盘接触中，法向应变为εfact，横向应变为εfact′，两者关系为

(24)

式中，Δl为实际轮盘接触单元的变形量.

在等效的虚拟材料层中，法向应变为εvirtual，横向应变为εvirtual′，两者关系为

(25)

(26)

3.3 虚拟材料层厚度

假定轮盘表面的微凸体高度值服从高斯正态分布，微凸体高度在[-3σ,3σ]的范围内能够涵盖99.7%的微凸体数量分布，当虚拟材料层厚度h>6σ时，可认为已包括所有的接触微凸体[10]. 而在实际的固体微观表层结构中，材料表面组织状态、形貌轮廓、油脂附着等均会产生一定厚度，因此虚拟材料层的厚度要结合轮盘微观表层结构进行选择[9].

虚拟材料层将覆盖从表面污染层到轻微变形层，总厚度表达式为

(27)

取左右轮盘虚拟材料层的厚度相等，单侧虚拟材料层厚度应大于lmax，故整体虚拟材料层厚度应大于2lmax.

3.4 虚拟材料层密度

根据接触面处两端的质量在等效为虚拟材料层前后保持不变的关系，推导得出虚拟材料层密度的表达式. 假设左右轮盘单元的长度分别为l1、l2，质量分别为m1、m2，虚拟材料层的长度为l，质量为m，各段质量为：

(28)

根据单元质量在等效为虚拟材料层前后保持不变的关系m=m1+m2，因假设l1=l2=l、ρ1=ρ2=ρ，整理可得虚拟材料层的密度ρvirtual为

(29)

4 算 例

本文建立了总长1 000 mm，具有5个压气机轮盘、1个涡轮盘的燃机拉杆转子模型，如图5所示. 表1给出了本文燃气轮机周向拉杆转子模型的各轮盘接触面积的参数.

图5 燃气轮机拉杆转子模型爆炸图

表1 接触面的尺寸参数

根据虚拟材料层参数的计算公式，结合所建立的燃机拉杆转子模型，可计算得到各虚拟材料层的参数. 其中，虚拟材料层的厚度和密度为与原转子尺寸材料相关的固定值，设虚拟材料层厚度l=0.005 m，密度ρvirtual=7 800 kg/m3.

根据图6结果可以发现，随着预紧力的增大，虚拟材料层弹性模量和泊松比均呈上升趋势. 虚拟层1、2、3、4由于轮盘接触面积较为接近，等效后参数值也较为接近，虚拟层6、7接触面积一致，等效后弹性模量和泊松比数值相等.

表2 预紧力F=50 kN时各轮盘虚拟材料层材料参数

(a)弹性模量变化关系

(b)泊松比变化关系

表3给出了在刚性支撑下不同预紧力下及实体转子的临界转速结果，可以发现，随着预紧力增大，转子一阶、三阶、四阶临界转速均出现上升趋势. 其中四阶临界转速数值变化波动最大，在一定范围内的预紧力变化导致的临界转速最大偏差达2.31%. 对于不考虑接触效应的实体转子，其整体刚度提升，各阶转速发生明显变化. 结果表明由拉杆预紧力引起接触效应对临界转速变化有明显影响.

表3 不同预紧力下临界转速

为进一步深入探究拉杆预紧力对轴系动态特性的影响，引入Capone轴承模型，考虑轴系在非线性油膜力作用下动态特性. 采用Newmark-β数值积分法求解动力学方程，并通过分岔图直观表现动力学行为. 表4给出了轴系的部分计算参数.

表4 转子轴系部分参数

选取三种预紧力大小(F=30 kN、F=40 kN、F=50 kN)与实体的轴段，探究一定转速范围内预紧力导致的接触效应对轴系动态特性的影响.

如图7所示给出了不同预紧力条件下转子的分岔图. 在考虑预紧力的情况下，随着转速的增加，转子从周期1运动进入周期2运动状态，之后出现了带有周期2特征的混动运动，转速达到10 000 r/min附近重新回到周期2运动，之后出现了混沌或概周期的运动状态. 分析不同预紧力大小可以发现，预紧力较小时，维持周期1运动的范围较大，且刚刚出现概周期运动时不宜转化为混沌状态. 因此，可根据转子实际工作转速调节拉杆预紧力，避开混沌或概周期运动状态. 在图7(d)实体接触状态的分岔图中，转子首先进入概周期分岔运动，之后便式中处于混沌状态，随着转速增加，混动运动的无量纲复杂呈现发散的状态，此时出现了油膜振荡失稳.

(a)F=30 kN分岔图

(c)F=50 kN分岔图

(b)F=40 kN分岔图

(d)不考虑接触实体转子分岔图

5 结 论

1)结合赫兹接触理论和变形光滑的边界条件，推导了弹塑性接触阶段接触面积、平均接触压力和接触载荷的表达式. 分析了不同塑性指数下法向接触刚度与变形量、接触载荷的关系，发现塑性指数较大的材料，法向接触刚度对接触间距的变化更为敏感.

2)引入虚拟材料层反映轮盘接触界面的接触效应，基于材料力学变形关系得到虚拟材料层的参数表征，分析结果表明虚拟材料层弹性模量和泊松比随预紧力的增大而增大，且对接触面积影响显著.

3)刚性支撑下预紧力变化导致的接触效应对第四阶临界转速的影响较为明显，转速最大偏差达2.31%，而不考虑接触效应的实体转子临界转速会出现较大偏差.

4)引入非线性油膜力，轴系动力学分岔图显示预紧力增大导致运动状态提前发生，实体转子与考虑预紧力接触的转子动态结果有较大分歧.