张占立 1，2，王 锋 1，王恒迪 1，熊明照 1
（1.河南科技大学 机电工程学院，河南 洛阳 471003；2.机械装备先进制造河南省协同创新中心，河南 洛阳 471003）
摘 要：为提高轴承试验机主轴部件临界转速设计计算的准确性和简便性，开发了基于传递矩阵理论和C#编程语言的临界转速求解模块。该模块可以将支撑轴承刚度、轴段剪切应力、回转效应和摆动惯性、轴承内圈及附加零件等因素加入到传递矩阵中，采用逐步扫描和“拟黄金分割法”自动切换的方法搜索临界转速，并且通过曲线动态显示搜索结果，搜索过程中求解速度可以控制。最后以轴承型号NU2209的试验机主轴为例进行求解，实现临界转速的快速精确求解，并将结果与ANSYS分析数据进行对比。实践表明：该模块计算准确、界面友好、操作简单、运行速度快、占用内存少，完全可以满足轴承试验机对临界转速的分析要求。
关键词：轴承试验机；C#；传递矩阵法；临界转速；转子动力学
轴承试验机是轴承性能试验的专用设备，其主轴部件的临界转速是轴承试验机的重要性能指标之一。传递矩阵法是计算主轴临界转速的重要方法，传递矩阵的阶数不会随着系统自由度的增加而增加，具有编程简单、运行速度快、占用资源少的特点，特别适合于链式系统[1-3]。如果将支撑轴承刚度、剪切效应、回转效应和摆动惯性等因素加入到传递矩阵中，可以更加接近真实工况。因此传递矩阵法在转子动力学的计算中占据重要地位[4-8]。C#是微软公司最新推出的一款面向对象的高级程序开发语言，吸收了VB简单可视化操作和C++运行高效的优点，以其强大的操作能力、创新的语言特性，成为.NET平台的旗舰语言。针对目前轴承试验机的设计需求，提出采用C#编制传递矩阵算法，基于传递矩阵理论，开发专用动力学分析软件，力求软件运行的稳定、快捷、操作方便，为试验机主轴设计提供技术支持。该设计方案对同类轴系部件的临界转速求解也具有借鉴意义。
文献[9]传递矩阵是一种广泛应用于计算转轴临界转速的数值计算法，经过多年的发展能将很多因素体现在计算程序中。其缺点是试算高阶模态时，计算精度降低，会出现数值不稳现象。该算法分析轴承试验机主轴的前三阶临界转速精确度较高，完全满足轴承试验机主轴设计的需要。该模块采用集总质量模型，将主轴转化为若干段具有无质量但有弹性的轴段和若干个集中质量组成的链状结构[10]，如图1所示。
图1 试验机主轴集总质量模型
Fig.1 Lumped Mass Model of Test Machine Spindle
通过对等截面轴的研究表明，若要求该简化模型带来的误差小于1%，传递矩阵分段数N应满足下式：
式中：R—需计算临界转速的最高阶数。
分段太多累计误差效应明显，同时降低了计算精度。其任一段的状态向量可以由截面的径向位移y、挠角θ、弯矩M、剪切力Q表示。记作D=（y θ M Q）T，根据力学平衡条件和材料力学理论，可以将轴段两侧状态向量的传递关系合写为矩阵（2），称为轴段传递矩阵。
式中
轴段截面积；Kt—截面系数（实心取0.886）；E—材料弹性模
量；I—轴的截面矩；L—轴段长度。
根据达朗伯原理将集中质量点两端状态向量传递关系写为矩阵（3），称为点传递矩阵。
式中：K—支撑刚度系数；ω—轴正进动速率；Id—直径转动惯量Id
则最右端和最左端传递关系可表示为式（4）：
根据式（2）可知，一个无长度的弹性轴段是一个单位矩阵，不影响计算结果，为了方便程序开发，在最左端添加一个长度为0的弹性轴段 L0，如图 1所示。将式（2）、式（3）代入式（4）得到添加该弹性轴段L0的总体传递关系，如式（5）所示。
添加L0后，两端都为自由端的边界条件均为（y θ 0 0）T。将边界条件代人式（5）后得式（6）。
将式（6）第三行、四行展开即可以得到转子同步正进动频率方程（7）。
图2 剩余量 △（ω2）曲线
Fig.2 Remaining Quantity Curve of △（ω2）
不同 ω代入式（7），求解函数 △（ω）的值，称为“剩余量”，使该剩余量为0的ω即为固有频率，剩余量△（ω2）值和ω函数曲线，如图2所示。
该模块选用微软公司的VS2010做为开发环境，设计模块的主界面，如图3所示。该模块将主轴各种尺寸、材料等参数汇集到数据库中，便于管理，界面简单友好、操作简单。用户只需要点击选择需要分析的主轴型号，系统自动联动搜索该主轴材料的相关参数；输入初始化步长；进行相关设置即可。
图3 模块主界面
Fig.3 The Main Interface Module
程序按照用户选择考虑的因素（包括刚度、剪切应力等因素）配置传递矩阵。按照边界条件建立频率方程，采用当前试算值加步长的方式进行搜索，当搜索到有根区间时，自动切换为“拟黄金分割法”进行缩小区间，直至搜索到满足精度的解。如果当前解不是所求临界转速的最高阶，继续按照原定步长进行下一个有根区间的搜索，直到完成所有解的搜索，且整个求解过程通过曲线动态显示。软件原理，如图4所示。
图4 流程及原理图
Fig.4 Process Flow Diagram and Schematic Diagram
为了使该模块操作简单、计算准确，需要具有独立运算传递矩阵、自动切换搜索方法、曲线显示搜索状态、准确判断解、求解速度可控等功能。下面分别进行阐述。
4.1 基于C#的传递矩阵运算技术
为了降低对第三方软件的依赖程度，保证模块的独立性和稳定性。矩阵算法采用C#来实现。采用长度为16（4行×4列）的三维数组记录每个轴段的参数，第一维为行数，第二维表为列数，用第三维表示轴段编号。实现从第N段向第1段的传递矩阵的相乘运算，从而避免调用matlab，提高了模块的独立性。代码如下所示：
代码中 b［i，j，n］为第 n 个轴段的传递矩阵，a［i，j，n+1］为第（n+1）到第N段总传递矩阵。将n从N循环到0，就实现了轴段振动状态从右端向左端的传递。
4.2 搜索算法切换技术
采用扫描法和“拟黄金分割法”相结合的搜索方式。从0开始，以指定初始化步长ω0开始试算。并且判断第n次剩余量的正负号，将该值赋值给数组C［n］（该数组用来记录每次搜索的状态），然后以 ωn+1=ω0+ωn（其中 ωn+1为第（n+1）次搜索值，ωn为第 n次搜索值）为下一次试算值进行第（n+1）次搜索，相邻两次剩余量异号则判断为有根区间，在曲线上表现为第一次穿越x轴。然后切换“拟黄金分割法”（“拟黄金分割法”指在已知区间（ωn，ωn+1）内用ωn+1-0.618ω0和频率方程该区间内理论解两个值做为插入点，通过判断两个插入值和两端的正负来消去一个区间的搜索方法），将 △ω0=-0.618ω0做为下一步长，以 ωn+1=ωn-0.618·ω0作为下一步试算值进行反向搜索，如果相邻两次剩余量仍同号继续按照△ω0=-0.618ω0反向搜索，如果出现异号，再次用“拟黄金分割法”以 △ω0=-0.618（-0.618ω0）作为步长正向搜索，以此类推不断反复缩小区间。在曲线上表现为第一次穿越x轴之后的无规律的反复振荡部分。若当前搜索步长小于指定精度，即 △ω0≤△，则认为搜索到第i阶临界转速D［i］=ωi。如果i小于指定搜索的阶数，系统扔采用ω=ω0+ωi进行下一有根区间搜索，直到将所有满足要求的临界转速搜索完毕。
4.3 剩余量曲线显示技术
Chart是用来显示图表和曲线的控件。每次搜索状态赋值给数组C［9999］，然后通过chart控件来显示。因此数组C［9999］的长度不能太长，否则影响曲线刷新效果，经多次调试将长度定为四位数，既满足搜索需要也节省计算机资源。剩余量是采用科学计数法进行计算的，变化幅度较大，控件需要不断调整y轴所需显示的最大值，有时甚至超过了控件的上下限以至于无法正常显示而报错。为了防止控件的数据溢出，对剩余量进行单位化处理。程序只关心剩余量的正负号并不关心具体大小，因此将所有大于0的剩余量用1表示，所有小于0的剩余量用-1表示，即可的用1、-1来定性的表示剩余量。显示效果，如图5所示。
图5 剩余量曲线
Fig.5 Remaining Quantity Curve
4.4 解的判断技术
该模块采用“拟黄金分割法”反复求解，第一次满足精度△ω0≤△的解有可能出现在x轴以上也可能在x轴以下，在曲线上表现为当前曲线转折点和第一次穿越x轴的转折点在x轴同侧或者两侧。第一种情况：如果在同侧说明理论解在当前ωn的前方区间（ωn，ωn+△ω0）内，如果此时判断为第 i阶解 ωi=ωn切换为逐步扫描法继续按照ωn+1=ω0+ωi开始搜索下一有根区间，因为该区间再次跨过第i阶理论解区间（ωn，ωn+△ω0），出现再次判断当前区间为有根区间而重复搜索的错误。第二种情况：如果在两侧说明理论解在当前 ωn的后方区间（ωn-△ω0，ωn）内，此时按照ωn+1=ω0+ωi可以顺利向下一个有根区间搜索。因此搜索步长第一次小于指定精度△时，如果出现第一种情况需要再次按照“拟黄金分割法”△ω0=-0.618（-0.618·ω0）反向缩小区间搜索，直到出现第二种情况可判断为搜索到解D［i］=ωi=ωn。该判断方式在曲线上表示为当曲线第一次穿越x轴是从上方穿越开始反复震荡，必须从下面结束震荡；反之，必须从上方结束。将曲线中间反复震荡部分看作只是一次穿越x轴，从整体上观察曲线，曲线必须是从左上方开始“穿针引线法”方式穿越x轴，才能保证搜索正常进行。即搜索曲线简化后必须，如图2所示。
4.5 循环速度的控制技术
求解过程占用内存虽然不是很多，但需要进行大量的数据计算。占用CPU比较严重，如果直接用for循环进行搜索，主界面会进入假死状态，求解过程中无法快速响应用户的其他操作。解决该问题常用操作是多线程编程，但多线程编程宏观上是多任务同时执行，微观上CPU还是一个一个执行任务，因此采用多线程编程后台进行大规模数据运算时，同样造成主界面运行迟钝。且多创建线程会带来争夺资源等不稳定因素。因此该模块采用另外一种编程思路，将求解主程序作为一个子函数，采用timer控件进行循环调用，用户通过TrackBar滑块控件来控制Timer的周期，间接控制CPU的占用率。每次搜索程序运行完成后触发一个全局事件，当Timer控件收到该全局事件且等到Timer下一个周期时，再次循环求解，以此来控制求解速度，该求解方案在配置一般的电脑上也可以很流畅进行求解。
剪切效应指计算挠度曲线变形时因剪切应力而发生剪切变形，其影响了主轴刚度，所以不能忽略该因素。另外轴通常情况下做正进动，即弓形回转方向与轴的转动方向相同时，回转效应和摆动惯性相当于增加了恢复力矩，影响了主轴刚度，因此也不能忽略。
以NU2209型号轴承试验机主轴为例，忽略倒角等因素，采用集总质量模型，将轴分为19段，支撑设置在第1段和第10段处。模型，如图1所示。分为四组求解，第一组考虑轴承刚度、剪切效应、回转效应和摆动惯性；第二组考虑轴承刚度、剪切效应；第三组考虑轴承刚度、回转效应和摆动惯性；第四组考虑轴承刚度、陀螺力矩。然后对前三组支撑刚度分别选取1E7、2E7……100E7N/m共100个刚度数据，对这三组共300组数据采用该模块对前三阶临界转速求解。对第四组支撑刚度分别取1E7、10E7、20E7……90E7、100E7N/m共11组数据采用该系统ANSYS二次开发模块对前三阶临界转速求解分析。将四组数据求解结果进行绘图，如图6所示。
图6 前三阶临界转速求解结果
Fig.6 The Result of Critical Speed
从图中曲线可以看出，该模块求解前两阶临界转速和ANSYS求解结果几乎重合，第三阶临界转速传递矩阵求解结果比ANSYS结果偏高，并随着刚度和转速升高误差增加。误差不大于15%。分别对比其他曲线可以得出：轴的的剪切应力降低了轴段的刚度，因此降低了主轴的临界转速，转速和阶数越高越明显。回转效应和摆动惯性提高主轴的刚度，从而提高了临界转速。
（1）采用C#根据传递矩阵理论开发临界转速求解模块，计算准确、界面友好、操作简单、运行速度快、占用内存少，可以快速求解轴承试验机主轴的临界值。
（2）该模块可以根据用户需要，将轴承支承刚度、剪切效应、回转效应和摆动惯性等因素中需要考虑的加入到传递矩阵中，方便工程对比分析。
（3）该临界转速求解模块开发方案对同类轴系的临界转速开发具有一定的借鉴意义。
参考文献
［1］李超，刘延峰，艾丽昆.基于混合模型的转子临界转速计算［J］.振动与冲击，2010，29（11）：245-248.（Li Chao，Liu Yan-feng，Ai Li-kun.Calculation of rotor critical speed based on hybrid model［J］.Vibration and Shock，2010，29（11）：245-248.）
［2］朱金虎，翁世修，蒋书运.高频电主轴临界转速计算及其影响参数分析［J］.机械设计与研究，2005，21（1）：28-30.（Zhu Jin-hu，Weng Shi-xiu，Jiang Shu-yun.Calculation of critical speed of high frequency electric spindle and its influence parameters of analysis［J］.Mechanical Design And Research，2005，21（1）：28-30.）
［3］张小龙，何洪庆.涡轮泵转子的临界转速研究（Ⅳ）分布质量轴的传递矩阵法［J］.推进技术，2000，21（2）：52-55.（Zhang Xiao-long，He Hong-qing.The study on critical speed of the rotor of a turbine pump（Ⅳ）transfer matrix method of distributed quality shaft［J］.Propulsion Technology，2000，21（2）：52-55.）
［4］柴山，刚宪约，姚福生.计算多转子系统临界转速的整体传递矩阵法［J］.上海理工大学学报，2002，24（1）：8-12.（Chai Shan，Gang Xian-yue，Yao Fu-sheng.Calculate rotor system more critical speed of the whole transfer matrix method［J］.Journal of University of Shanghai for Science and Technology，2002，24（1）：8-12.）
［5］冀成，杨兆建，宋高峰.多轮盘转子系统临界转速的计算方法分析［J］.机械设计与制造，2012（12）：28-30.（Yi Chen，Yang Zhao-jiang，Song Gao-feng.Calculation method of critical speed of multi disc rotor system［J］.Mechanical Design and Manufacturing，2012（12）：28-30.）
［6］陈建平，崔永鹏.基于传递矩阵法的冷却塔风机传动轴系模态分析［J］.机械设计与制造，2015（7）：219-222，226.（Chen Jiang-ping，Cui Yong-peng.Based on transfer matrix method of cooling tower fan shaft system modal analysis［J］.Mechanical Design and Manufacturing，2015（7）：219-222，226.）
［7］孟杰，陈小安.电主轴动力学分析的传递矩阵法［J］.机械设计，2008，25（7）：37-40.（Meng Jie，Chen Xiao-an.Dynamic analysis of motorized spindle of transfer matrix method［J］.Machine Design，2008，25（7）：37-40.）
［8］杨光，王玉丹.传递矩阵法进行电主轴的动力学特性分析［J］.机械设计与制造，2001（3）：50-51.（Yang Guang，Wang Yu-dan.The transfer matrix method TCo analyze dynamic characteristics of motorized spindle［J］.Mechanical Design and Manufacturing，2001（3）：50-51.）
［9］李日朝，郝东旭.燃气涡轮机高速转子临界转速仿真［J］.鱼雷技术，2014，22（6）：457-460.（Li Ri-chang，Hao Dong-xu.Simulation of critical speed of high speed rotor for gas turbine［J］.Torpedo Technology，2014，22（6）：457-460.）
［10］罗忠辉，薛晓宁，何真.高速转轴多阶临界转速的精确计算［J］.机床与液压，2004（12）：81-82，75.（Luo Zhong-hui，Xue Xiao-ning，He Zheng.High speed shaft order critical speed more precise calculation［J］.Machine Toolamp;Hydraulics，2004（12）：81-82，75.）
Development of the Rotor System Dynamics Module Based on Transfer Matrix
ZHANG Zhan-li1，2，WANG Feng1，WANG Heng-di1，XIONG Ming-zhao1
（1.School of Mechatronics Engineering，Henan University of Science and Technology，He’nan Luoyang 471003，China；2.The Collaborative Innovation Center of Mechanical Equipment and Advanced Manufacture in He’nan Province，He’nan Luoyang 471003，China）
Abstract：In order to improve the design quality of the spindle component of bearing machine and shorten the cycle of design，The design transfer matrix algorithm is put forward by Using C#softward，The critical module scheme of speed spindle bearing testing machine based on the transfer matrix theory is developed.The module can added support bearing stiffness，shear effect in shaft segment，inertia in the process of Swing，gyroscopic effect，bearing inner ring and additional parts to the transfer matrix，using stepwise scanning andquot;Golden section lawquot;automatic switching method to search the exact solutions.Through dynamic curve displays the results in search and the loop speed can be controlled in the process of solving.As an illustration，the machine spindle of bearing model NU2209 finally is solved to implement of the critical speed quickly and accurately，and the results were compared with ANSYS analysis and the test data.It shows that The module has the advantages of accurate calculation，friendly interface，simple operation，fast running speed and less memory occupation，which can completely meet the requirements of bearing tester to analysis of the critical speed.
Key Words：Bearing Testing Machine；C#；Transfer Matrix Method；Critical Speed；Rotor Dynamics
中图分类号：TH16；TP315；TH133
文献标识码：A
文章编号：1001-3997（2017）11-0122-04
来稿日期：2017-05-08
基金项目：国家科技部国防科工委基金项目（JPT-125-171）；河南省杰出人才创新基金项目（144200510020）；广东省省级科技计划项目（2013B011301020）
作者简介：张占立，（1965－），男，河南孟州人，博士研究生，硕士生导师，教授，主要研究方向：轴承设计、特种加工
联系客服
