刚性接触网-受电弓系统快速仿真方法

陈龙; 刘志刚; 段甫川; 胡泽尧; 徐钊; 陈可

doi:10.3969/j.issn.0258-2724.20230206

刚性接触网-受电弓系统快速仿真方法

doi: 10.3969/j.issn.0258-2724.20230206

陈龙^1,,
刘志刚^1, ,,
段甫川¹,
胡泽尧¹,
徐钊²,
陈可¹

1.
西南交通大学电气工程学院，四川成都 611756
2.
同济大学交通运输工程学院，上海 201804

基金项目: 国家自然科学基金项目（52172408，52102478，U2468229）

详细信息

作者简介:
陈龙（1995—），男，博士研究生，研究方向为刚性接触网-受电弓系统动力学，E-mail：longchen@my.swjtu.edu.cn

通讯作者:
刘志刚（1975—），男，教授，博士，研究方向为高速铁路弓网动力学、状态检测与评估，E-mail：liuzg_cd@126.com

中图分类号: U225.1
计量
- 文章访问数: 102
- HTML全文浏览量: 81
- PDF下载量: 21
- 被引次数: 0
出版历程
- 收稿日期: 2023-05-10
- 修回日期: 2023-09-11
- 网络出版日期: 2024-11-12

Fast Simulation Method for Pantograph and Overhead Conductor Rail System

1.
School of Electrical Engineering, Southwest Jiaotong University, Chengdu 611756, China
2.
School of Transportation Engineering, Tongji University, Shanghai 201804, China

摘要

摘要:
针对当前刚性接触网-受电弓系统有限元模型仿真速度慢、计算时间成本高的问题，对采用三维接触算法的弓网仿真方法与流程进行改进. 首先，采用中心差分思想，将求解弓网接触副相对运动速度时需迭代计算的方程转换为可直接计算的显式方程；然后，将刚性接触网在静平衡处进行线性化处理，以避免刚度矩阵组装耗时，并加快刚性网内力计算；其次，对弓网接触状态进行惰性判断以减少计算量；最后，对本文所提快速仿真方法在不同情况下的计算效率与精度进行分析. 研究结果表明：在30跨8 m跨距的刚性接触网-受电弓仿真算例中，快速仿真方法相比标准仿真方法节省97.67%的仿真时间，且接触力结果最大偏差仅为0.48%；随着模型规模的增大，其节省的时间迅速增加，计算效率优势愈发显著，同时接触力结果偏差均小于1.0%；随着运行速度的提高，所节省的时间占比基本不变，接触力结果偏差略有增大趋势，在230 km/h以下的速度工况中，接触力标准差偏差均小于1.0%.
- 刚性接触网 /
- 有限元方法 /
- 计算效率 /
- 绝对节点坐标法 /
- 三维接触
Abstract:
Finite element model simulation of the pantograph and overhead conductor rail (OCR) system is slow, and the time cost of calculation is high. Therefore, the simulation method and process of the pantograph and OCR system using three-dimensional contact formulation were improved. Firstly, the equation that needs to be calculated iteratively when solving the relative velocity of the contact pair between the pantograph and the OCR was replaced with an explicit equation that can be calculated directly based on the central difference method. Then, the OCR model was linearized at the static equilibrium state to avoid the time-consuming stiffness matrix assembly procedure and increase the efficiency when calculating the internal force of the OCR. Next, a lazy judgment strategy was used to estimate the contact state of the pantograph and the OCR to reduce the computational load. Finally, the computational efficiency and accuracy of the fast simulation method in different cases were analyzed. The results show that compared with the standard simulation method, the proposed fast simulation method can save 97.67% of computational time in the example of pantograph and OCR with 30 spans of 8 m, and the maximum error of contact force results is only 0.48%. With the increase in the model scale, the time saved by the fast simulation method increases sharply, and its computational efficiency advantage becomes more and more significant. Meanwhile, the errors of contact force results are all less than 1.0%. With the increase in the operation speed, the proportion of time saved by the fast simulation method remains stable, and the errors of contact force results increase slightly. At speeds under 230 km/h, the standard deviation errors of contact force are all less than 1.0%.
- overhead conductor rail /
- finite element method /
- computational efficiency /
- absolute nodal coordinate formulation /
- three-dimensional contact

HTML全文

磁悬浮铣削电主轴是以磁悬浮轴承为支承的高速机床主轴部件，其将磁悬浮轴承的主动可控性、无摩擦、免润滑、低噪声、长寿命等优点与电主轴零传动、超高速的特点结合在一起，满足了高速切削加工的要求^[1-2]. 但由于磁悬浮铣削电主轴在高速切削时，切屑会沿着容屑槽进入刀具一起旋转，进而改变“主轴-刀具”系统的质量大小及分布，并由此引发“主轴-刀具”系统的惯性特性、模态参数以及系统动态响应随时间的非线性变化, 进而影响铣削加工的质量. 因此，进行磁悬浮铣削电主轴变质量系统的旋转特性与动态响应研究，以提高铣削加工的质量，具有重要的理论和应用价值.

现有国内外对磁悬浮转子系统动力学建模方法的研究，主要将磁悬浮轴承视为约束，将转子等效为盘-轴模型的Jeffcott转子，或是将其简化为质量不变且连续分布的梁-轴模型，如欧拉-伯努利梁、瑞利梁、铁木辛科梁等. 在具体研究时，针对细长轴，通常基于欧拉-伯努利梁理论进行建模；对于非细长轴，若同时考虑梁的剪切变形和梁弯曲变形所引起的转动惯量，则采用铁木辛科梁理论进行建模；若只考虑梁弯曲变形所引起的转动惯量，则基于瑞利梁理论进行建模. 主要研究成果有：钟志贤等^[3]为研究以主动磁悬浮轴承支承的高速转子的裂纹故障特征，依据Jeffcott转子理论建立系统的运动微分方程；吴超等^[4]以磁悬浮汽轮机多盘转子系统为研究对象，利用基于欧拉-伯努利梁理论的有限单元法，建立考虑盘轴耦合特性的磁悬浮转子系统的运动微分方程；Lei等^[5]利用基于瑞利梁振动理论的有限单元法，并考虑陀螺力矩的激励效应，建立磁悬浮柔性转子的运动微分方程；KUNG等^[6]基于欧拉-伯努利梁理论，利用Hamilton原理和Galerkin理论推导出由磁悬浮轴承支承的裂纹主轴系统的运动方程；Bouaziz等^[7]采用基于铁木辛科梁理论的有限单元法，建立磁悬浮铣削电主轴系统的运动微分方程等.

另外，有部分学者在磁悬浮转子动力学模型的基础上，围绕外部激励对系统动态响应的影响进行了深入研究. 如：Zhang等^[8]通过Newmark积分法求解磁悬浮柔性转子的运动微分方程，研究磁悬浮转子在移动载体上工作时，其基底运动激励引起转子振动的动态变化规律；陈小安等^[9-10]基于电磁学和机械系统动力学基本理论，建立各种偏心状态下高速电主轴的广义不平衡力表达式，分析了高速电主轴转子在电磁不平衡拉力与离心力作用下，系统在偏心状态下振动的幅值及频率特性；宋春生等^[11]针对磁悬浮柔性转子的解耦控制问题，利用有限元法建立磁悬浮柔性转子的动力学模型; 王艺宇等^[12]针对磁悬浮叶轮机械，采用有限元方法，基于铁木辛科梁单元模型建立考虑界面接触影响的磁悬浮轴承转子动力学模型.

上述研究成果都是以质量不变的转子系统为研究对象，而对质量时变转子系统的动力学建模与动态特性研究，目前国内外的成果主要集中在纺织机械和机床加工主轴领域. 如：彭超英等^[13-14]建立纺织机械变质量转子的动力学模型，并采用龙格库塔法进行求解，以分析变质量效应对转子振动响应的影响；袁龙翔等^[15]针对车削过程中工件质量时变导致工件的质量矩阵、陀螺矩阵、刚度矩阵同步时变的现象，采用转子动力学有限单元法，建立工件的时变动力学模型，并引入动态切削力，分析工件质量时变对工件系统的临界转速、车削加工稳定性区域的影响；Cao等^[16]针对磁悬浮车削电主轴在加工过程中由于工件材料的去除导致系统质量时变的现象，以车削加工三维移动切削力、切屑流失造成的反推力以及振动量变化引起的时变电磁力为激励，基于瑞利连续梁振动理论，利用假设模态法和拉格朗日方程建立磁悬浮车削电主轴系统动力学方程，分析质量时变对系统动态特性的影响；欧阳智海等^[17]针对磁悬浮铣削电主轴在切削过程中因切屑不断进入和离开刀具容屑槽中而导致“刀具-主轴”系统质量变化的现象，通过建立“切屑-刀具-主轴”系统的变质量动力学方程，分析不平衡磁拉力、惯性力等对系统动态特性的影响. 但这些研究成果都没有考虑当系统质量改变时，主轴系统在高速旋转下惯性载荷的非线性变化规律，以及由此所引起的系统新的动态特性.

基于此，本文首先利用连续梁振动理论的有限单元法建立磁悬浮铣削电主轴变质量系统的动力学模型，分析切屑从进入到离开容屑槽的整个过程，时变切屑质量对系统临界转速的影响规律；然后探索由时变切屑质量所引起的旋转惯性载荷，并以初始不平衡离心力、切削力、磁悬浮轴承电磁力以及由切屑质量引起的旋转惯性载荷等为激励，利用MATLAB求解系统的振动响应，探索由时变切屑质量所引起的系统旋转惯性载荷对变质量系统动态特性的影响规律，进而为类似质量时变系统的动力学建模与动态分析提供理论支撑与方法借鉴.

1. 磁悬浮铣削电主轴变质量系统动力学建模

以磁悬浮铣削电主轴系统的整体结构为依据，忽略系统前后保护轴承和传感器的分布位置，获得“磁悬浮轴承-主轴-刀具”系统的结构简图如图1（a）所示. 图中：1为径向磁悬浮轴承，2为电机转子，3为径向磁悬浮轴承，4为铣刀，5为切屑 6为工件.

图 1 磁悬浮铣削电主轴系统模型示意

Figure 1. System model of motorized spindle for magnetic suspension milling

下载: 全尺寸图片幻灯片

图1（a）中：F_Ax、F_Ay、F_Bx、F_By分别为上、下2个径向磁轴承在 $x$ ， $y$ 方向的电磁力；F_cx、F_cy分别为切屑旋转产生的旋转惯性力在在 $xO{\textit{z}}$ 平面和 $yO{\textit{z}}$ 平面内的投影；F_q为铣削力.

当磁悬浮铣削电主轴加工时，切屑会不断进入铣刀容屑槽，并随主轴同步旋转，由此引发2个后果：一是“主轴-刀具”系统的质量随着切屑的不断进入和离开而不断发生变化，使系统质量具有时变特性；二是随着切屑的不断进入，会使系统不平衡质量的大小和位置非线性变化，并在高速下产生旋转惯性载荷，从而引发系统新的动态特性. 为分析高速旋转下质量时变系统的非线性动态过程，探索系统动力学模型的构建，是进行动态分析的前提.

1.1 “主轴-刀具”变质量系统有限元模型的构建

采用有限元方法建立“主轴-刀具”系统质量随切屑进入和流出而不断变化的转子动力学模型. 在建模过程中，假设刀具与主轴为刚性连接的整体，将切屑视为一形状不变的刚性质量块并随铣刀同频旋转，忽略主轴切削过程中剪切变形的影响. 将“主轴-刀具-切屑”系统沿轴线划分为由“弹性轴 + 刚性圆盘”组成的等效单元，各单元之间通过节点联结，使系统简化为具有 $N$ 个节点的模型，并且将切屑等效到最后一个节点上，如图1（b）所示.

建模时，以电主轴转子的末端中心为坐标原点，主轴轴心线为z轴，如（a）所示的O-xyz坐标系；将主动磁悬浮轴承的支承特性以电磁力的形式施加在主轴转子的相应节点上；铣削产生的切削力和由切屑所导致的旋转惯性载荷作用在最后一个节点上，以及将切屑等效为在最后一个节点上的圆盘. 每个节点有 $x$ 、 $y$ 方向上的2个平移自由度和2个转动自由度，因此，对于具有 $N - 1$ 个轴段、 $N$ 个节点联结而成的“主轴-刀具”系统的x、y方向位移向量分别表示为

$\left\{ \begin{gathered} {{\bf{U}}_1} = {\left[ {{x_1},{\theta _{y1}},{x_2},{\theta _{y2}}, \cdots ,{x_N},{\theta _{yN}}} \right]^{\mathrm{T}}}, \\ {{\bf{U}}_2} = {\left[ {{y_1}, - {\theta _{x1}},{y_2}, - {\theta _{x2}}, \cdots ,{y_N}, - {\theta _{xN}}} \right]^{\mathrm{T}}}. \\ \end{gathered} \right.$

(1)

忽略系统内阻和磁悬浮轴承的阻尼，依据有限元建模理论，得到系统的运动方程为

$\left\{ \begin{gathered} {{\boldsymbol{M}}_1}{{{\boldsymbol{\ddot U}}}_1} + \Omega {{\boldsymbol{G}}_1}{{{\boldsymbol{\dot U}}}_2} + {{\boldsymbol{K}}_1}{{\boldsymbol{U}}_1} = {{\boldsymbol{Q}}_1} \\ {{\boldsymbol{M}}_1}{{{\boldsymbol{\ddot U}}}_2} - \Omega {{\boldsymbol{G}}_1}{{{\boldsymbol{\dot U}}}_1} + {{\boldsymbol{K}}_1}{{\boldsymbol{U}}_2} = {{\boldsymbol{Q}}_2} \\ \end{gathered} \right.$

(2)

式中： ${{\boldsymbol{M}}_1}$ 为质量矩阵，由于切屑质量时变导致质量矩阵 ${{\boldsymbol{M}}_1}$ 是时变的； ${{\boldsymbol{G}}_1}$ 为陀螺矩阵； ${{\boldsymbol{K}}_1}$ 为刚度矩阵，是 $2N \times 2N$ 阶对称矩阵，； $\varOmega$ 为“主轴-刀具”系统的自转角速度， ${{\boldsymbol{Q}}_1}$ 和 ${{\boldsymbol{Q}}_2}$ 均为 $2N \times 1$ 维广义力向量，包括初始不平衡力、磁悬浮轴承的电磁力、切屑旋转产生的惯性力和惯性力矩以及铣削力.

式（2）中的质量矩阵、陀螺矩阵和刚度矩阵可通过每个轴段的质量矩阵、陀螺矩阵和刚度矩阵组装起来. 每个轴段的移动惯性矩阵 ${\boldsymbol{M}}_{{\mathrm{sT}}}$ 、转动惯性矩阵 ${\boldsymbol{M}}_{{\mathrm{sR}}}$ 、陀螺矩阵 ${\boldsymbol{J}}_{{\mathrm{s}}}$ 和刚度矩阵 ${\boldsymbol{K}}_{{\mathrm{s}}}$ 表示^[18]为

$\left\{\begin{aligned} & \boldsymbol{M}_{\mathrm{sT}}=\frac{\mu l}{420}\left[\begin{array}{cccc} 156 & 22 l & 54 & -13 l \\ 22 l & 4 l^2 & 13 l & -3 l \\ 54 & -3 l & 156 & -22 l \\ -13 l & -3 l & -22 l & 4 l^2 \end{array}\right], \quad & \boldsymbol{M}_{\mathrm{sR}}=\frac{\mu r^2}{120 l}\left[\begin{array}{cccc} 36 & 3 l & -36 & 3 l \\ 3 l & 4 l^2 & -3 l & -l^2 \\ -36 & -3 l & 36 & -3 l \\ 3 l & -l^2 & -3 l & 4 l^2 \end{array}\right], \\ & \boldsymbol{J}_{\mathrm{s}}=\frac{\mu r^2}{60 l}\left[\begin{array}{cccc} 36 & 3 l & -36 & 3 l \\ 3 l & 4 l^2 & -3 l & -l^2 \\ -36 & -3 l & 36 & -3 l \\ 3 l & -l^2 & -3 l & 4 l^2 \end{array}\right], \quad & \boldsymbol{K}_{\mathrm{s}}=\frac{E I}{60 l}\left[\begin{array}{cccc} 12 & 6 l & -12 & 6 l \\ 6 l & 4 l^2 & -6 l & 2 l^2 \\ -12 & -6 l & 12 & -6 l \\ 6 l & 2 l^2 & -6 l & 4 l^2 \end{array}\right], \end{aligned}\right.$

(3)

式中： $\mu$ 为轴段单位长度的质量，l为轴段的长度；r为轴段的半径，E为轴段弹性模量，I为轴段截面的惯性矩.

以组装系统的整体质量矩阵为例进行分析. 首先将轴段的质量矩阵写成分部矩阵的形式，即

${\boldsymbol{M}}_{{\mathrm{s}}}^{(i)} = {\boldsymbol{M}}_{{{\mathrm{sT}}}}^{(i)} + {\boldsymbol{M}}_{{{\mathrm{sR}}}}^{(i)} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{m}}_{{\mathrm{11}}}^{({\mathrm{s}}i)}&{\boldsymbol{m}}_{{\mathrm{12}}}^{({\mathrm{s}}i)} \\ {\boldsymbol{m}}_{{\mathrm{21}}}^{({\mathrm{s}}i)}&{\boldsymbol{m}}_{{\mathrm{22}}}^{({\mathrm{s}}i)} \end{array}} \right].$

(4)

则磁悬浮铣削电主轴时变质量系统的整体质量矩阵M₁可表示为

${{{{\boldsymbol{M}}}}_1} = \left[ {\begin{array}{*{20}{c}} {{m}}_{{\mathrm{11}}}^{({\mathrm{s}}1)}&{{m}}_{{\mathrm{12}}}^{({\mathrm{s}}1)}&0& \cdots &0&0 \\ {{m}}_{{\mathrm{21}}}^{({\mathrm{s}}1)}&{{{m}}_{{\mathrm{22}}}^{({\mathrm{s}}1)} + {{m}}_{{\mathrm{11}}}^{({\mathrm{s}}1)}}&{{m}}_{{\mathrm{12}}}^{({\mathrm{s}}2)}& \cdots &0&0 \\ 0&{{m}}_{{\mathrm{21}}}^{({\mathrm{s}}2)}&{{{m}}_{{\mathrm{22}}}^{({\mathrm{s}}2)} + {{m}}_{{\mathrm{11}}}^{({\mathrm{s}}3)}}& \cdots &0&0 \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ 0&0&0& \cdots &{\left[ {{{{m}}_{22}}} \right]_s^{N - 2} + \left[ {{{{m}}_{{{11}}}}} \right]_s^{N - 1}}&{\left[ {{{{m}}_{12}}} \right]_s^{N - 1}} \\ 0&0&0& \cdots &{\left[ {{{{m}}_{21}}} \right]_s^{N - 1}}&{\left[ {{{{m}}_{22}}} \right]_s^{N - 1} + {{{M}}_{\mathrm{d}}}} \end{array}} \right],$

(5)

式中： ${{M}}_{{\mathrm{d}}}=\left[\begin{array}{cc}m\left(t\right)& 0\\ 0& {J}_{{\mathrm{d}}}\end{array}\right]$ ， $m\left(t\right)$ 为随时间变化的切屑质量， ${J}_{{\mathrm{d}}}$ 为时变切屑质量引起的直径转动惯量的变化量.

整体陀螺矩阵J₁和整体刚度矩阵K₁也可用类似方法形成.

1.2 铣削时电主轴变质量模型的建立

根据（a），立铣刀圆柱面上分布的切削刃呈螺旋线型，在铣削过程中，螺旋切削刃逐渐进入或退出工件，导致切削刃中切屑的厚度不断发生变化,即在切削过程中的任一时刻，沿切削刃方向分布的切屑厚度各不相同. 因此，当切屑不断进入容屑槽时，其瞬时切屑质量由切屑的瞬时切削厚度和几何边界确定；而切屑的瞬时体积为从切削刃开始接触工件的时刻 ${t_0}$ 到当前时刻 $t$ 的所有瞬时切削截面面积的积分，由此可推出单个切屑的任意时刻的质量 $\Delta m$ . 具体计算过程如下：

引用文献[19]的瞬时切削厚度模型：

$h\left(\theta \right)=\left\{\begin{array}{l}{f}_{{\mathrm{z}}}\mathrm{sin}\theta ,\quad0 < \theta < {\phi }_{{\mathrm{e}}},\\ 0,\quad 其他.\end{array}\right.$

(6)

式中： ${f_{\mathrm{z}}}$ 为每齿进给量； ${\varphi _j}\left( t \right)$ 为第 $j$ 条刀齿接触工件表面至某一瞬时时对 $t$ 时的转角，如式（7）； ${\varphi _{\mathrm{e}}}$ 为一条切削刃完成切削时的转角. 假设铣刀完全轴对称，在（b）中刀具末端单元的质心 ${o_1}$ 点建立固定坐标系 ${o_1}{x_1}{y_1}{{\textit{z}}_1}$ ，铣刀端面坐标系及相关参数如所示. 图中， ${a_{\mathrm{e}}}$ 为径向切削深度.

图 2 铣削过程示意

Figure 2. Milling process

下载: 全尺寸图片幻灯片

根据图2可得

${\varphi _{\mathrm{e}}} = {\cos ^{ - 1}}\left(\frac{{r_{\mathrm{a}} - {a_{\mathrm{e}}}}}{r_{\mathrm{a}}}\right)$

(7)

${\varphi _j}(t) = \varOmega t - (j - 1)\frac{{2{\text{π}} }}{{\textit{z}}}$

(8)

式中： $r_{\mathrm{a}}$ 为铣刀半径.

根据所示，四刃铣刀加工时，通过选择合适的铣刀规格和切削参数，可以确保每次只有一条切削刃参与切削^[20]，如所示，任一切削刃从开始切削到切削完成的时间 $\tau {\text{ = }}{{{\varphi _{\mathrm{e}}}} /\varOmega }$ ，相邻两切屑的质量变化规律在时间上相差 $\Delta \tau {\text{ = }}{ {T /k}}$ ，其中， $T$ 为主轴旋转周期， $k$ 为铣刀齿数.

图 3 切屑质量变化示意

Figure 3. Chip quality change

下载: 全尺寸图片幻灯片

由于相邻切削刃之间存在 $\Delta \tau$ 的时间间隔，任意一条切屑的质量随时间 $t$ 的变化函数为

${m}_{j}(t)=\left\{\begin{array}{l}\rho {a}_{{\mathrm{p}}}{\displaystyle {\int }_{{\phi }_{j}({t}_{0})}^{{\phi }_{j}(t)}h\left(\theta \right){\mathrm{d}}\theta },\quad{t}_{0} < t\leqslant {t}_{0} + \tau, \\ 0,\quad 其他,\end{array} \right.$

(9)

式中： $j=1、2、3、4$ ， $\rho$ 为切屑密度， ${a_{\mathrm{p}}}$ 为轴向切削深度， ${t_0} = KT + \left( {j - 1} \right)\Delta \tau$ ， $K=0、1、2\cdots$ .

当切屑 $\Delta m$ 进入容屑槽，其质心会产生一个从 ${o_1}$ 到切屑进入系统后新的质心位置 $G$ 的微小偏移量 $e$ ，以 $G$ 为坐标原点建立平行于坐标轴 ${O_1}\text{-}{x_1}{y_1}{{\textit{z}}_1}$ 的质心坐标系 $G\text{-}{x'_1}{y'_1}{{\textit{z}}'_1}$ ；此外，由于切屑的进入，原本与几何对称轴相一致的惯性主轴也会相应的产生角度为 $\varepsilon$ 的偏移，进而建立过质心的惯性主轴坐标系 $G\text{-}{x_2}{y_2}{{\textit{z}}_2}$ .

由（a）中切削区域可知，每条切削刃都是从点 $a$ 开始切削，切屑开始形成，到点 $c$ 结束切削，切屑质量达到最大. 假设切削完成时切屑即脱离容屑槽，此时切削过程中的切屑质心随着切削的进行，从点 $a$ 逐渐移动到点 $c$ . 而（a）中 $\Delta m$ 即表示切削过程中某一时刻切屑的位置；因此，切屑质心在坐标系 ${O_1}\text{-}{x_1}{y_1}{{\textit{z}}_1}$ 下的位置坐标可近似表示为

图 4 带有切屑的转动单元

Figure 4. Rotational unit with chips

下载: 全尺寸图片幻灯片

$\left\{ \begin{gathered} {x_{\mathrm{q}}} = r\cos\;\; \varOmega t \\ {y_{\mathrm{q}}} = r\sin\;\; \varOmega t \\ {{\textit{z}}_{\mathrm{q}}} = \frac{l_{\mathrm{d}}}{2} - \frac{{{a_{\mathrm{p}}}}}{\tau }t, \\ \end{gathered} \right.$

(10)

式中： $l_{\mathrm{d}}$ 为铣刀最后一个单元的长度.

由此可得到包含单个切屑的转动单元，在坐标系 ${O_1}\text{-}{x_1}{y_1}{{\textit{z}}_1}$ 下的惯性矩阵为

$\begin{split} &\boldsymbol{I}_{O_1}=\left[\begin{array}{ccc} I_{x 1} & 0 & 0 \\ 0 & I_{y 1} & 0 \\ 0 & 0 & I_{{\textit{z}} 1} \end{array}\right]+ \\ &\quad m_j\left[\begin{array}{ccc} y_{\mathrm{q}}^2+{\textit{z}}_{\mathrm{q}}^2 & -x_{\mathrm{q}} y_{\mathrm{q}} & -{\textit{z}}_{\mathrm{q}} x_{\mathrm{q}} \\ -x_{\mathrm{q}} y_{\mathrm{q}} & {\textit{z}}_{\mathrm{q}}^2+x_{\mathrm{q}}^2 & -y_{\mathrm{q}} {\textit{z}}_{\mathrm{q}} \\ -{\textit{z}}_{\mathrm{q}} x_{\mathrm{q}} & -y_{\mathrm{q}} {\textit{z}}_{\mathrm{q}} & x_{\mathrm{q}}^2+y_{\mathrm{q}}^2 \end{array}\right] \end{split},$

(11)

式中： ${I}_{x_1}、{I}_{y_1}$ 和 ${I_{{\textit{z}}_1}}$ 为未考虑切屑质量时，该单元在坐标系 ${o_1}\text{-}{x_1}{y_1}{{\textit{z}}_1}$ 下绕 $x_1、y_1、{\textit{z}}_1$ 的转动惯量， ${m_j}$ 为第 $j$ 条切削刃中的切屑质量.

由于切屑的进入，系统质心会产生一个从 ${o_1}$ 到 $G$ 的微小偏移量 $e$ ，如式（12）.

$\begin{split} &{\boldsymbol{e}} = {\left[ {{x_G}}\quad{{y_G}}\quad{{\textit{z}_G}} \right]^T} =\\ &\quad \left[ {\frac{{{m_j}}}{{{M_{\mathrm{a}}} + {m_j}}}{x_{\mathrm{q}}}}\quad{\frac{{{m_j}}}{{{M_{\mathrm{a}}} + {m_j}}}{y_{\mathrm{q}}}}\quad{\frac{{{m_j}}}{{{M_{\mathrm{a}}} + {m_j}}}{{\textit{z}}_{\mathrm{q}}}} \right] \end{split}$

(12)

式中： ${x_G}、{{y_G}}、{{\textit{z}_G}}$ 为点G坐标值， $M_{\mathrm{a}}$ 为转动单元的质量.

又因为坐标系 $G\text{-}{x'_1}{y'_1}{{\textit{z}}'_1}$ 与坐标系 ${O_1}\text{-}{x_1}{y_1}{{\textit{z}}_1}$ 平行，可根据平行轴定理，得到包含单个切屑的转动单元在坐标系 $G\text{-}{x'_1}{y'_1}{{\textit{z}}'_1}$ 下的惯性矩阵为

$\begin{split}&{{\boldsymbol{I}}_G} = {{\boldsymbol{I}}_{o_1}} - (m + {m_j}) \times \\ &\quad \left[\begin{array}{ccc} y_G^2+{\textit{z}}_G^2 & -x_G y_G & -{\textit{z}}_G x_G \\ -x_G y_G & {\textit{z}}_G^2+x_G^2 & -y_G {\textit{z}}_G \\ -{\textit{z}}_G x_G & -y_G {\textit{z}}_G & x_G^2+y_G^2 \end{array}\right]\end{split}.$

(13)

此外，由于切屑的进入，原本与几何对称轴相一致的惯性主轴也会产生倾斜，如图4（b）所示. 依据惯性矩阵的定义，总能找到一个特殊的变换矩阵 ${\boldsymbol{T}}$ 使式（13）变换为对角阵，从而获得以系统质心为原点的惯性主轴坐标系 $G\text{-}{x_2}\text{-}{y_2}{{\textit{z}}_2}$ 下的惯性矩阵：

${{\boldsymbol{T}}^{\mathrm{T}}}{{\boldsymbol{I}}_G}{\boldsymbol{T}} = \left[ {\begin{array}{*{20}{c}} {{I_x}}&0&0 \\ 0&{{I_y}}&0 \\ 0&0&{{I_{\textit{z}}}} \end{array}} \right],$

(14)

式中： ${I_x}$ 和 ${I_y}$ 为铣刀末端单元-切屑系统的x、y轴直径转动惯量， ${I_{\textit{z}}}$ 为铣刀末端单元-切屑系统的极转动惯量.

切屑随铣刀同频旋转所引起的旋转惯性力和惯性力矩可分别表示^[21]为

$\left\{ \begin{aligned} & F_{{\mathrm{c}} x}=m r \varOmega^2 \cos\;\; \varOmega t \\ & F_{{\mathrm{c}} y}=m r \varOmega^2 \sin\;\; \varOmega t \end{aligned}\right.$

(15)

$\left\{ \begin{aligned} & M_{{\mathrm{c}} x}=\left(I_x-I_{\textit{z}}\right) \varepsilon\varOmega ^2 \cos\;\; \varOmega t \\ & M_{{\mathrm{c}} y}=\left(I_y-I_{\textit{z}}\right) \varepsilon\varOmega ^2 \sin\;\; \varOmega t \end{aligned}\right.$

(16)

式中： ${F_{{\mathrm{c}}x}}$ 和 ${F_{{\mathrm{c}}y}}$ 分别为切屑运动所产生的旋转惯性力在 $xO{\textit{z}}$ 平面和 $yO{\textit{z}}$ 平面内的投影， ${M_{{\mathrm{c}}x}}$ 和 ${M_{{\mathrm{c}}y}}$ 分别为切屑运动所产生的旋转惯性力矩在 $xO{\textit{z}}$ 平面和 $yO{\textit{z}}$ 平面内的投影； $\varepsilon$ 其值可根据惯性矩阵 ${{\mathbf{I}}_G}$ 求得^[22].

1.3 铣削力模型

磁悬浮电主轴在铣削加工中会受到铣削力载荷的作用，在此引用文献[]的铣削力经验公式，可得到 $x$ 、 $y$ 方向上的铣削力均为

${F_0} = 0.85{C_{\mathrm{F}}}a_{\mathrm{e}}^{0.88}f_{\mathrm{z}}^{0.72}{a_{\mathrm{p}}}{d^{ - 0.88}}Z,$

(17)

式中： ${C_{\mathrm{F}}}$ 为铣削力系数，根据机械加工工艺设计实用手册^[24]取 ${C_{\mathrm{F}}} = 68.2$ ； $d$ 为铣刀直径.

另外，由于主轴旋转速度极快，主轴旋转一周所经历的时间极少，为简化运算，在小周期内将铣削力 ${F_0}$ 的大小视为恒定的常量，其作用时间即刀具切削刃与工件的接触时间，即产生切屑的时间，由此可建立第 $j$ 条切削刃切削时所受铣削力的模型为

${F}_{{\mathrm{q}}j}=\left\{\begin{array}{l}{F}_{0},\quad{t}_{0} < t < {t}_{0} + \tau , \\ 0,\quad 其他.\end{array}\right.$

(18)

1.4 磁悬浮轴承电磁力模型

根据磁悬浮轴承理论，并假设磁悬浮轴承各向同性，（a）中上、下两个径向磁悬浮轴承A、B在 $x$ 、 $y$ 方向产生的线性电磁力^[16]为

$\left\{ \begin{gathered} {F_{{\mathrm{A}}x}} = {k}_1{i_{ax}}{\text{ + }}{k_{\textit{z}}}{x_a} \\ {F_{{\mathrm{B}}x}} = {k}_1{i_{bx}}{\text{ + }}{k_{\textit{z}}}{x_b} \\ {F_{{\mathrm{A}}y}} = {k}_1{i_{ay}}{\text{ + }}{k_y}{y_a} \\ {F_{{\mathrm{B}}y}} = {k}_1{i_{by}}{\text{ + }}{k_y}{y_b} \\ \end{gathered} \right.$

(19)

式中： ${k_i} = \dfrac{{{\mu _0}{B^2}i_0^2A\cos {\alpha }/{2}}}{{g_0^2}}$ ， ${k_x} =k_y - \dfrac{{{\mu _0}{B^2}i_0^2A\cos {\alpha }/{2}}}{{g_0^3}}$ ，为x、y向刚度， ${i}_{{\mathrm{A}}x} 、{i}_{{\mathrm{B}}x}$ 和 ${i}_{{\mathrm{A}}y}、{i_{{\mathrm{B}}y}}$ 为磁悬浮轴承A（B） x、y的控制电流， ${\mu _0}$ 为磁导率， $B$ 为线圈匝数， ${i_0}$ 为偏置电流， $A$ 为磁路有效横截面积， ${g_0}$ 为气隙长度， $\alpha$ 为两相邻磁极之间的夹角.

2. 动态响应的计算

2.1 “主轴-刀具-切屑”质量时变系统的固有特性求解

为求解系统的固有特性，根据式（2）系统的运动微分方程，可得其齐次表达式如下：

$\left\{ \begin{gathered} {{\boldsymbol{M}}_1}{{{\boldsymbol{\ddot U}}}_1} + \varOmega {{\boldsymbol{G}}_1}{{{\boldsymbol{\dot U}}}_2} + {{\boldsymbol{K}}_1}{{\boldsymbol{U}}_1} = {\boldsymbol{0}}, \\ {{\boldsymbol{M}}_1}{{{\boldsymbol{\ddot U}}}_2} - \varOmega {{\boldsymbol{G}}_1}{{{\boldsymbol{\dot U}}}_1} + {{\boldsymbol{K}}_1}{{\boldsymbol{U}}_2} = {\boldsymbol{0}}. \\ \end{gathered} \right.$

(20)

式（19）中，用复数向量表示系统各节点的位移，即令 ${\boldsymbol{z}} = {{\boldsymbol{U}}_1} + i{{\boldsymbol{U}}_2}$ ，则式（19）可写为

${{\boldsymbol{M}}_1}{\boldsymbol{\ddot z}}{\text{ - }}{\mathrm{i}}\varOmega {{\boldsymbol{G}}_1}{\boldsymbol{\dot z}} + {{\boldsymbol{K}}_1}{\boldsymbol{z}} = {\boldsymbol{0}}$

(21)

式（20）是 $2N$ 个具有复系数的二阶微分方程组. 设其解为 ${\boldsymbol{z}} = {{\boldsymbol{z}}_0}{{\mathrm{e}}^{{\mathrm{i}}\omega t}}$ ，代入式（20）可得

$\left( { - {{\boldsymbol{M}}_1}{\omega ^2} + {{\boldsymbol{G}}_1}\varOmega \omega + {{\boldsymbol{K}}_1}} \right){{\boldsymbol{z}}_0} = {\boldsymbol{0}} ,$

(22)

式中： $\omega$ 为电主轴涡动角速度.

因此可以得到考虑陀螺力矩的频率方程为

${\text{|}} - {{\boldsymbol{M}}_1}{\omega ^2} + {{\boldsymbol{G}}_1}\varOmega \omega + {{\boldsymbol{K}}_1}{\text{|}} = 0.$

(23)

求解其特征值和特征向量，即可求出在给定自转角速度 $\varOmega$ 时，系统的涡动角速度. 为简化计算，可将式（22）改写为

${\text{|}} - \left( {{{\mathbf{M}}_1} - \frac{\varOmega }{\omega }{{\mathbf{G}}_1}} \right){\omega ^2} + {{\mathbf{K}}_1}{\text{|}} = 0.$

(24)

此时，在给定比值 $\varOmega/\omega$ 的情况下，通过求解式（22）的特征值以及特征向量就可以得到系统的模态频率和振型.

2.2 运动微分方程的求解

为方便描述，采用广义变量将式（2）合并为

${\boldsymbol{M\ddot U}} + \varOmega {\boldsymbol{G\dot U}} + {\boldsymbol{KU}} = {\boldsymbol{Q}},$

(25)

式中：质量矩阵 ${\mathbf{M}}$ 、陀螺矩阵 ${\mathbf{G}}$ 及刚度矩阵 ${\mathbf{K}}$ 为 $4N \times 4N$ 阶对称矩阵； ${\mathbf{U}}$ 为 $4N \times 1$ 维位移向量； $\boldsymbol{Q}$ 为 $4N \times 1$ 维广义力向量，其中磁轴承节点处分别为磁轴承电磁力F_Ax、F_Ay、F_Bx、F_By，末端节点处为铣削力F_q和切屑引起的旋转惯性载荷F_cx、F_cy、M_cx、M_cy，此外，其他元素为0.

由于式（24）为二阶非齐次微分方程组，求解时，首先进行降阶，然后采用四阶龙格库塔法求解该二阶微分方程组.

令 ${\boldsymbol{S}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{S}}_{\mathrm{a}}}} \\ {{{\boldsymbol{S}}_{{{\mathrm{b}}}}}} \end{array}} \right]$ ，代入式（26）整理后可得

${\boldsymbol{\dot S}} = \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\dot S}}}_{\mathrm{a}}}} \\ {{{{\boldsymbol{\dot S}}}_{\mathrm{b}}}} \end{array}} \right] = \left[ \begin{gathered} {{\boldsymbol{S}}_{\mathrm{b}}} \\ {{\boldsymbol{M}}^{ - 1}}\left[ { - \left( {\varOmega {\boldsymbol{G}}{{\boldsymbol{S}}_{\mathrm{b}}} + {\boldsymbol{K}}{{\boldsymbol{S}}_{\mathrm{a}}}} \right) + {\boldsymbol{Q}}} \right] \\ \end{gathered} \right],$

(26)

式中： ${{\boldsymbol{S}}_{\mathrm{a}}} = {\boldsymbol{U}}$ ， ${{\boldsymbol{S}}_{\mathrm{b}}} = {\boldsymbol{\dot U}}$ .

四阶龙格库塔法的计算式为

$\left\{ \begin{gathered} {{\boldsymbol{K}}_1} = F\left( {{t_n},{{\boldsymbol{S}}_n}} \right), \\ {{\boldsymbol{K}}_2} = F\left( {{t_n} + \frac{1}{2}h,{{\boldsymbol{S}}_n} + \frac{1}{2}h{{\boldsymbol{K}}_1}} \right), \\ {{\boldsymbol{K}}_3} = F\left( {{t_n} + \frac{1}{2}h,{{\boldsymbol{S}}_n} + \frac{1}{2}h{{\boldsymbol{K}}_2}} \right), \\ {{\boldsymbol{K}}_4} = F\left( {{t_n} + h,{{\boldsymbol{S}}_n} + h{{\boldsymbol{K}}_3}} \right), \\ {{\boldsymbol{S}}_{n + 1}} = {{\boldsymbol{S}}_n} + \frac{h}{6}\left( {{{\boldsymbol{K}}_1} + 2{{\boldsymbol{K}}_2} + 2{{\boldsymbol{K}}_3} + {{\boldsymbol{K}}_4}} \right). \\ \end{gathered} \right.$

(27)

式中：h为时间步长.

由式（25）可知，当已知时刻 ${t_n}$ 的 ${{\boldsymbol{S}}_n}$ 值时，给定时间步长 $h$ ，则可求出下一时刻 ${t_{n{\text{ + }}1}}$ 时的 ${{\boldsymbol{S}}_{n{\text{ + }}1}}$ 值. 因此，计算时，给定时间变量 $t$ 的范围、时间步长 $h$ 以及 ${\boldsymbol{S}}$ 的初始值，代入式（25）和式（26），即可逐步求出每个节点任意时刻的位移、转角以及速度、角速度等.

3. 动态特性分析

为分析磁悬浮铣削电主轴系统在质量非线性变化下的动态特性，根据前述磁悬浮铣削电主轴变质量系统的动力学模型，并设在额定工作转速下铣刀每齿进给量为0.06 mm ，径向切削深度为0.04 mm，轴向切削深度为1 mm，工件材料为45钢. 同时取0时刻的位移初始值 ${{\boldsymbol{S}}_0} = {\boldsymbol{0}}$ ，计算总时间设为0.5 s，仿真时利用Maltab中的ode45求解器进行，且系统具体的结构计算参数如表1所示

表 1 磁悬浮电主轴-刀具系统的物理几何参数

Table 1. Physical geometric parameters for motorized spindle of magnetic suspension–tool system

参数	数值
转子总长/mm	500
转子质量/kg	6.427
转子材料	30CrNiMo8
转子弹性模量/Pa	2 × 10¹¹
刀具总长/mm	80
刀具直径/mm	16
刀具材料	W6Mo5Cr4V2Co8
工作转速/（r•min⁻¹）	15000
磁轴承等效刚度/（N•m⁻¹）	−2 × 10⁷

下载: 导出CSV

| 显示表格

3.1 “主轴-刀具-切屑”变质量系统的固有特性分析

根据式（24）的分析和转子动力学理论，令 $\varOmega/\omega$ 为定值时，求解式（24）的特征值和特征向量，即可获得系统的模态频率和振型. 同时，为了分析系统质量变化时，转子涡动频率变化导致系统固有频率的变化趋势，以 $\varOmega/\omega$ 为横坐标，系统模态频率为纵坐标，绘制磁悬浮电主轴-刀具系统的坎贝尔图，以得到整个转速范围内系统模态频率的变化特征，如图5所示.

图 5 “主轴-刀具-切屑”变质量系统的坎贝尔图

Figure 5. Campbell diagram of “spindle–tool–chip” time-varying mass system

下载: 全尺寸图片幻灯片

由图5知，系统各阶模态频率随转子自转速度的增加而升高. 这是因为系统在切削过程中，由于切屑不断进入刀具容屑槽，不仅改变了刀具系统的质量大小和分布，同时，因切屑质量偏心所引起的陀螺力矩在正进动状态下，将对转轴的横向变形起到抑制作用，从而增强了转轴的弯、曲刚度，使得转子的正进动模态频率升高.

以图5为基础做进一步分析可知，当系统自转速度不变时，随着切屑不断进入刀具容屑槽，即系统质量增加时，系统任意阶的模态频率呈下降趋势，如图6所示. 但阶数越高，质量对系统模态频率的影响越小. 综合对比图5和图6，可以得到高转速下陀螺力矩对系统固有频率的影响要大于系统质量对固有频率的影响.

图 6 切屑质量对“主轴-刀具-切屑”变质量系统临界转速的影响

Figure 6. Influence of chip quality on critical speed of “spindle–tool–chip” time-varying mass system

下载: 全尺寸图片幻灯片

在获得“主轴-刀具-切屑”变质量系统在不同自转速度下的固有频率和变化趋势的前提下，以 $\Delta m = 2.08 \times {10^{ - 5}}\;{\text{kg}}$ 时为计算依据，求得该状态下转速比为1时系统前三阶的模态振型如图7所示.

图 7 前三阶模态振型图

Figure 7. The first three mode shapes

下载: 全尺寸图片幻灯片

由图7可知：系统前三阶振型都为弯曲振型；再结合图1（b），在磁悬浮轴承处，即节点8和节点18处，由于位移刚度的影响，系统振幅最小；在前两阶振型下，由于电机定/转子之间电磁力和电磁力矩的影响，系统振幅均接近于0；在磁悬浮轴承的两端，即系统的悬伸部分，系统不受到任何约束，类似于两悬臂梁，因此振动振幅较大，且距离轴承越远，振幅越大；在第三阶振型下，磁悬浮轴承的支承处的振幅较前两阶振型稍有增大，且两磁悬浮轴承之间的电机转子部分，出现非常明显的弯曲振动特征，容易让电机定/转子之间产生碰摩；同时，刀具部分（节点22～26）的振幅较前两阶更大. 因此，磁悬浮电主轴高速运行和进行切削加工时应当避开临界转速，防止系统失稳.

3.2 振动响应分析

在获得“磁悬浮电主轴-刀具-切屑”变质量系统固有频率、振型等固有特征随系统质量变化趋势的基础上，为进一步分析系统在运行和切削过程中的受迫振动响应特性，以系统在切削时的初始不平衡离心力、磁轴承电磁力、切削力为激励和切削过程中切屑进入容屑槽后由切屑引起的旋转惯性载荷为激励，对系统的强迫振动响应特性进行详细的仿真分析.

仿真时，令转子初始偏心量为0.2 g•mm；依据上述切削参数，获得系统的单个切屑质量 $m_j$ = 0～2.08×10⁻⁵ kg，进而可得到旋转惯性力和力矩的变化范围依次为 ${F_{{\mathrm{c}}x}} =$ 0～0.40 N， ${F_{{\mathrm{c}}y}} =$ 0～0.12 N， ${M_{{\mathrm{c}}x}} =$ 0～3.6×10⁻³ N•m， ${M_{{\mathrm{c}}y}} =$ 0～1.1×10⁻³ N•m. 最后，在Matlab中运用ode45求解器，对式（24）进行求解，得到系统的受迫振动响应如图8所示

图 8 磁悬浮电主轴-刀具-切屑系统振动响应

Figure 8. Vibration response of motorized spindle of magnetic suspension–tool–chip system

下载: 全尺寸图片幻灯片

图8（a）和图8（b）分别表示磁悬浮铣削电主轴系统各节点x、y向的振动响应；图8（c）和图8（d）分别表示系统各节点x、y向的角向振动响应. 当沿时间轴方向进行分析时，在0～0.1 s内，表示电主轴启动并旋转至工作转速. 在此期间，电主轴处于空载状态，系统只受到由初始不平衡质量引起的离心力和磁悬浮轴承电磁力的作用；从0.1 s开始，系统进入切削状态，此时由切削所引起的周期性恒定切削力和因切屑进入容屑槽产生质量偏心所引起的惯性力和惯性力矩，将使系统各节点的振动幅值迅速增大，并导致各节点平衡位置相对于中心轴有不同程度的偏移，但最后各节点的振动逐步稳定在一定范围，且呈周期性变化趋势.

当沿电主轴轴向（即节点序号）方向分析时，两磁悬浮轴承位置（节点8～18）处的振动幅值最小，其振动平衡位置基本在中心轴附近. 而在磁悬浮轴承的两端，即从节点8～1，以及从节点18～26，振动幅值越来越大；其中，节点8～1的振动平衡位置基本与中心轴一致，但节点18～26的振动平衡位置偏离中心轴的距离也越来越大，并与图7一阶模态振型的振动趋势相符. 因此，图8中系统振动的趋势与图7相同，且振动的幅值与式（15）一致，呈正弦规律周期性变化. 又因切削力施加在节点26处，故该节点及其周边处的振动幅值最大，如图8所示.

3.3 切屑引起的旋转惯性载荷对系统振动响应的影响

为进一步直观分析系统在切屑持续进入和流出容屑槽所造成的旋转惯性载荷的作用下，其振动随时间的变化规律. 选择切削加工点处以沿 $x$ 方向的径向位移和绕 $x$ 方向的角向位移为例进行分析，然后分别计算不计及旋转惯性载荷和计及旋转惯性载荷时的系统振动响应，所得到的仿真结果如图9所示.

图 9 磁悬浮电主轴-刀具-切屑系统切削加工点的振动响应

Figure 9. Vibration response of motorized spindle of magnetic suspension–tool–chip system at cutting point

下载: 全尺寸图片幻灯片

对比图9（a）、（b）可知，旋转惯性载荷的引入基本不会改变切削点处径向振动和角向振动的振动形式，但是会使径向振动响应的幅值出现0～9.7×10⁻⁷ m不等的增大；使角向振动响应的幅值出现0～2.5×10⁻⁵ rad不等的增大. 此外，旋转惯性载荷还会使加工点处径向振动和角向振动平衡位置的偏移距离分别增加约5.1～×10⁻⁷ m和9.3×10⁻⁶ rad. 因此，尽管高速铣削时由切屑引起的旋转惯性载荷不会改变系统的振动形式，但会增大系统振动的幅值，进而影响系统运行的稳定性，并降低系统切削加工的精度.

4. 结　论

1）针对磁悬浮铣削电主轴在切削过程中，因切屑不断进入容屑槽而导致“磁悬浮轴承-主轴-刀具-切屑”系统质量不断变化的情况，基于金属切削原理，建立了时变切屑质量引起的旋转惯性载荷模型，并基于转子动力学知识和有限元思想建立了磁悬浮铣削电主轴时变质量系统的动力学模型，利用Matlab软件进行求解，得到了系统固有频率的变化趋势、振型和动态响应.

2）由于陀螺力矩的影响会使系统的模态频率随着自转角速度的增大而增大；另外，切屑的进入会导致系统临界转速的下降，且对临界转速的影响随阶数增大而降低；并且，高转速下陀螺力矩对系统固有频率的影响要大于系统质量对固有频率的影响.

3）在高速切削阶段，在初始不平衡离心力、磁悬浮轴承电磁力、切削力以及切屑引起的旋转惯性载荷的共同影响下，系统振动呈现周期性；并且切屑引起的旋转惯性载荷对系统的动态特性有一定影响，特别是对刀具部分的振动影响较大，会严重影响系统的切削精度.

致谢：湘潭市科技局科技攻关项目（GX-YB20221004）的支持.

图 1 全局三维接触示意

Figure 1. Global three-dimensional contact

下载: 全尺寸图片幻灯片

图 2 刚性接触网-受电弓系统整体仿真流程

Figure 2. Simulation flowchart of pantograph and OCR system

下载: 全尺寸图片幻灯片

图 3 替换后总接触力结果对比

Figure 3. Comparison between total contact forces after replacing

下载: 全尺寸图片幻灯片

图 4 替换后仿真方法各环节耗时占比

Figure 4. Time proportion of each procedure in simulation method after replacing

下载: 全尺寸图片幻灯片

图 5 跨中点振动轨迹图

Figure 5. Vibration trajectory of mid-span point

下载: 全尺寸图片幻灯片

图 6 刚性接触网-受电弓系统快速仿真流程

Figure 6. Fast simulation flowchart of pantograph and overhead conductor rail system

下载: 全尺寸图片幻灯片

图 7 快速仿真方法与标准方法结果对比

Figure 7. Comparison of results from fast simulation method and standard method

下载: 全尺寸图片幻灯片

图 8 快速仿真方法环节耗时占比

Figure 8. Time proportion of each procedure in fast simulation method

下载: 全尺寸图片幻灯片

图 9 不同规模模型中快速仿真方法节省时间对比

Figure 9. Time saved for models in different scales by fast simulation method

下载: 全尺寸图片幻灯片

图 10 不同规模模型中结果偏差对比

Figure 10. Result error comparison of models in different scales

下载: 全尺寸图片幻灯片

图 11 不同速度下快速仿真方法节省时间对比

Figure 11. Time saved for models at different speeds by fast simulation method

下载: 全尺寸图片幻灯片

图 12 不同速度下快速仿真方法结果偏差对比

Figure 12. Result error comparison for models at different speeds by fast simulation method

下载: 全尺寸图片幻灯片

表 1 受电弓参数

Table 1. Parameters of pantograph

自由度	质量/kg	刚度/（N·m⁻¹）	阻尼/（N·s·m⁻¹）
自由度 1	7.12	9430.00	20
自由度 2	6.00	14100.00	20
自由度 3	5.80	0.01	70

下载: 导出CSV

表 2 各改进措施效果

Table 2. Effect of different improvement measures

改进方法	总接触力			仿真时间
改进方法	最大相对偏差/%	平均相对偏差/%	标准差相对偏差/%	耗时/s	节省时间比/%
标准方法				2092.89
①	3.50 × 10⁻⁴	4.80 × 10⁻⁵	5.90 × 10⁻⁵	960.13	54.12
① + ②	0.48	0.10	0.06	496.64	76.27
① + ② + ③	0.48	0.10	0.06	116.83	94.42
① + ② + ③ + ④	0.48	0.10	0.06	56.53	97.30
① + ② + ③ + ④ + ⑤	0.48	0.10	0.06	48.75	97.67

下载: 导出CSV

参考文献(34)

[1]	杜智恒. 北京大兴国际机场线160 km/h刚性接触网系统弓网耦合受流质量分析[J]. 城市轨道交通研究,2019,22(12): 157-160. DU Zhiheng. Current collection quality of pantograph-catenary coupling for 160 km/h rigid catenary system on Beijing Daxing international airport express[J]. Urban Mass Transit, 2019, 22(12): 157-160.
[2]	何国军. 160 km/h架空刚性接触网技术研究进展及展望[J]. 工程建设与设计,2022,19:131-134. HE Guojun. Research progress and prospect of 160 km/h overhead rigid catenary[J]. Construction & Design for Engineering, 2022, 19: 131-134.
[3]	梅桂明,张卫华. 刚性悬挂接触网动力学研究[J]. 铁道学报,2003,25(2): 24-29. doi: 10.3321/j.issn:1001-8360.2003.02.006 MEI Guiming, ZHANG Weihua. Study on dynamics of rigid suspension catenary[J]. Journal of the China Railway Society, 2003, 25(2): 24-29. doi: 10.3321/j.issn:1001-8360.2003.02.006
[4]	关金发,田志军,吴积钦. 基于弓网动力学仿真的160 km/h刚柔过渡系统方案研究[J]. 铁道学报,2018,40(9): 48-56. doi: 10.3969/j.issn.1001-8360.2018.09.007 GUAN Jinfa, TIAN Zhijun, WU Jiqin. Research of 160 km/h transition structure proposal between overhead conductor rail and contact line based on dynamic simulation[J]. Journal of the China Railway Society, 2018, 40(9): 48-56. doi: 10.3969/j.issn.1001-8360.2018.09.007
[5]	SONG Y, ANTUNES P, POMBO J, et al. A methodology to study high-speed pantograph-catenary interaction with realistic contact wire irregularities[J]. Mechanism and Machine Theory, 2020, 152: 103940.1-103940.18.
[6]	GREGORI S, GIL J, TUR M, et al. Analysis of the overlap section in a high-speed railway catenary by means of numerical simulations[J]. Engineering Structures, 2020, 221: 110963.1-110963.14.
[7]	周宁,邹欢,邹栋,等. 城市轨道交通弓网系统仿真模型适应性研究[J]. 西南交通大学学报,2017,52(2): 408-415,423. doi: 10.3969/j.issn.0258-2724.2017.02.026 ZHOU Ning, ZOU Huan, ZOU Dong, et al. Investigation on the applicability of pantograph and catenary model urban railway system[J]. Journal of Southwest Jiaotong University, 2017, 52(2): 408-415,423. doi: 10.3969/j.issn.0258-2724.2017.02.026
[8]	SONG Y, ZHANG M J, ØISETH O, et al. Wind deflection analysis of railway catenary under crosswind based on nonlinear finite element model and wind tunnel test[J]. Mechanism and Machine Theory, 2022, 168: 104608.1-104608.22.
[9]	吴凡平,徐钊,刘志刚,等. 山区峡谷地形风场下柔性接触网风振特性研究[J]. 铁道学报,2021,43(5): 47-61. doi: 10.3969/j.issn.1001-8360.2021.05.006 WU Fanping, XU Zhao, LIU Zhigang, et al. Research on wind-induced vibration of flexible catenary under canyon terrain wind[J]. Journal of the China Railway Society, 2021, 43(5): 47-61. doi: 10.3969/j.issn.1001-8360.2021.05.006
[10]	熊嘉铭,徐钊,鲁小兵,等. 曲线区段受电弓-接触网系统建模及动态性能分析[J]. 铁道学报,2022,44(8): 25-35. doi: 10.3969/j.issn.1001-8360.2022.08.003 XIONG Jiaming, XU Zhao, LU Xiaobing, et al. Modeling and dynamic performance analysis of pantograph-catenary system in curved railway track[J]. Journal of the China Railway Society, 2022, 44(8): 25-35. doi: 10.3969/j.issn.1001-8360.2022.08.003
[11]	BAUTISTA A, MONTESINOS J, PINTADO P. Dynamic interaction between pantograph and rigid overhead lines using a coupled FEM—multibody procedure[J]. Mechanism and Machine Theory, 2016, 97: 100-111. doi: 10.1016/j.mechmachtheory.2015.10.009
[12]	VERA C, SUAREZ B, PAULIN J, et al. Simulation model for the study of overhead rail current collector systems dynamics, focused on the design of a new conductor rail[J]. Vehicle System Dynamics, 2006, 44(8): 595-614. doi: 10.1080/00423110500165499
[13]	CHOVER J A, SUAREZ B, RODRIGUEZ P. Simulation techniques for design of overhead conductor rail lines for speeds over 140 km/h[C]//Proceedings of the 22nd International Symposium on Dynamics of Vehicles on Roads and Tracks (IAVSD 2011), Manchester: [s.n.], 2011: 14-19.
[14]	CHEN L, DUAN F C, SONG Y, et al. Assessment of dynamic interaction performance of high-speed pantograph and overhead conductor rail system[J]. IEEE Transactions on Instrumentation and Measurement, 2021, 71: 9001914.1-9001914.14.
[15]	周宁,王俊东,刘跃平,等. 基于图像处理技术检测弓网接触力的新方法[J]. 西南交通大学学报,2023,58(1): 1-8,57. doi: 10.3969/j.issn.0258-2724.20210509 ZHOU Ning, WANG Jundong, LIU Yueping, et al. Image processing based method for measuring contact force in pantograph-catenary system[J]. Journal of Southwest Jiaotong University, 2023, 58(1): 1-8,57. doi: 10.3969/j.issn.0258-2724.20210509
[16]	关金发. 受电弓与刚性接触网动力相互作用研究[D]. 成都: 西南交通大学,2016.
[17]	唐浩. AC25kV刚性接触网若干关键技术研究[D]. 成都: 西南交通大学,2018.
[18]	田珂,高仕斌,于金鑫. 基于ADAMS的地铁弓网耦合仿真分析[J]. 电气化铁道,2018,29(5): 51-54.
[19]	代洪宇. 200 km/h交流刚性接触网方案研究[D]. 成都: 西南交通大学,2019.
[20]	MASSAT J P, LAINE J P, BOBILLOT A. Pantograph-catenary dynamics simulation[J]. Vehicle System Dynamics, 2006, 44(S1): 551-559.
[21]	KULKARNI S, PAPPALARDO C M, SHABANA A A. Pantograph/catenary contact formulations[J]. Journal of Vibration and Acoustics, 2017, 139(1): 011010.1-011010.12.
[22]	CHEN L, DUAN F C, SONG Y, et al. Three-dimensional contact formulation for assessment of dynamic interaction of pantograph and overhead conductor rail system[J]. Vehicle System Dynamics, 2023, 61(9): 2432-2455. doi: 10.1080/00423114.2022.2112607
[23]	PAPPALARDO C M, PATEL M D, TINSLEY B, et al. Contact force control in multibody pantograph/catenary systems[J]. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-Body Dynamics, 2016, 230(4): 307-328. doi: 10.1177/1464419315604756
[24]	SONG Y, LIU Z G, WANG H R, et al. Nonlinear modelling of high-speed catenary based on analytical expressions of cable and truss elements[J]. Vehicle System Dynamics, 2015, 53(10): 1455-1479. doi: 10.1080/00423114.2015.1051548
[25]	COLLINA A, BRUNI S. Numerical simulation of pantograph-overhead equipment interaction[J]. Vehicle System Dynamics, 2002, 38(4): 261-291. doi: 10.1076/vesd.38.4.261.8286
[26]	CARPENTER N J, TAYLOR R L, KATONA M G. Lagrange constraints for transient finite element surface contact[J]. International Journal for Numerical Methods in Engineering, 1991, 32(1): 103-128. doi: 10.1002/nme.1620320107
[27]	ZHANG W H, MEI G M, WU X J, et al. Hybrid simulation of dynamics for the pantograph-catenary system[J]. Vehicle System Dynamics, 2002, 38(6): 393-414. doi: 10.1076/vesd.38.6.393.8347
[28]	FACCHINETTI A, GASPARETTO L, BRUNI S. Real-time catenary models for the hardware-in-the-loop simulation of the pantograph-catenary interaction[J]. Vehicle System Dynamics, 2013, 51(4): 499-516. doi: 10.1080/00423114.2012.748920
[29]	GIL J, TUR M, CORRECHER A, et al. Hardware-in-the-loop pantograph tests using analytical catenary models[J]. Vehicle System Dynamics, 2022, 60(10): 3504-3518. doi: 10.1080/00423114.2021.1962538
[30]	JIMENEZ-OCTAVIO J R, CARNICERO A, SANCHEZ-REBOLLO C, et al. A moving mesh method to deal with cable structures subjected to moving loads and its application to the catenary-pantograph dynamic interaction[J]. Journal of Sound and Vibration, 2015, 349: 216-229. doi: 10.1016/j.jsv.2015.03.051
[31]	SONG Y, LIU Z G, XU Z, et al. Developed moving mesh method for high-speed railway pantograph-catenary interaction based on nonlinear finite element procedure[J]. International Journal of Rail Transportation, 2019, 7(3): 173-190. doi: 10.1080/23248378.2018.1532330
[32]	GREGORI S, TUR M, NADAL E, et al. Fast simulation of the pantograph-catenary dynamic interaction[J]. Finite Elements in Analysis and Design, 2017, 129: 1-13. doi: 10.1016/j.finel.2017.01.007
[33]	MEI G M. Tribological performance of rigid overhead lines against pantograph sliders under DC passage[J]. Tribology International, 2020, 151: 106538.1-106538.9.
[34]	BERZERI M, SHABANA A A. Development of simple models for the elastic forces in the absolute nodal co-ordinate formulation[J]. Journal of Sound and Vibration, 2000, 235(4): 539-565. doi: 10.1006/jsvi.1999.2935