具有内阻的旋转锥形复合材料刀杆铣削过程稳定性的仿真分析

任勇生，张玉环，张金峰

(山东科技大学机械电子工程学院，青岛 266590)

高速切削技术由于具有可提高加工表面质量和生产率、降低制造成本、缩短产品开发周期等优势，目前在航空航天、汽车和军工等领域得到了广泛的应用。其中，高速镗削和铣削是实现高速切削的最主要的方法[1]。高速切削系统包含机床、刀杆、工件、夹具等。从动力学的观点，刀杆作为机床最柔软的部件，其动力学特性决定着切削系统的稳定性。切削振动主要有强迫振动和自激振动[2-6]。切削过程的自激振动也称之为颤振，是由于刀具和工件之间存在交变的相互作用力而诱发产生的一种不稳定的振动[7]。

适用于深孔加工、颤振性能优异的大长径比刀杆是人们追求的目标。然而，由于材料本身的局限性，采用金属刀杆的高速铣削和镗削过程，往往由于颤振而无法进行[8-12]。动刚度和固有频率是度量刀杆在切削过程中动力学特性的两个重要的指标，不发生颤振的极限切削深度与动刚度成正比[13]。动刚度是静刚度和阻尼比的乘积，动刚度越大，颤振抑制能力越强。因此，为了提高切削系统的稳定性，刀杆应当具有较高的静刚度和阻尼。此外，增加刀杆的固有频率能够提高切削速度，从而最大限度地发挥高速切削的优越性。

采用传统的各类阻尼器增加刀杆的结构阻尼，长径比可以达到6[3,7,14]。然而，由于受几何和安装条件的限制，施加阻尼器通常不太适合于镗削系统。硬质合金刀杆的长径比虽然可以达到6 以上，但它具有较高的密度，因而固有频率增加的效果微乎其微。事实上，传统金属镗杆的动刚度和切削速度往往难以同时得到改善。智能材料刀杆虽然一定程度上增加了颤振抑制的可靠性和效果[8,15-17]，但主动与半主动控制系统结构相对复杂、设计成本昂贵，势必进一步提高深孔加工的成本。

纤维增强复合材料是由纤维和基体复合而成的一种力学性能优良的先进材料。碳纤维的模量可达933 GPa，是钢丝模量的5 倍左右[18]；聚酯树脂复合材料基体，由于高分子结构内部晶体内摩擦，显示出高阻尼特性，此外，复合材料密度低、质量轻，由此获得的高模量、高阻尼、轻质复合材料，加之力学性能的可剪裁性，这些显然是常规金属材料无法媲美的。因此，近30 年来，复合材料在高性能、高精度和高效率机床的关键部件，包括机床框架结构[19]、进给系统[20]和主轴[21]，以及刀杆[13,22-24]的研发过程中，得到了越来越广泛的应用。

借助于复合材料动刚度、比强度高的特点，可以采用更大长径比的复合材料刀杆取代金属刀杆。轻量化无疑有助于获得更高的固有频率，从而提高切削速度和加工效率。特别是，复合材料刀杆的力学特性还能通过组分材料选择、材料混杂以及优化铺层方式进行设计，因而可进一步改善深孔切削加工的能效。研究发现，两种构型的复合材料刀杆长径比可分别达到5.6 和7[13,23]，而最大切削深度是金属镗杆的5 倍[13]；采用环氧花岗岩填充的复合材料刀杆与传统金属刀杆相比，阻尼比、固有频率和加工表面粗糙度，分别提高了300%、15%和30%[25]；与硬质合金刀杆相比，复合材料刀杆的第一阶固有频率、阻尼比和动刚度，也分别提高72%、168%和28%以上，旋转复合材料刀杆的长径比可达到10.7[22]。

多年来，一些学者致力于研究复合材料刀杆切削系统稳定性问题。例如，LEE 和SUH[13]基于切削实验和冲击振动测试的方法，研究了不旋转复合材料刀杆切削系统的稳定性；NAGANO 等[23]基于有限元分析并结合切削实验，研究了不旋转复合材料刀杆的固有频率和切削极限；ZHANG 等[26]提出了一个具有约束层阻尼的复合材料刀杆切削系统理论模型，分析了切削过程的稳定性。

为了研究旋转效应对复合材料刀杆切削系统稳定性的影响，LEE 等[22]采用振动和切削实验对旋转复合材料刀杆进行了动力学性能优化；KIM等[27]研究了旋转复合材料刀杆在切削力作用下的受迫振动与颤振特性；MA 和REN[28]研究了旋转复合材料薄壁刀杆的自由振动和颤振特性。

根据转子动力学理论[29]，旋转系统的阻尼(即内阻或者旋转阻尼)不仅依赖其能量耗散特性，而且也依赖旋转速度。事实上，由于复合材料的高阻尼特性，在旋转速度超过失稳阈的条件下，复合材料内阻不但不能提高稳定性，反而会诱发颤振，削弱系统的稳定性。这与外阻(即不旋转阻尼)总是抑制颤振的作用，明显不同。因此，当刀杆的旋转速度大于失稳阈时，内阻将会对切削系统的稳定性产生负面的影响。

因此，首先有必要精确地预测复合材料刀杆的阻尼耗散特性，并以此为基础，研究材料内阻对切削稳定性的影响。然而，迄今为止，涉及内阻对旋转复合材料刀杆切削系统稳定性影响的研究极为少见，而关于复合材料刀杆阻尼建模和分析的研究也较为缺乏。文献[27 - 28]虽然考虑了陀螺效应、锥度比和内阻的影响，然而，其中内阻参数是根据经验确定的。REN 和ZHANG[30]研究了考虑内、外阻的旋转锥形刀杆的切削过程的颤振稳定性，但刀杆仅限于各向同性金属材料，并且外阻和内阻参数也是人为选定的。事实上，上述研究均缺乏有关复合材料刀杆的阻尼耗散机理的研究，此外，涉及复合材料刀杆临界速度和失稳阈等转子动力学特性也有待研究，特别是旋转复合材料刀杆的内阻失稳与切削系统颤振稳定性之间的联系，尚未被揭示和阐明。

在铣削过程中，切削力的大小和方向会随着刀杆的旋转产生周期变化，颤振数学模型是带有周期变化系数的时滞偏微分方程。按照周期系统的Floquet 定理，如果所有的特征乘数都位于复平面上的单位圆之内，则系统稳定。由于时滞系统有无穷多个特征乘数，并且特征方程不能表示为封闭形式，于是，周期时滞系统的稳定性的理论解析将变得十分困难。为此，铣削系统的稳定性分析通常采用数值的或者半解析的方法进行。KOENIGSBERGER 和TLUSTY[31]通过假定刀齿平均角以及平均切削力方向，将时滞偏微分方程中的周期系数用一个常数近似代替，在一定程度上简化了铣削系统的颤振求解过程。在半解析法中，主要应用的有时域半离散法[32-33]、时域全离散法[34]、频域零阶近似解[35]和多频解[36]。时域半、全离散法不仅能适用于时不变时滞系统，而且也适用于时变时滞系统的稳定性分析。在频域法中，刀杆的频率响应函数(FRF)可直接用于求解过程，刀杆动力学的时不变性(即具有不随时间变化的FRF)，是频域法得以应用的基本假设。研究表明，在许多实际应用中，采用频域零阶近似解，可以获得足够精确的稳定性分析结果[2]。

现有的铣削系统的稳定性研究，无论是否考虑刀杆的旋转，绝大多数是在惯性坐标系下进行的[2-3,27-28,30]。如果旋转刀杆具有轴对称特征，则无论是在惯性坐标系还是在旋转坐标系，刀杆-工件接触区的频率响应函数(FRF)均为时不变；如果刀杆为非对称，则FRF 在旋转坐标系为时不变[37-39]，而在固定坐标系为时变[37]。此时，如果依然假定FRF 为时不变，将会导致切削稳定性预测结果出现较大的误差[37-38]。

本文重点研究具有内阻的旋转锥形复合材料刀杆铣削系统的再生颤振稳定性。首先，不同于现有的文献中通常采用的基于二自由度集中质量的颤振分析模型[35-40]，从各向异性连续介质力学出发，基于复合材料粘弹性本构关系、Hamilton原理和Rayleigh 梁理论，并结合再生时滞铣削力模型，分别导出旋转坐标和惯性坐标表示的连续分布切削系统颤振偏微分方程。采用Galerkin 法对颤振偏微分方程进行离散化。其次，为了精确预测阻尼耗散特性，基于能量法计算复合材料刀杆的模态损耗因子；基于特征值分析技术计算旋转复合材料刀杆的临界速度和失稳阈；最后，分别采用时域半离散法和频域零阶近似解法预测切削系统的稳定性。分析锥度比、长径比、铺层角、铺层方式以及内、外阻对切削稳定性的影响。

1 具有内阻的旋转复合材料锥形刀杆铣削系统理论模型

1.1 运动方程

1.1.1 固定坐标表示的运动方程

图1 所示为一根转速为Ω的旋转复合材料锥形刀杆。其横截面随刀杆轴向位置x的变化按照如下规律：R(x)=[1-(1-σ)x/l]RT。式中：l是刀杆长度；σ=RR/RT表示锥度比；RT和RR分别为固定端和自由端的外径。

图1 旋转复合材料锥形刀杆结构示意图Fig.1 Schematic of the rotating tapered composite cutter bar

粘弹性阻尼耗散复合材料单层平面应力状态下的应力-应变关系为[40]：

式中：σ=[σ1σ2τ23τ13τ12]T和ε=[ε1ε2γ23γ13γ12]T分别为应力和应变； ε˙为应变对时间的一阶导数；Q为单层正轴刚度矩阵；Qd=Qη为考虑阻尼的单层刚度矩阵， η为单层阻尼矩阵，它与阻尼矩阵 ψ之间的关系如式(2)所示，阻尼矩阵 ψ表示复合材料单层的阻尼耗散特性，具体表达式为[41]：

假定单层为平面应力状态，则Q为5×5 阶方阵，其中各元素可以表示为Q11=E1/(1-ν12ν21)，Q12=Q21=ν12E2/(1-ν12ν21)，Q22=E2/(1-ν12ν21)，Q33=G23，Q44=G13，Q55=G12。其中：E1和E2为正轴弹性模量；G23、G13和G12为横向剪切模量。

式中： ψ1、 ψ2和 ψ12为与正轴薄膜影响相关的阻尼耗散系数； ψ23和 ψ13为与横向剪切影响相关的阻尼耗散参数。

本构方程式(1)是由正轴坐标(1, 2)表示的，如果在刀杆固定坐标系(x，y，z)下表示，则有偏轴刚度矩阵和偏轴阻尼刚度矩阵分别为：=T-1QT-T和=T-1QdT-T，其中，T为坐标变换矩阵：

式中， θ为纤维定向与刀杆轴线x之间的夹角。

旋转复合材料锥形刀杆的动能[30]：

式中：uy和uz分别为刀杆的沿y和z方向的横向位移； ψz和 ψy分别为刀杆的绕y和z轴的转角；m(x)和Im(x)分别为单位长度的质量和质量惯性矩。

对于如图2 所示的含M个单层的复合材料刀杆，m(x) 和Im(x) 能够表示如下[42]：

图2 层合复合材料刀杆横截面示意图Fig.2 M-layered laminate of the cutter bar

式中：rk(x)、rk+1(x)为第k层外径和内径； ρk为量密度。

复合材料锥形刀杆的势能为：

式中，随位置变化的弯曲刚度D11(x)表达如下[42]：

非保守力的虚功能够分解为外阻和内阻的虚功：

假设刀杆受到分别外阻的作用，外阻与弯曲振动速度成正比，则外阻的虚功为：

式中，CE为粘滞外阻系数。

根据文献[42]，刀杆的弯矩是由弹性应力和耗散应力共同产生：

式中：CI为Kelvin-Voigt 内阻系数；I(x)为几何惯性矩。

式(11)和式(12)中的第一项代表保守力矩，已经体现在势能表达式(7)中；其余项代表由内阻产生的非保守力的矩。

内阻的虚功表示如下：

切削力的虚功为：

式中：

式中： δ为变分符号； δD为狄拉克函数；Fy和Fz为作用于刀杆尖端的集中切削力。

基于Rayleigh 梁理论，截面转角和位移之间满足下列关系：

为了建立切削系统的运动方程，采用Hamilton原理如下：

将式(5)、式(7)、式(9)、式(14)和式(16)代入式(17)，可以导出在固定坐标系下的弯-弯耦合运动方程：

将力矩表达式(11)和式(12)代入式(18)，得：

式中：变量右上方的两撇表示对x的二阶偏导数；Ly和Lz为固定坐标系下的分布切削力。

式(19)包含2 个反对称矩阵，其中一个矩阵具有系数2Im(x)Ω，它与陀螺力有关，与旋转速度、质量惯性矩成正比，称为陀螺矩阵；另一个矩阵具有系数CII(x)Ω，称为循环矩阵，与旋转速度、内阻以及几何惯性矩成正比。显然，如果不考虑内阻，即CI=0 或者不考虑刀杆的转动惯量，即Im(x)=0，再或者当旋转速度为 0，循环矩阵矩阵将会消失，不存在内阻效应；如果转动惯量为 0，则陀螺矩阵消失，此时，不存在陀螺效应。

1.1.2 旋转坐标表示的运动方程

采用下列变换，能够将旋转坐标系转换到固定坐标系：

式中，U(t)=[uy uz]T、U(t)=[uy uz]T分别为刀杆的固定与旋转坐标系下的位移矢量。

坐标变换矩阵为：

对式(20)进行微分，得速度和加速度矢量：

式中：

将式(22)、式(23)代入式(19)，并左乘R(t)，同时考虑式(24)，可以导出在旋转坐标系下的运动方程：

式中，Ly和Lz为旋转坐标系下的分布切削力。由式(25)可见，在旋转坐标系下，陀螺矩阵与质量m(x)有关，而与质量惯性矩Im(x)无关；循环矩阵与外阻Ce有关，而与内阻CI无关。与Ω2成正比的矩阵是由离心力引起的，对旋转刀杆系统将产生负刚度。由于具有对称性，刀杆的动力学参数，如质量和质量惯性矩、内阻和外阻、刚度系数等，在旋转坐标系和在固定坐标系下是相同的。

1.2 切削力模型

1.2.1 固定坐标表示的切削力

具有旋转复合材料锥形刀杆的铣削系统横截面如图3 所示。假定刀杆具有S个刀齿。

图3 固定坐标系的切削力示意图Fig.3 Illustration of the cutting force in fixed coordinate frame

当刀杆旋转时，每个刀齿也随刀杆一起转动。将刀齿受到的再生切削力投影到固定坐标系(Y,Z)，并对此求和，得到总的再生切削力[35]：

式中：F=[Fy Fz]T；∆U=[∆Uy∆Uz]T；b为切削深度。

随时间变化的切削力方向系数为：

式中：

τ为前后两齿之间的延迟时间；Kt和Kr分别为切向和径向切削系数；ayy、ayz、azy和azz为方向系数，它们依赖于时间以及刀齿j的瞬时浸入角；ϕj=jϕP+Ωt，ϕp=2π/S为螺距角。

单位阶跃函数g(ϕj)用于确定刀齿是否处在切削状态：

式中， ϕst和 ϕex分别为刀齿的切入角和切出角。

需要说明的是，由于切削力作用于刀杆的自由端，式(19)中的矢量 ∆U具有形式：

根据频域零阶近似解法[35]，将切削力表达式(26)傅里叶级数展开并且只保留零阶项， 切削力表达式(26)可以简化为：

式中，矩阵A0是不随时间变化的，形式如下：

式(30)中的各元素为：

1.2.2 旋转坐标表示的切削力

假设刀具有N个齿，转速为Ω，旋转坐标系下的铣削过程模型图如图4 所示，

图4 旋转坐标系的切削力示意图Fig.4 Illustration of the cutting force in rotating coordinate frame

刀齿j引起的动态切削厚度可以在局部旋转坐标系中表示为[37]：

式中：c为每转一圈每齿的进给量；ϕj=ϕ-jϕpϕc为刀齿j的瞬时浸入角， ϕc为局部旋转坐标系与主旋转坐标系之间的夹角，ϕ=Ωt为刀杆的旋转角；{uyj(t)uz j(t)}T和{uy(j-1)(t-τ)uz(j-1)(t-τ)}T分别为当前刀齿和前一个刀齿的振动位移。

刀齿j在径向方向和切向方向上的切削力可以表示为：

将式(32)代入式(33)，可得：

式中，Tj为齿j的正交坐标变换矩阵，表示局部旋转坐标系与整体旋转坐标系之间的关系为：

1.3 模型求解

1.3.1 Galerkin 法

为了求解偏微分方程式(19)，首先采用Galerkin方法对其进行离散化。

假设：

式中：Uy j(t)和Uz j(t)(j=1, 2, ···,n)为广义坐标；ψj(x)为振型函数。旋转Rayleigh 梁的振型函数受转速和转动惯量参数的影响，其计算涉及非线性特征方程的求解，过程复杂[43]；而不旋转Euler-Bernoulli 梁振型函数与转速无关，计算过程简单。研究表明[43]：在Galerkin 近似计算中，采用不旋转Euler-Bernoulli 梁振型函数能够得到与旋转Rayleigh 梁的振型函数非常接近的计算结果。

因此，为了简化Galerkin 求解过程，选取悬臂Euler-Bernouli 梁的振型函数，ψj(x)(j=1, 2, ···,n)，它满足边界条件：

将式(38)代入式(19)，利用Galerkin 法，并且考虑边界条件式(39)，可以得到下列方程：

其中：

在式(41)中，矩阵中的各项元素的计算表达式分别为：

从式(42)中可以看出，阻尼矩阵C中包含陀螺效应和内阻和外阻，陀螺效应只与阻尼矩阵相关。刚度中的交叉耦合刚度系数可能会导致内阻失稳。

1.3.2 基于能量法的阻尼预测

为了建立锥形复合材料刀杆的内阻预测模型，需要研究如何求解得到振动方程中的内阻尼系数CI和外阻尼系数CE。

单自由度系统阻尼比的定义为[30]：

式中，ζ为切削系统的总阻尼比，可以表示成内阻尼比与外阻尼比之和，即ζ=ζE+ζI。

内阻尼项为：

当给定了系统的总阻尼比ζ 和内阻尼比 ζI，通过式(43)和式(44)，可求得内阻项和外阻项，再由式(42)中的第四、第五式，可求得内阻尼系数CI和外阻尼系数CE。模型中的内阻尼比ζI，可以通过锥形复合材料刀杆的损耗因子 ηI来求解得到。

对于弱阻尼情况，模态阻尼比与损耗因子之间存在下列关系[44]：

式中， ηI为锥形复合材料刀杆的损耗因子，其数值可以通过能量法求解得到[45]：

式中：U为不旋转复合材料刀杆的无阻尼挠度；和分别为不旋转锥形复合材料刀杆的刚度矩阵和阻尼矩阵，其计算表达式分别为：

其中，刚度矩阵元素计算表达式与方程组式(42)中的第三式相同。

阻尼矩阵的元素计算表达式为：

1.3.3 复合材料刀杆转子动力学

为了对旋转复合材料锥形刀杆的转子动力学特性进行分析，令式(40)右端的切削力为0，由此导出自由振动方程，其自由振动解可其转化为求解下列特征问题：

式中，0 和I分别为零矩阵和单位矩阵。

特征方程式(49)的根可以表示成λk=σk+iωdk,(k=1,2,3,···)，实部 σk为系统的衰减率，虚部ωdk为系统的涡动频率，它们与转速有关。如果衰减率是负的，那么该系统是稳定的；如果衰减率是正的，那么该系统是不稳定的。

1.3.4 铣削过程稳定性

1.3.4.1 频域零阶近似解

对固定坐标系下的离散颤振方程式(40)进行拉氏变换，令s=λ+iω，s的取值与铣削系统的稳定性存在以下关系：如果λ>0，那么铣削系统是不稳定的；如果λ<0，那么铣削系统是稳定的；如果λ=0，那么铣削系统处在稳定与不稳定之间的临界状态。

令λ=0，将s=iωc代入式(40)拉氏变换后的方程中，可以得到：

式中，ωc为颤振频率。

铣削系统的稳定性问题等价于下列特征值问题：

式中，I为单位矩阵。

复合材料刀杆传递函数：

铣削过程的极限切削深度可以通过(52)中的第一式求得如下：

将复特征值Λ=ΛR+iΛI以及欧拉公式e-iωcT=cos(ωcT)-isin(ωcT)代入式(54)，并且令其虚部等于0，可得：

由式(55)，可得：

注意，κ=ΛI/ΛR=arctanψ ， ψ为相位角。

由式(56)，进一步可得：

式中，ε=π-2ψ表示当前刀齿和前一个刀齿振痕间的相移。于是，旋转速度为Ω=60/S T。

最后，将式(56)代入式(54)的实部，可得极限切削深度：

由于旋转复合材料刀杆的传递函数是随转速变化的，ALTINTAS 和BUDAK[35]提出的稳定性极限图计算方法和步骤不能直接应用。为此，将采用如下迭代过程计算稳定性极限：

1) 根据刀具的材料参数和几何参数计算方向系数矩阵。

2) 给定转速Ω，并开始循环。

3) 计算刀杆的一阶固有频率。

4) 参考计算得到的一阶固有频率，设定扫描频率的范围，计算刀杆结构的传递函数。

5) 根据特征方程式(51)求解特征根的实部和虚部，将特征根代入式(58)，得到极限切削深度及相对应的转速。

6) 如果步骤5)中的转速与给定转速的差值满足给定的精度要求，那么以转速为横轴、以极限切削深度为纵轴绘制稳定性极限图；如果两者的差值不满足精度要求，继续循环步骤2)～步骤6)，直到两者的差值达到要求的精度。

7) 选择新的j并计算相邻的叶瓣。

1.3.4.2 时域半离散解

将采用时域半离散法，对旋转坐标系下的铣削系统(式(25))进行颤振稳定性分析。对于固定坐标下的铣削系统(式(19))而言，时域半离散法的求解过程可以采用完全类似的计算格式和计算步骤。

采用Galerkin 法对方程式(25)进行离散化，得：

切削力表达式为：

式中，矩阵Z中的元素表示为Zk f=ψk(l)ψf(l),(k,f=1,2,···,n)。

根据柯西变换，颤振方程式(59)可以用下列一阶微分方程表示：

式中：

将延迟时间τ划分成n个时间间隔 ∆t，q(ti-τ)可近似为：

因此，动力系统(式(62))在离散时域内有下列形式：

在时间间隔 ∆t内，微分方程式(65)的解可以表示为齐次通解和特解之和[2]：

C0依赖于初始条件，其值可通过在t=ti+1时求解式(61)得到：

当t=ti+1时，

从式(68)可以看出，求解qi+1时，需要用到qi、qi-n+1和qi-n。

离散时域内的离散映射可以表示成：

式中：

在时间间隔Δt内，通过求解式(69)的离散解集，可以模拟时变铣削过程。

由于铣削过程在刀齿周期内是周期性变化的，所以，只需要求解m个时间间隔内的微分方程就足够了。

在m个时间间隔内运用方程式(69)，可以计算铣削系统的稳定性

根据 Floquet 理论，铣削系统的稳定性可以通过矩阵T的特征值进行评价[32]：如果矩阵T的特征值小于1，系统稳定，如果矩阵T的特征值等于1，系统处于临界状态，如果矩阵T的特征值大于1，系统不稳定。

时域半离散法计算复合材料刀杆铣削系统稳定性极限图的主要骤如下：

1) 定义齿数、切向切削力系数、径向切削力系数等切削参数，根据刀杆的性能给出y方向和z方向上刀杆的质量、阻尼、刚度等参数；

2) 给定系统的起始转速与结束转速，设定切削深度的扫描范围；

3) 定义离散时间间隔的个数与步长，分别计算y方向和z方向上的动态切削力；

4) 设定转速与切削深度的扫描步长，并开始循环；

5) 计算传递矩阵T；

6) 计算传递矩阵的特征值，并分别保存使得矩阵特征值小于1 的转速值和切削深度值；

7) 循环计算，直至完成给定的转速范围和切削深度范围内的计算。

将步骤6)中满足条件的转速值和切削深度值分别取出，以转速为横轴、以切削深度为纵轴便可以绘制出稳定性极限图。

2 结果与讨论

2.1 阻尼耗散因子

为了检验本文复合材料刀杆内阻模型和求解方法的正确性，下面针对圆形空心截面复合材料悬臂梁，其平均直径d=352 mm，厚度h=10.16 mm，铺层数为16 层，复合材料的性能参数，见文献[45]。

表1 给出对应于不同铺层方式和长径比的模态损耗因子，并且与文献[45]中的结果进行比较。

表1 复合材料圆形截面梁的模态损耗因子Table 1 Modal loss factors of circular composite beam

从表1 中可以看出，采用本文模型计算得到的不同长径比、铺层方式的模态损耗因子的数值结果，与文献[45]的有限元计算数值结果较为吻合，说明本文阻尼模型与求解方法是正确的。

在下面的算例中对复合材料刀杆的损耗因子进行数值计算，刀杆的材料为碳纤维环氧树脂复合材料，其材料性能参数见表2，刀杆的外圆直径为0.06 m、内圆直径为0.04 m，铺层方式为[±θ]8，计算时选取的振型函数个数为3。

表2 碳纤维环氧树脂复合材料性能参数[41]Table 2 Mechanical properties of carbon/epoxy composite material[41]

图5 表示复合材料刀杆的模态损耗因子随铺层角的变化图，从图中看出，当铺层角从0°增加到81°，损耗因子随之增加。当铺层角等于0°时，损耗因子最小(约为0.45%)；当铺层角等于81°时，损耗因子最大(约为4.28%)；当铺层角从81°增加到90°时，损耗因子略微减小。铺层角较大时的耗散因子较大的原因，是纤维横向阻尼能力大于其纵向阻尼能力。

图5 模态损耗因子随铺层角变化图(l/d=10，σ=1)Fig.5 Variation of modal loss factor with ply angle(l/d=10，σ=1)

图6 表示模态损耗因子随长径比的变化图，从图中可以清晰地看出，长径比对损耗因子的影响非常小，可以忽略不计，这与文献[45]中的结论是一致的。图7 表示模态损耗因子随锥度比的变化图，从图中可以看出，锥度比的变化对损耗因子的影响很小。

图6 模态损耗因子随长径比变化图(σ=1)Fig.6 Variation of modal loss factor with aspect ratio(σ=1)

图7 模态损耗因子随锥度比变化图Fig.7 Variation of modal loss factor with taper ratio

2.2 复合材料刀杆试件的锤击实验

实验测试平台以及测试现场如图8 所示。试件锤击模态测试系统包括LC02 型号力锤、3A102型号力传感器、DH187E 型号IEPE 压电式加速度传感器、DH5857 电荷适调器以及DH5923N 动态信号测试分析系统等。试件根部加持固定在车床卡盘上，试件长度450 mm、外径22 mm、内径6 mm，铺层方式 [0°/±45°]7，由碳纤维HR40 与环氧树脂采用缠绕成型工艺制作而成，纤维含量46%。试件组分的力学性能参数，分别如表3 和表4 所示。

表3 碳纤维HR40 的材料性能参数Table 3 Mechanical properties of carbon fiber HR40

表4 环氧树脂的材料性能参数Table 4 Mechanical properties of epoxy resin

图8 实验测试台及测试现场Fig.8 Experiment test bench and experiment test site

表5 给出刀杆前三阶固有频率和阻尼比的实验和模型分析结果。当悬伸长度增加时，刀杆固有频率随之减小，而阻尼比几乎不受悬伸长度变化的影响，从理论和实验数据中反映出的上述变化趋势是相同的。但是，由于刀杆的根部与卡盘间存在间隙，可能会影响实验结果。此外，试件加工误差也会影响其力学性能参数的准确性，从而导致理论与实验之间的偏差。

表5 不同悬伸长度下刀杆1 的前三阶固有频率与阻尼比Table 5 The first three natural frequencies and damping ratios of cutter bar 1 under different overhang lengths

2.3 临界转速与失稳阈

利用旋转复合材料刀杆的特征方程，可以计算旋转复合材料刀杆的涡动频率、衰减率、临界转速和失稳阈等振动特性(如不特殊指明，刀杆的几何尺寸和材料性能参数都与2.2 节相同)。

表6 表示铺层角对不旋转复合材料刀杆的前三阶固有频率和衰减率的影响，其中，分别列出采用本文模型以及ANSYS 有限元软件的计算结果；利用方程式(49)计算时选取的振型函数个数为3。图9 给出锥形复合材料刀杆的ANSYS 三维有限元模型，采用六面体扫描法对模型进行网格划分，共生成8032 个单元，设置瑞利阻尼，通过QR阻尼特征值提取法可以得到模型的模态分析结果。

表6 不同铺层角的不旋转复合材料刀杆的前三阶固有频率与衰减率(ζ=0.01, l/d=10, σ=1)Table 6 The first three natural frequencies and decay rates of the non-rotating composite cutter bar with different ply angles(ζ=0.01, l/d=10, σ=1)

图9 锥形复合材料刀杆的三维有限元结构模型Fig.9 3D FEM structural model of the tapered composite cutter bar

从表6 中可以看出，随着铺层角的增加，刀杆的前三阶固有频率减小，衰减率的绝对值增加，本文模型结果与ANSYS 分析结果之间较为接近。

表7 给出了具有不同锥度比的旋转复合材料刀杆在不同转速下的第一阶正进动涡动频率(FW1)和第一阶反进动涡动频率(BW1)。

表7 不同转速下旋转锥形复合材料刀杆的FW1 和BW1(ζ=0.01, θ=30°, l/d=10)Table 7 The first forward and backward whirl frequenciesof the tapered rotating composite cutter bar with different rotating speeds (ζ=0.01, θ=30°, l/d=10)

从表7 中可以看出，采用本文提出的求解方法得到FW1 和BW1 与采用ANSYS 有限元仿真得到的结果较为吻合。

旋转复合材料刀杆的临界转速可通过直线和正进动涡动曲线的交点获得；正进动衰减率由负变为正的转速代表失稳阈。

表8 表示铺层角对临界转速和失稳阈的影响(铺层方式[±θ]8)。结果表明：临界转速和失稳阈随着铺层角的增加而减小。此外，从表8 还可以看出，在考虑内外阻的情况下，失稳阈总是大于临界转速。

表8 铺层角对临界转速和失稳阈的影响(ζ=0.01, l/d=10, σ=1)Table 8 Effect of ply angle on critical rotating speed and instability threshold (ζ=0.01, l/d=10, σ=1)

表9 表示3 种不同的对称铺层对复合材料刀杆临界转速和失稳阈的影响。对称铺层方式包含6 个单层，具有0°和90°铺层角，结果表明，0°的层数越多，临界转速越大，失稳阈也越大，这是由于纤维角为0°的层数越多，刀杆的纵向刚度越大。

表9 对称铺层方式对临界转速和失稳阈的影响(ζ=0.01, l/d=10, σ=1)Table 9 Effect of symmetric laminate on critical rotating speed and instability threshold (ζ=0.01, l/d=10, σ=1)

表10 表示5 种不同的非对称铺层对临界转速和失稳阈的影响。铺层方式包含10 层，从刀杆内部开始铺设，其中0°铺层数为6，90°铺层数为2，45°和-45°铺层数均为1。结果表明，0°铺层越接近复合材料刀杆的外表面，临界转速和失稳阈越大。由于靠近复合材料刀杆的外表面的层具有更大的周长和体积，它们承受更多的弹性应力和耗散应力产生的弯矩，因此，外表面层控制着复合材料刀杆的刚度以及阻尼。与其他铺层相比，0°铺层具有更高的刚度和更小的阻尼，因此将0°铺层设置在复合材料刀杆的外表面，可以提高刚度和稳定性。

表10 非对称铺层方式对临界转速和失稳阈的影响(ζ=0.01, l/d=10, σ=1)Table 10 Effect of asymmetric laminate on critical rotating speed and instability threshold (ζ=0.01, l/d=10, σ=1)

2.4 稳定性叶瓣曲线

为了对旋转复合材料刀杆铣削系统的稳定性进行数值计算，考虑一个半侵入逆铣加工过程，对应的切削参数均取自文献[35]。铣削力刚度为Kt=1500MPa，Kr=450MPa，刀齿数N=8，切入角和切出角ϕst=0◦，ϕex=90◦。

2.4.1 模型验证

首先要验证本文提出的铣削系统颤振稳定性时域半离散求解方法的正确性，数值算例模型是旋转坐标系下的铣削系统模型，刀杆结构参数与切削参数均与文献[37]中的一致：

图10 为旋转坐标下铣削系统颤振稳定性叶瓣图，由图10 可以看出，采用本文模型得到的叶瓣图，位于文献[37]的叶瓣图的上方，这是由于本文模型同时考虑了内、外阻的影响，而文献[37]的模型中仅包含外阻。如果在本文模型中不考虑内阻，那么，采用本文的叶瓣图，将退化为文献[37]的叶瓣图。

图10 旋转坐标系下的颤振稳定性叶瓣图Fig.10 Chatter stability diagrams in rotating coordinate frame

2.4.2 基于不同方法的稳定性预测结果对比

图11 表示采用频域零阶近似法和时域半离散法获得的铣削系统稳定性叶瓣图。结果表明：基于两种方法的叶瓣图是基本重合的。此外，从图11还可以看到，除了位于低转速区的稳定性叶瓣图之外，在高转速区，旋转刀杆铣削系统出现了新的不稳定区域，这是由于考虑了复合材料的阻尼、刀杆旋转以及再生切削力综合作用的结果。在随后的计算过程中，将对这种新的颤振不稳定现象进行分析与讨论。

图11 不同铺层方式下频域零阶近似法与半离散结果比较(ζ=0.01, σ=0.75, l/d=15, [90°]6)Fig.11 Comparison between zero-order frequency domain method and semi-discretization method under different laminates (ζ=0.01, σ=0.75, l/d=15, [90°]6)

需要指出的是，在对颤振方程进行Galerkin离散化，绘制稳定性叶瓣图的过程中，分别取刀杆的振型函数个数n=1, 2, 3，并将计算结果进行了比较，发现n=2 和3 的叶瓣图基本重合，这说明取3 项振型函数即可得到收敛的结果。因此，在随后的计算中，如不加以说明，振型函数的个数均取n=3。

2.4.3 基于不同坐标系的稳定性预测结果对比

图12 给出了固定坐标系和旋转坐标系下，铣削系统的稳定性叶瓣图。可以看出，这些叶瓣图是重合的。

图12 对称铺层方式[0°]6 对两种坐标系下铣削系统稳定性的影响(ζ=0.01, σ=0.75, l/d=10)Fig.12 Effect of symmetric laminate [0°]6 on the stability of milling system in two coordinate frames(ζ=0.01, σ=0.75, l/d=10)

文献[37]和文献[38]基于简单的二自由度铣削动力学模型，研究了铣削系统在固定坐标和旋转坐标系下的颤振稳定性，结果表明，当刀杆具有对称动力学特性时(在两个自由度方向上的集中质量、刚度和阻尼参数分别相等)，在两个坐标系下的叶瓣图是一致的。本文从更为复杂的复合材料刀杆连续分布的铣削动力学模型出发，预测了铣削系统的颤振稳定性，研究发现，具有轴对称刀杆的铣削系统切削的稳定性分析结果，不会因动力学模型参照系选取的不同而有所不同。因此，相对于文献[37 - 38]而言，本文的研究结果显然更具有一般性。

2.4.4 颤振稳定性的参数影响分析

图13(a)表示旋转锥形复合材料刀杆铣削系统的稳定性叶瓣图，图13(b)表示旋转锥形复合材料刀杆的衰减率曲线，图14 表示与图13(a)中的4 个参数点对应的时间响应曲线。这些响应曲线是采用4 阶龙格-库塔法对方程式(40)进行数值积分获得。表11 给出了四个点的转速值和切削深度。在图13(a)中，深色区域代表不稳定，其它区域代表稳定，因此，所选取的点1 和点3 对应稳定切削状态，点2 和点4 对应不稳定切削状态。从图14中可以看到，点1 和点3 的时间响应曲线是收敛的，说明这两个点对应的是稳定切削状态，点2 和点4 的时间响应曲线是发散的，说明这两个点对应的是不稳定切削状态，这与图13(a)中的预测结果是一致的。图13(a)显示的新的颤振不稳定启动速度等于64 890 rad/s，这与旋转锥形复合材料刀杆转子系统的失稳阈是相等的，如图13(b)所示。

表11 图13(a)中参数点对应的切削参数Table 11 Process parameters according to the points in Fig.13(a)

图13 旋转锥形复合材料刀杆铣削系统稳定性叶瓣图与衰减率曲线(ζ=0.01, θ=0°, σ=0.75, l/d=10)Fig.13 Stability lobe diagram and decay rate curve of the rotating tapered composite cutter bar milling system (ζ=0.01,θ=0°, σ=0.75, l/d=10)

图14 时间响应曲线Fig.14 Time domain response curves

2.4.4.1 铺层角的影响

图15 表示铺层角对铣削系统稳定性的影响。从图15 中可以看出，随着铺层角的增加，新的颤振不稳定启动速度随着铺层角的增加而减小，这是由于旋转复合材料刀杆的内阻会随着铺层角的增加而增加。

图15 铺层角对铣削系统稳定性的影响(ζ=0.01, σ=0.75, l/d=10)Fig.15 Effect of ply angle on the stability of milling system(ζ=0.01, σ=0.75, l/d=10)

2.4.4.2 长径比的影响

图16 表示长径比对铣削系统稳定性区域的影响。由图16 可以看出，随着长径比的增加，新的颤振不稳定启动速度减小，因此，细长刀杆相对更容易引起切削颤振。

图16 长径比对铣削系统稳定性的影响(ζ=0.01, σ=0.75, θ=30°)Fig.16 Effect of aspect ratio on the stability of milling system (ζ=0.01, σ=0.75, θ=30°)

2.4.4.3 内阻与外阻的影响

图17 表示内、外阻对铣削系统稳定性的影响。图17 中的三种情形分别表示只考虑内阻(图17(a))、考同时虑内、外阻(图17(b))和只考虑外阻(图17(c))。结果表明：无论是否考虑外阻，只要考虑将内阻的影响，高转速区必定出现新的不稳定区域，不稳定启动速度与旋转复合材料刀杆的失稳阈相同；如果只考虑外阻，由于失稳阈为无穷大，在高转速区将不会出现新的不稳定区域；如果不考虑外阻，不稳定启动速度则为临界转速(远低于考虑外阻时的失稳阈)，这意味着铣削系统的稳定性将会减弱。综上所述，如果忽略内阻的影响，将会高估颤振稳定性；反之，如果忽略外阻的影响，将会低估颤振稳定性。内阻在低转速区起到提高铣削系统稳定性的作用，但在高转速区会引起新的不稳定现象。

图17 内阻和外阻对铣削系统稳定性的影响(σ=0.75, θ=90°, l/d=10)Fig.17 Effect of internal and external damping on the stability of milling system (σ=0.75, θ=90°, l/d=10)

2.4.4.4 锥度比的影响

图18 表示锥度比对铣削系统稳定性区域的影响，其中包含锥形复合材料刀杆的体积与等截面复合材料刀杆的体积相等(图18(a))和体积不相等(图18(b))两种情形。图18(a)表明：在低转速区，随着锥度比的减小，稳定性叶瓣图的位置升高；在高转速区，新不稳定区域开始出现的起始转速值增大，铣削系统的稳定性提高。图18(b)中的结果是与图18(a)中的结果刚好相反的。说明，只有当锥形与等截面刀杆的体积相等时，减小锥度比，才能提高系统的切削稳定性。

图18 锥度比对铣削系统稳定性的影响(ζ=0.01, θ=90°, l/d=10)Fig.18 Effect of taper ratio on the stability of milling system(ζ=0.01, θ=90°, l/d=10)

2.4.4.5 铺层方式的影响

图19 表示对称铺层方式对铣削系统稳定性的影响。3 种对称铺层方式均为6 层，两种铺层角为0°和90°。在低转速区，0°的铺层越多，叶瓣图的位置越高，稳定性越好，这是由于，0°铺层越多，刀杆的刚度越大；在高转速区，0°的铺层越多，颤振不稳定启动速度越大，这是因为，0°的铺层越多，刀杆的阻尼越小。

图19 对称铺层方式对铣削系统稳定性的影响(ζ=0.01, σ=0.75, l/d=10)Fig.19 Effect of different symmetric laminates on the stability of milling system (ζ=0.01, σ=0.75, l/d=10)

图20 表示非对称铺层方式对铣削系统稳定性的影响。五种铺层方式均包含10 层，从内部开始铺设，其中，0°铺层共6 层，90°铺层共2 层，45°和-45°铺层各1 层。由于靠近外表面的铺层具有更大的周长和体积，它们承受更多的弹性应力和耗散应力产生的弯矩，所以，外表面层控制着刀杆的刚度及阻尼。与其他角度的铺层相比，0°铺层具有更高的刚度和更小的阻尼，因此，将0°铺层放在外表面的位置，可以提高刀杆的刚度和稳定性。此外，从图20 还可以看出，五种不同的铺层方式中，在低转速区，铺层方式[90°/45°/-45°/90°/0°6]对应的稳定性叶瓣图位置最高，铺层方式[0°6/90°/45°/-45°/90°]对应的稳定性叶瓣图位置最低；在高转速区，铺层方式[90°/45°/-45°/90°/0°6]对应的不稳定启动速度最大，系统稳定性最好。

图20 非对称铺层方式对铣削系统稳定性的影响(ζ=0.01, σ=0.75, l/d=10)Fig.20 Effect of different asymmetric laminates on the stability of milling system (ζ=0.01, σ=0.75, l/d=10)

表12 和表13 分别表示对称铺层复合材料刀杆铣削系统的颤振启动速度随长径比和锥度比的变化规律。表14 和表15 分别表示长径比和锥度比对非对称铺层复合材料刀杆铣削系统的颤振启动速度的影响。

表12 长径比对三种对称铺层方式的不稳定启动速度的影响(ζ=0.01, σ=0.75)Table 12 Onset rotating speeds of new unstable region for three symmetric laminates with variable aspect ratios(ζ=0.01, σ=0.75)

表13 锥度比对三种对称铺层方式的不稳定启动速度的影响(ζ=0.01, l/d=10)Table 13 The onset rotating speeds of new unstable region for three symmetric laminates with variable taper ratios(ζ=0.01, l/d=10)

表14 长径比对五种非对称铺层方式的不稳定启动速度的影响(ζ=0.01, σ=0.75)Table 14 The onset rotating speeds of new unstable region for five asymmetric laminates with variable aspect ratios(ζ=0.01, σ=0.75)

表15 锥度比对五种非对称铺层方式的不稳定启动速度的影响(ζ=0.01, l/d=10)Table 15 The onset rotating speeds of new unstable region for five asymmetric laminates with variable taper ratios(ζ=0.01, l/d=10)

表12 和表13 表明，随着长径比/锥度比的增加，不同对称铺层方式的不稳定启动速度之间的差基本不变。例如，当长径比分别为10、12、15 时，铺层方式A 与C 之间的差分别是179%、180%、180%；当锥度比分别为0.5、0.75 和1 时，铺层方式A 与C 之间的差分别是178%、179%、181%。但是，从表14 和表15 中可以看出，随着长径/锥度比的增加，不同非对称铺层方式之间的差是变化的。

3 结论

由于复合材料的高阻尼特性，当转速大于失稳阈时，旋转复合材料轴转子系统将会发生内阻失稳。为了揭示内阻对高速旋转复合材料刀杆铣削系统颤振稳定性的影响，本文分别在固定和旋转坐标系下，推导出具有旋转锥形复合材料刀杆的铣削系统的偏微分运动模型。基于复合材料粘弹性本构关系并结合能量法，对复合材料刀杆的阻尼耗散特性进行理论分析。采用Galerkin 法对偏微分方程进行离散化。采用时域半离散法和频域零阶近似解法，预测铣削系统的颤振稳定性。研究结果如下：

(1) 铺层角和铺层方式影响阻尼耗散因子从而对铣削系统的内阻产生重要影响；只有在考虑旋转的情况下，才能捕捉内阻对铣削系统稳定性的影响，如果不考虑旋转，颤振只存在于低速区，而在高速区是稳定的；在同时考虑旋转和内阻的情况下，在高速区将出现一个新的不稳定区域，新的不稳定区域的起始转速等于旋转复合材料刀杆的失稳阈。

(2) 颤振起始速度随着长径比、铺层角或者锥度比的增加而减小；在刀杆外表面一侧尽可能多地铺设具有较高轴向刚度和较低轴向阻尼的层，能够提高铣削系统的稳定性。

(3) 在固定和旋转坐标系下，轴对称复合材料刀杆铣削系统颤振稳定性预测结果是相同的。

(4) 基于Galerkin 法的稳定性叶瓣图振型截断近似计算过程，具有较好的收敛性；基于时域半离散法和频域零阶近似解法的稳定性叶瓣图的一致性；极限切削深度-转速参数平面上点的时间响应特征(收敛、发散)与稳定性叶瓣图的预测结果具有一致性。

(5) 加工制造复合材料刀杆时，将更小角度的铺层铺设在刀杆外表面，有助于提高刀杆的刚度和稳定性。