﻿ 基于虚拟材料层燃机转子接触面建模方法
 哈尔滨工业大学学报  2021, Vol. 53 Issue (1): 117-123  DOI: 10.11918/202006087 0

### 引用本文

ZHAO Runchao, JIAO Yinghou, QU Xiuquan, ZHANG Sai, WU Xianglin. Modeling method of gas turbine rotor contact surface based on virtual material layer[J]. Journal of Harbin Institute of Technology, 2021, 53(1): 117-123. DOI: 10.11918/202006087.

### 文章历史

Modeling method of gas turbine rotor contact surface based on virtual material layer
ZHAO Runchao, JIAO Yinghou, QU Xiuquan, ZHANG Sai, WU Xianglin
Lab of Vibration and Noise Control Laboratory(Harbin Institute of Technology), Harbin 150001, China
Abstract: To improve the accuracy of gas turbine pull-rod rotor modeling, a virtual material layer was introduced at the contact surface of the wheel to characterize the contact effect. By the mathematical statistical model and gaussian distribution function, the mechanical model of the normal contact stiffness of contact surface was obtained, and the contact characteristics under different parameters were analyzed. Based on the constitutive relation of the deformation of material mechanics, the elastic modulus and poisson′s ratio of the virtual material layer were derived, and the contact layer was modeled. The results of dynamic calculation show that, under rigid support, the contact effect caused by the change of preload has a significant influence on the fourth-order critical speed, and the speed deviation can reach 2.31%. The results of introducing nonlinear oil film force show that the increase of preloading force leads to the advance of all motion states, and the dynamic results of solid rotor and rotor considering the contact of preloading force are quite different. The introduction of virtual material layer simplifies the modeling process of gas turbine rotor considering the contact effect and reflects the weakening effect of the contact effect on the transverse vibration stiffness of the rotor.
Keywords: virtual material layer    contact effect    gas turbine    rotor    modeling

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

1 轮盘表面单峰接触模型

 图 1 接触状态过渡示意图 Fig. 1 Transition diagram of contact state
1.1 单峰弹性接触阶段

 ${{a_e} = {\rm{ \mathsf{ π} }}{R^\prime }\omega ,}$ (1)
 ${{f_{\rm{e}}} = \frac{4}{3}{E^\prime }{R^{\prime \frac{1}{2}}}{\omega ^{\frac{3}{2}}},}$ (2)
 ${{p_{\rm{e}}} = \frac{{{f_e}}}{{{a_e}}} = \frac{4}{{3{\rm{ \mathsf{ π} }}}}{E^\prime }{R^{\prime - \frac{1}{2}}}{\omega ^{\frac{1}{2}}}.}$ (3)

1.2 单峰塑性接触阶段

 ${{a_{\rm{p}}} = 2{\rm{ \mathsf{ π} }}{R^\prime }\omega ,}$ (4)
 ${{f_{\rm{p}}} = H{a_{\rm{p}}},}$ (5)
 ${{p_{\rm{p}}} = \frac{{{f_{\rm{e}}}}}{{{a_{\rm{e}}}}} = H.}$ (6)

1.3 单峰弹塑性接触理论

 $\;\;\;\;\;\;\;{{a_{{\rm{ep}}}}(\omega ) = {\rm{ \mathsf{ π} }}R\varGamma (\omega ),}\\ {\varGamma (\omega ) = {a_1}{\omega ^3} + {a_2}{\omega ^2} + {a_3}\omega + {a_4}.}$ (7)

 ${{a_1} = \frac{{({\omega _{\rm{e}}} + {\omega _{\rm{p}}})}}{{{{({\omega _{\rm{e}}} - {\omega _{\rm{p}}})}^3}}},}$
 ${a_2} = - \frac{{2(\omega _{\rm{e}}^2 + {\omega _{\rm{e}}}{\omega _{\rm{p}}} + \omega _{\rm{p}}^2)}}{{{{({\omega _{\rm{e}}} - {\omega _{\rm{p}}})}^3}}},$
 ${{a_3} = \frac{{2\omega _{\rm{e}}^3 - 2\omega _{\rm{e}}^2{\omega _{\rm{p}}} + 7{\omega _{\rm{e}}}\omega _{\rm{p}}^2 - \omega _{\rm{p}}^3}}{{{{({\omega _{\rm{e}}} - {\omega _{\rm{p}}})}^3}}},}$
 ${{a_4} = - \frac{{2\omega _{\rm{e}}^2\omega _{\rm{p}}^2}}{{({\omega _{\rm{e}}} - {\omega _{\rm{p}}})(\omega _{\rm{e}}^2 - 2{\omega _{\rm{e}}}{\omega _{\rm{p}}} + \omega _{\rm{p}}^2)}}.}$

 ${r_{{\rm{ep}}}} = {(\delta {R^\prime }\omega )^{\frac{1}{2}}},({\omega _{\rm{e}}} \le \omega \le {\omega _{\rm{p}}})$ (8)

 ${p_{{\rm{ep}}}} = {m_1} + {m_2}\ln \frac{\omega }{{{r_{{\rm{ep}}}}}},$ (9)

 ${p_{{\rm{ep}}}} = {m_3} + {m_4}\ln \omega .$ (10)

 $\left\{ {\begin{array}{*{20}{l}} {{m_3} + {m_4}\ln {\omega _{\rm{e}}} = kH,}&{\omega = {\omega _{\rm{e}}};}\\ {{m_3} + {m_4}\ln {\omega _{\rm{p}}} = H;}&{\omega = {\omega _{\rm{p}}}.} \end{array}} \right.$ (11)

 ${p_{{\rm{ep}}}} = H\left[ {1 - (1 - k)\frac{{\ln {\omega _{\rm{p}}} - \ln \omega }}{{\ln {\omega _{\rm{p}}} - \ln {\omega _{\rm{e}}}}}} \right].$ (12)

 $\begin{array}{l} {f_{{\rm{ep}}}} = {a_{{\rm{ep}}}} \cdot {p_{{\rm{ep}}}} = \\ {\rm{ \mathsf{ π} }}R\varGamma (\omega ) \cdot H\left[ {1 - (1 - k)\frac{{\ln {\omega _{\rm{p}}} - \ln \omega }}{{\ln {\omega _{\rm{p}}} - \ln {\omega _{\rm{e}}}}}} \right]. \end{array}$ (13)
2 轮盘表面实际接触分析

 ${N = \eta {A_n},}$ (14)
 ${{n_d} = N\int_d^\infty \varPhi (z){\rm{d}}z = \eta {A_n}\int_d^\infty \varPhi (z){\rm{d}}z.}$ (15)

 $\varPhi (z) = \frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}} \sigma }}\exp \left( { - \frac{{{{(z - \mu )}^2}}}{{2{\sigma ^2}}}} \right).$ (16)

 $\begin{array}{*{20}{l}} {{A_{{\rm{ real }}}} = {A_{\rm{e}}} + {A_{{\rm{ep}}}} + {A_{\rm{p}}} = \eta {A_n}[\int_d^{d + {\omega _{\rm{e}}}} {{a_{\rm{e}}}} \varPhi (z){\rm{d}}z + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \int_{d + {\omega _{\rm{e}}}}^{d + {\omega _{\rm{p}}}} {{a_{{\rm{ep}}}}} \varPhi (z){\rm{d}}z + \int_{d + {\omega _{\rm{p}}}}^{ + \infty } {{a_{\rm{p}}}} \varPhi (z){\rm{d}}z].} \end{array}$ (17)

 $\begin{array}{*{20}{l}} {{F_{{\rm{ real }}}} = {F_{\rm{e}}} + {F_{{\rm{ep}}}} + {F_{\rm{p}}} = \eta {A_n}[\int_d^{d + {\omega _{\rm{e}}}} {{f_{\rm{e}}}} \varPhi (z){\rm{d}}z + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \int_{d + {\omega _{\rm{e}}}}^{d + {\omega _{\rm{p}}}} {{f_{{\rm{ep}}}}} \varPhi (z){\rm{d}}z + \int_{d + {\omega _{\rm{p}}}}^{ + \infty } {{f_{\rm{p}}}} \varPhi (z){\rm{d}}z].} \end{array}$ (18)

 ${k_c} = \frac{{\partial f}}{{\partial \omega }}.$ (19)

 $\begin{array}{*{20}{l}} {{K_{{\rm{ real }}}} = {K_{\rm{e}}} + {K_{{\rm{ep}}}} + {K_{\rm{p}}} = \eta {A_n}[\int_d^{d + {\omega _{\rm{e}}}} {{k_{\rm{e}}}} \varPhi (z){\rm{d}}z + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \int_{d + {\omega _{\rm{e}}}}^{d + {\omega _{\rm{p}}}} {{k_{{\rm{ep}}}}} \varPhi (z){\rm{d}}z + \int_{d + {\omega _{\rm{p}}}}^{ + \infty } {{k_{\rm{p}}}} \varPhi (z){\rm{d}}z].} \end{array}$ (20)

 图 2 多峰接触分析 Fig. 2 Multipeak contact analysis

3 虚拟材料层参数计算

 图 3 含虚拟材料的接触界面示意图 Fig. 3 Schematic diagram of contact interface with virtual materials
3.1 虚拟材料层弹性模量

 图 4 虚拟材料层等效变形量示意图 Fig. 4 Schematic diagram of equivalent deformation of virtual material layer

 $\varDelta {l_{{\rm{fact}}}} = \varDelta {l_1} + \varDelta {l_2} + \varDelta {l_3} = \frac{{F{l_1}}}{{EA}} + \frac{F}{K} + \frac{{F{l_3}}}{{EA}},$ (21)
 $\varDelta {l_{{\rm{ virtual }}}} = \frac{{F{l_{{\rm{ virtual }}}}}}{{\bar EA}}.$ (22)

 $\bar E = \left( {\frac{{2lK}}{{2lK + EA}}} \right)E = \left( {\frac{{2lk}}{{2lk + E}}} \right)E.$ (23)
3.2 虚拟材料层泊松比

 $\begin{array}{*{20}{l}} {{\varepsilon _{{\rm{ fact }}}}^\prime = \mu \cdot {\varepsilon _{{\rm{ fact }}}} = \mu \frac{{\varDelta l}}{{({l_1} + {l_3})}} = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mu \left( {\frac{{F{l_1}}}{{EA}} + \frac{F}{K} + \frac{{F{l_3}}}{{EA}}} \right)/({l_1} + {l_3}),} \end{array}$ (24)

 ${\varepsilon _{{\rm{virtual }}}}^\prime = \bar \mu \cdot {\varepsilon _{{\rm{virtual }}}} = \frac{{F\bar \mu }}{{\bar EA}}.$ (25)

 $\bar \mu = \left( {\frac{{\bar E}}{E} + \frac{{\bar EA}}{{2lk}}} \right)\mu .$ (26)
3.3 虚拟材料层厚度

 ${l_{\max }} = \sum\limits_{i = 1}^n {\max } ({l_i}).$ (27)

3.4 虚拟材料层密度

 $\left\{ {\begin{array}{*{20}{l}} {{m_1} = {\rho _1} \cdot V({l_1}),}\\ {{m_2} = {\rho _2} \cdot V({l_2}),}\\ {m = {\rho _{{\rm{ virtual }}}} \cdot V(l).} \end{array}} \right.$ (28)

 ${\rho _{{\rm{ virtual }}}} = \frac{{\rho \cdot V(l) + \rho \cdot V(l)}}{{V(2l)}} = \rho .$ (29)
4 算例

 图 5 燃气轮机拉杆转子模型爆炸图 Fig. 5 Explosion diagram of gas turbine pull-rod rotor model

 图 6 虚拟材料层参数与拉杆预紧力关系 Fig. 6 Relationship between parameters of virtual material layer and preload force

 图 7 不同预紧力转子分岔图 Fig. 7 Bifurcation diagram of rotors under different preload forces
5 结论

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

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

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

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

 [1] GREENWOOD J A, WILLIAMSON J B P. Contact of nominally flat surfaces[J]. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences (1934-1990), 1966, 295(1442). [2] PULLEN J, WILLIAMSON J B P. On the plastic contact of rough surfaces[J]. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences (1934-1990), 1972, 327(1569). [3] 汪光明, 饶柱石, 夏松波, 等. 拉杆转子力学模型的研究[J]. 航空学报, 1993, 14(8): 419. WANG Guangming, RAO Zhushi, XIA Songbo, et al. The analysis of mechanical model of rod fastening rotor[J]. Acta Aeronautica ET Astronautica Sinica, 1993, 14(8): 419. DOI:10.3321/j.issn:1000-6893.1993.08.017 [4] 李辉光, 刘恒, 虞烈. 粗糙机械结合面的接触刚度研究[J]. 西安交通大学学报, 2011, 45(6): 69. LI HuiGuang, LIU Heng, YU Lie. Contact stiffness of rough mechanical joint surface[J]. Journal of Xi'an Jiaotong University, 2011, 45(6): 69. [5] 程礼, 钱征文, 陈卫, 等. 结构参数对拉杆转子双稳态振动特性的影响[J]. 振动测试与诊断, 2012, 32(5): 767. CHENG Li, QIAN Zhengwen, CHEN Wei, et al. The influence of structural parameters on the bistable vibration characteristics of the tie rod rotor[J]. Journal of Vibration, Measurement & Diagnosis, 2012, 32(5): 767. [6] 石坤, 宋俐, 师俊平. 机械结合部等效材料参数建立与试验[J]. 农业机械学报, 2014, 45(2): 297. SHI Kun, SONG Li, SHI Junping. Establishment and test of equivalent material parameters of mechanical joint[J]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(2): 297. [7] 赵永武, 吕彦明, 蒋建忠. 新的粗糙表面弹塑性接触模型[J]. 机械工程学报, 2007, 43(3): 95. ZHAO Rongwu, LÜ Yanming, JIANG Jianzhong. A new elastic-plastic contact model for rough surfaces[J]. Journal of Mechanical Engineering, 2007, 43(3): 95. DOI:10.3321/j.issn:0577-6686.2007.03.016 [8] 李玲, 蔡安江, 蔡力钢, 等. 螺栓结合面微观接触模型[J]. 机械工程学报, 2016, 52(7): 205. LI Ling, CAI Anjiang, CAI Ligang, et al. Micro-contact model of bolted-joint interface[J]. Journal of Mechanical Engineering, 2016, 52(7): 205. [9] 田红亮.机械结构固定结合部虚拟材料的动力学建模[D].武汉: 华中科技大学, 2011 TIAN Hongliang. Dynamic modeling on fixed joint interface virtual material in mechanical structure[D]. Wuhan: Huazhong University of Science and Technology, 2011 [10] 肖会芳, 孙韵韵, 徐金梧, 等. 螺栓固结界面接触刚度的虚拟材料建模[J]. 中南大学学报(自然科学版), 2019, 50(9): 2100. XIAO Huifang, SUN Yunyun, XU Jinwu, et al. Virtual material based modeling for contact stiffness characteristics of bolted joint interface[J]. Journal of Central South University (Science and Technology), 2019, 50(9): 2100.