接触应力对飞轮转子动力特性的影响

任正义，杨立平，曹志好，曹 飞

(1.哈尔滨工程大学工程训练国家级实验教学示范中心，黑龙江 哈尔滨150001；2.哈尔滨工程大学机电工程学院，黑龙江 哈尔滨150001)

1 引言

环境恶化和资源枯竭是现如今人类面临的巨大挑战，人类在发展的道路上，正在寻找节能环保高效的能源来代替传统能源。飞轮储能系统作为新能源领域的代表，越来越受到人们的重视［1］。飞轮储能系统是将电能转化为机械能再到电能的一种能量转化装置，主要由转子系统、高速永磁无刷电机、控制箱、支架和冷却系统等组成［2］。如今，人们对飞轮储能系统的研究日渐成熟。韩邦成等对工作转速为5000rpm的飞轮转子进行了模态分析，并通过试验验证有限元计算的结果［3］。陈强等对飞轮储能系统有限元分析，求解其临界转速和模态振型［4］。刘佩等对1MW/60MJ轴系转轴进行优化分析和模态分析，并未考虑过盈配合处应力的影响［5］。付雄新等提出了将安装基础运动转化为对转子激励的动力学求解方法来分析飞轮电池转子在路面输入引起安装基础振动变化情况下的振幅变化［6］。以上论文都是在无接触应力的基础上分析的，实际情况是轴系结构主要部件之间都是过盈装配，存在接触应力。因此建立有接触应力的高精度飞轮储能系统的动力学模型显得尤为重要。

2 飞轮储能系统

2.1 飞轮储能系统结构组成

该飞轮储能系统采用立式左右对称结构布局，主要包括机械部分和控制部分两大模块。机械部分的结构示意图，如图1所示。其中金属飞轮转动惯量0.59kg.m^2，最大额定转速为16000rpm，储能量为600Wh。电机采用电动/发电一体的永磁无刷高速直流电机，额定功率为7.4Kw，额定转速为16000rpm。该系统主要由循环水冷系统进行冷却。另外，电机的转子安装在转轴上，使结构更加紧凑，但也会使转子轴系结构过长。推力盘和轴向电磁轴承安装在转轴上部分，主要作用是使转子系统悬浮在空中。上下两端各有电磁轴承和保护轴承作为支撑，对转子系统的动力学特性至关重要［7］。

图1 600Wh飞轮储能系统结构布局Fig.1 Structure Layout of 600Wh Flywheel Energy Storage System

2.2 飞轮转自三维模型

飞轮转子系统是整个储能系统的关键部分，大约95%的储能量储存在这里。利用Soildworks将飞轮转子零件装配，得到的三维图，如图2所示。

图2 飞轮转子三维图Fig.2 Graphic Model of Flywheel Rotor

转子系统主要部分的材料及属性，如表1所示。

表1 轴系部件材料属性Tab.1 Material Properties of Shafting Components

3 接触应力分析

3.1 静态过盈量计算

该600Wh飞轮储能系统设计的额定转速为16000rpm，不平衡质量产生的离心力也会很大。为了保证轴上过盈配合处零件在高速旋转的时候与转轴之间不发生松脱，合适的过盈装配量显得尤为重要。轴上过盈配合处的结构简图，如图3所示。a表示轴径，b表示配合面的半径，c表示与轴过盈配合处部件外径，p表示过盈配合面之间的静压力，E1表示轴的弹性模量，E2表示配合件的弹性模量，x、y分别为转轴和配合件的径向位移，μ1表示轴的泊松比，μ2表示配合件的泊松比，r1表示轴的半径，r2表示配合件半径，Δ表示配合面的过盈量。

图3 过盈配合处的结构简图Fig.3 Structural Sketch of Interference Fit

根据弹性力学理论［8］得：

根据上式三个方程可以求解出r1=r2=b处配合面的接触压力为：

由于该飞轮转轴为实心，即a=0时，则配合处的接触压力为：

要满足使用要求，需满足接触压力值小于材料的屈服强度即：

根据关系式（4），以配合处最低屈服强度的材料属性参数代入，求得各过盈配合处最大过盈量，如表2所示。

表2 各配合处最大静态过盈量Tab.2 Maximum Static Interference at Each Mating Point

3.2 动态过盈量计算

当转子以角速度ω旋转时，由F=mrω2得到，飞轮毂、轴承转子电机转子和推力盘的离心力大于转轴的离心力，会使上述各部分的内径变形比转轴的外径变形大，两者之间存在一个变形差值。随着转速的增大，变形差越大，当差值大于两者之间的过盈量时，转轴和各部分之间会发生脱动。转动时与轴过盈配合处各部件的径向位移为y1，轴的径向位移为x1，配合面的过盈量减少值为Δ1，则［8］

动态过盈量为：

式中：ρ—材料密度。

以额定转速16000rpm为例，由表1查得各个值代入式（5）和式（6），计算出各过盈配合处不发生脱动的最小动态过盈量，如表3所示。

表3 各配合处最小动态过盈量Tab.3 Minimum Static Interference at Each Mating Points

4 转子系统有限元仿真与分析

4.1 转子系统有限元模型的建立

过盈配合是一种高度的非线性接触，利用Ansys Workbench进行模态分析时，Modal模块会自动的将非线性接触默认为线性接触，导致计算结果不准确。考虑轴系装配过盈量，需要先在Static Structural模块进行静应力分析，再将静应力分析结果导入到Modal模块进行模态求解。在进行有接触应力的仿真过程中将默认的线性接触改为非线性接触的方法是在过盈配合处创建接触对，在选择目标面时需要遵循的原则是：当凸面与平面或凹面接触时，应选择平面或凹面为目标面［9］。因此在静力学分析设置时，设置转轴的表面为接触面，设置飞轮毂、推力盘、电机转子和磁轴承转子的内表面为目标面。

为了精确计算模型，必须在设置接触对的过程中设置合理的FKN（接触刚度因子）值。对于绑定和不分离的线性接触类型的问题，默认FKN=10；对于以体积为主的非线性接触问题，默认FKN=1或选择“Program Controlled”；对于以弯曲问题为主的非线性接触，可以设置为FKN=0.01～0.1［9-10］。飞轮储能系统转子主要以弯曲为主，因此在设置FKN时大小设置为0.07。得到的600Wh飞轮转子网格划分结果，如图4所示。

图4 转子有限元模型Fig.4 Finite Element Model of Rotor

4.2 转子系统的有限元仿真结果分析

由于飞轮毂承载主要的能量，产生的力学变形主要发生在飞轮毂处，因此仿真时保持其他配合处过盈量不变，只改变飞轮毂和轴配合处的过盈量，得到的最大等效应力和前两阶临界转速结果，如表4所示。

表4 有限元仿真结果Tab.4 Finite Element Simulation Results

从表4可以看出，过盈量每增大0.01mm，等效应力大约增大16%左右，一阶和二阶临界转速也随之增大。其中对第一阶临界转速影响程度大于对第二阶临界转速的影响程度，过盈量每增加0.1mm，一阶临界转速增大12%左右。整体趋势是，过盈量增大，等效应力和前两阶临界转速呈增大趋势。当过盈量超过0.2mm时，最大等效应力超出了飞轮毂材料的屈服强度，飞轮毂开始变形，因此飞轮毂和转轴过盈配合处过盈量最佳取值范围是0.07mm～0.19mm。

4.3 应力理论值与仿真值对比分析

根据应力计算公式求得轴毂配合处各过盈量下应力的理论值，得到轴毂处应力的理论值和仿真值对比表，如表5所示。从表5中可以看出理论值与仿真值基本一致。

表5 应力理论值与仿真值对比分析Tab.5 Contrastive Analysis of Stress Theory Value and Simulated Value

4.4 转子系统的模态分析

在转轴与飞轮毂配合处设置过盈量为0.09mm，求解得到的等效应力图和模态振型图，如图5、6、7所示。

由等效应力图5可知，轴系的最大等效应力出现在飞轮毂与轴配合出最大值为187.41MPa，小于材料的屈服强度。由振型云图6、7可知，转子系统的一阶固有频率为64.62Hz，其振幅表现为从中心飞轮毂处向两端减小，最大相对位移为6.2436mm。这是一阶平动状态。转子系统的二阶固有频率为189.31Hz，其振幅表现为两端开始变大，中间飞轮毂处出现分层现象，最大相对位移为4.0498mm。这是二阶锥动状态。

图5 等效应力图Fig.5 Equivalent Stress Diagram

图6 一阶模态振型Fig.6 First Mode Shapes

图7 二阶模态振型Fig.7 Second Mode Shapes

5 结论

通过理论计算公式计算出轴上各主要部件与转轴过盈配合处过盈量的范围，并且建立了考虑接触应力的飞轮储能系统有限元模型，然后进行有限元分析。分析发现，接触应力对转子系统的临界转速有很大的影响。通过改变飞轮毂和轴的过盈配合量，得出该配合处过盈量的大致范围为0.07mm～0.19mm，为相同类型的飞轮储能系统设计提供了参考依据。