陈霁恒， 刘 星， 谭经松

(海军士官学校 机电系，安徽 蚌埠 230012)

0 引 言

薄壁件如叶轮叶片等在航空航天领域具有广泛应用，该类零件结构复杂、刚度低，刀具与工件之间的动态相互作用易产生振动失稳等问题，其高品质加工一直是制造领域的难题。由于刀具切削位置的改变及材料去除现象，工件在铣削加工中呈现结构模态时变特性，工件的实际切削量与名义值之间产生偏差，因此加工精度和表面质量难以控制。针对薄壁件加工过程结构模态的变化进行研究，分析材料去除对薄壁件模态的影响规律，对掌握薄壁件加工状态、保证其加工质量具有重要意义，为薄壁件铣削加工稳定性的深入分析，例如振动预测和颤振预报等奠定基础[1]。

在材料去除过程的动力学建模分析中，常规研究手段是利用相关有限元分析软件[2]，例如ANSYS和Abaqus等，进行模态仿真分析[3-4]；对薄壁件模型进行修改，改变其结构参数，即质量矩阵、刚度矩阵，模拟材料去除过程[5]。但在薄壁件加工的每个刀具切削位置均进行薄壁件几何模型重建与有限元模型网格重划，必然导致工作量增加。因此，为提高动力学分析计算效率，学者们从多个角度进行改进研究。

万敏等[6]使用径向切削量名义值，对原始有限元网格模型进行修正，得到工件材料去除后的动态特性，无须考虑刀具进给的中间过程。李江威[7]综合考虑刀具和工件子系统的动态特性，建立相应的铣削加工模型，并以延时微分方程描述系统工作状态。杨东辉等[8]从薄壁件加工模型理论入手，建立不同切削参数计算模型，并利用有限元软件Abaqus对模型进行仿真分析。何恩元等[2]基于三维弹性理论，提出薄板铣削加工的材料去除特性模型，并将该模型推广至整体框结构加工动态分析。MODAK[9]将模型扩展至具有弹性结构的振-声系统，利用直接矩阵修正公式，更新模型弹性结构与声学部分的质量矩阵和刚度矩阵。TUYSUZ等[10]提出一种基于改进的降阶子结构时域动态参数更新模型，降低结构矩阵维数，将该模型的预测结果与其先前建立的频域模型[11]进行比较，计算效率提升4倍。AHMADI[12]针对型腔铣削过程中的薄壁件结构动力学分析，提出一种基于有限条带建模(Finite Strip Modeling，FSM)的新方法，可同时结合多跨板(Multi-Span Plate，MSP)模型的高计算效率和有限元方法的通用性。

由于薄壁件结构模态的变化造成加工参数选择和稳态极限区域的变化，因此通过模态试验和仿真规避加工共振区域是研究结构动态特性和保证加工稳定性的主要手段之一[13-16]，除将材料去除过程划分为几段离散的加工状态进行分析外，DANG等[17]进一步提出减少质量/刚度矩阵维数与分析模态数量，通过确定主导模态，在理论上得到被加工工件的频响函数，以有效预测薄壁件铣削加工动态特性变化。TIAN等[18]将矩阵摄动法与有限元法相结合，简化求解特征值过程，通过简单的数值计算得到不同切削状态的近似模态参数。KIRAN等[19]基于逆容差耦合确定固定边界条件容差，在此基础上分析预测在铣削加工过程中材料去除时工件固有频率的变化。

薄壁件加工是在基体材料上逐步去除材料的过程，在分析材料去除过程薄壁件结构模态变化时，矩阵摄动法在原模型基础上增加摄动项并通过迭代求解，具有效率高的特点。但用于薄壁件加工过程模态分析的精度和适用性仍不明确，因此针对薄壁件加工建立模态分析的矩阵摄动方法模型，并与商业有限元软件的分析结果进行对比和验证。

1 基于摄动法的薄壁件模态参数计算

在已知每个工艺状态薄壁件的整体刚度矩阵与质量矩阵前提下，快速求解特征值和特征向量是动态模态计算的关键。但在物理空间表达的工件动力学方程的阶数较高，重复求解工件动力学方程的广义特征值存在计算量大和耗时的问题[1]。因此基于摄动理论，从初始工件结构开始，将每道工序的材料去除部分作为整个结构发生的微小改动(质量和刚度摄动)(见图1)，基于初始工件结构的模态参数，直接获取修正的模态参数，无须针对每个材料去除步骤重复求解式(1)。

图1 加工步骤i时的薄壁件结构变化示例

假设薄壁件在加工步骤i-1时的广义特征值可被描述为

Kwi-1φj,wi-1=λj,wi-1Mwi-1φj,wi-1

(1)

式中：Kwi-1和Mwi-1分别为工件的整体刚度矩阵与质量矩阵；φj,wi-1为相应的特征向量；λj,wi-1为薄壁件wi-1的j阶特征值，λj,wi-1=(ωj,wi-1)2，且ωj,wi-1为固有频率。

式(1)满足如下模态正交条件：

(2)

式中：δjk为克罗内克函数，其形式为

(3)

经1个加工步骤的材料去除，薄壁件特征值变为

Kwiφj,wi=λj,wiMwiφj,wi

(4)

但仍受到模态正交条件(式(2))的约束。由于经1个加工步骤的材料去除，系统结构参数发生小幅变动[20]，因此可将材料去除部分rsi的质量和刚度作为摄动项，则薄壁件wi的质量和刚度矩阵以摄动形式可表示为

(5)

式中：ΔKrsi和ΔMrsi为由材料去除引起的工件系统刚度和质量矩阵变化量；ε为偏差参数。

根据摄动理论，加工中的工件wi的特征值λj,wi和特征向量φj,wi可按小参数ε展开为幂级数[21]：

(6)

式中：Δλq,j,rsi和Δφq,j,rsi(q=1,2,3,…,n)分别为由于加工中的工件wi-1受到材料去除部分rsi的摄动所产生的有关j阶特征值和特征向量的q阶摄动项。由于摄动理论计算精度较高，因此可在二阶摄动项后截断式(6)，略去高阶次项O(ε3)。

将式(5)和式(6)代入式(4)，合并等式两边ε的同次幂因数，忽略高阶次项，整理可得：

ε(Κwi*1Δφ1,j,rsi+ΔKrsiΔφj,wi-1)=ε(λj,wi-1Μwi-1Δφ1,j,rsi+λj,wi-1ΔMrsiφj,wi-1+Δλ1,j,rsiMwi-1φj,wi-1)

(7)

ε2(Κwi-1Δφ2,j,rsi+ΔΚrsiΔφ1,j,rsi)=ε2(λj,wi-1Mwi-1Δφ2,j,rsi+λj,wi-1ΔMrsiΔφ1,j,rsi+Δλ1,j,rsiMwi-1Δφ1,j,rsi+

Δλ1,j,rsiΔMrsiφj,wi-1+Δλ2,j,rsiMwi-1φj,wi-1)

(8)

引入加工中的工件wi-1的未摄动特征向量，则一阶和二阶振型摄动项可表示为

(9)

式中：ajk和bjk分别为工件wi-1的k阶振型对j阶模态的一阶和二阶振型摄动项贡献因数。

将式(9)中的Δφ1,j,rsi表达式代入式(7)，分别乘以(φj,wi-1)T和(φk,wi-1)T(k≠j)，可得式(10)与式(13)：

(10)

由式(2)的模态正交条件可得：

(11)

结合式(10)与(11)，可得j阶特征值的一阶摄动公式为

Δλ1,j,rsi=(φj,wi-1)T(ΔKrsi-λj,wi-1ΔMrsi)φj,wi-1

(12)

同理，

λj,wi-1ΔMrsiφj,wi-1+Δλ1,j,rsiMwi-1φj,wi-1)

(13)

由式(2)的模态正交条件可得：

(14)

结合式(13)与(14)，可得贡献因数ajk的表达式为

(15)

类似地，将式(9)中的Δφ2,j,rsi表达式代入式(8)，分别乘以(φj,wi-1)T和(φk,wi-1)T(k≠j)，可得式(16)与式(17)：

Δλ1,j,rsiMwi-1Δφ1,j,rsi+Δλ1,j,rsiΔMrsiφj,wi-1+Δλ2,j,rsiMwi-1φj,wi-1)

(16)

Δλ1,j,rsiMwi-1Δφ1,j,rsi+Δλ1,j,rsiΔMrsiφj,wi-1+

Δλ2,j,rsiMwi-1φj,wi-1)

(17)

同理，由式(2)的模态正交条件根据式(16)可得j阶特征值的二阶摄动公式为

Δλ2,j,rsi=(φk,wi-1)T(ΔKrsi-λj,wi-1ΔMrsi-Δλ1,j,wi-1Mwi-1)Δφ1,j,rsi-(φk,wi-1)TΔλ1,j,rsiΔMrsiφj,wi-1

(18)

由式(2)的模态正交条件根据式(17)可得贡献因数bjk的表达式为

(19)

为获得j阶振型对其一阶和二阶振型摄动项的贡献因数(ajj和bjj)，取薄壁件在加工步骤i(wi)时的模态正交条件表达式为

(φk,wi)TMwiφj,wi=1

(20)

将式(5)中的质量矩阵表达式Mwi及式(6)中的j阶振型表达式φwi代入式(20)，可得：

(φj,wi-1+εΔφ1,j,rsi+ε2Δφ2,j,rsi)T(Mwi-1+εMrsi)(φj,wi-1+εΔφ1,j,rsi+ε2Δφ2,j,rsi)=1

(21)

将式(21)展开，合并ε的同次幂因数，并忽略高阶次项，整理可得：

(22)

将式(9)中的一阶和二阶振型摄动项Δφ1,j,rsi和Δφ2,j,rsi的表达式分别代入式(22)的一次和二次幂因数，并令j=k时的ε和ε2的因数为0，可得式(23)与式(25)：

(23)

由式(2)的模态正交条件根据式(23)可得：

(24)

同理，

(25)

由式(2)的模态正交条件根据式(25)可得：

(26)

将所推导的贡献因数表达式(15)、(19)、(24)和(26)代入式(9)，可得j阶振型的一阶和二阶摄动公式为

(27)

(28)

由一阶摄动公式(式(12)和式(27))与二阶摄动公式(式(18)和式(28))可知：若提取材料去除过程某个加工步骤的薄壁件动态特性，仅需要知道前一步的动态参数即可直接代入上述公式进行求解，从而避免重复求解特征方程，不仅可提升计算效率，而且可简化实际模态试验的复杂过程，便于后续的工程应用。基于该摄动理论，对上述材料去除过程进行同样的求解。

2 模态参数计算及误差分析

基于薄壁件有限元数值建模方法，组装薄壁件材料去除过程不同加工步骤的系统整体刚度和质量矩阵，并调用MATLAB中的eig函数求解广义特征值，得到其不同工序阶段的前五阶模态参数。基于上述摄动理论所推导的一阶和二阶摄动公式，同样在MATLAB中编程对材料去除过程进行求解。与Abaqus模态仿真分析所获取的固有频率进行分析对比，对摄动法求解的可靠性进行验证，可将有限元结构数值建模后的特征值求解作为精确解；将摄动公式解作为近似解，且Abaqus模态仿真分析结果较为准确，因此假定仿真解为实际值，则所设定的误差η1和η2可表示为

(29)

式中：ωFEM为薄壁件在单边固支条件下的有限元结构数值建模获得的固有频率预测值；ωFER为基于摄动理论推导的固有频率近似解；ωAbaqus为Abaqus模态仿真分析得到的固有频率。

材料去除过程薄壁件固有频率预测及误差分析如表1～表5所示。

表1 材料去除过程薄壁件固有频率预测及误差分析(初始工件)

表2 材料去除过程薄壁件固有频率预测及误差分析(工序1)

表3 材料去除过程薄壁件固有频率预测及误差分析(工序2)

表4 材料去除过程薄壁件固有频率预测及误差分析(工序3)

表5 材料去除过程薄壁件固有频率预测及误差分析(工序4)

根据表1～表5的数据，利用MATLAB绘制摄动求解误差与有限元数值建模分析误差分布图。不同工序阶段一阶和二阶摄动求解误差η1分布如图2所示。有限元数值建模分析误差η2分布如图3所示。由图3可知：一阶和二阶摄动公式对于不同工序阶段的薄壁件频率求解具有较好的近似，误差分布大约以2%作为中值。一阶摄动解已具备较好的计算精度，二阶解与一阶解相比精度提升有限，因此一阶摄动公式的求解完全可满足薄壁件材料去除过程的结构动态分析与预测。有限元数值建模求解与Abaqus模态仿真分析结果接近，尤其对于初始工件结构动态特性分析而言精度较高。观察不同工序阶段的误差随模态阶次的变化可知：随着模态阶次的提升，各工序误差变化呈现相同的趋势。

图2 不同工序阶段一阶和二阶摄动求解误差η1分布

图3 有限元数值建模分析误差η2分布

对基于一阶摄动公式求解得到的前五阶固有频率展开分析，绘制薄壁件材料去除过程不同工艺阶段的频率变化趋势。薄壁件材料去除过程模态演变趋势如图4所示。由结构动力学理论可知：模态频率变化趋势取决于模态刚度与模态质量二者谁的变化占主导地位。在模态质量下降占主导时，模态频率处于上升趋势；在模态刚度下降占主导时，模态频率处于下降趋势。对于薄壁件铣削加工过程，工件材料不断被切除，工件的模态刚度与模态质量均呈整体下降趋势。因此，由图4的频率变化趋势可知：在工件加工初期，即初始工件至工序1阶段，前三阶频率曲线呈一定上升趋势，模态质量下降占主导；随着加工的进行，各阶频率均明显下降，模态刚度下降更占优；对于频率变化幅度，在整个固有频率范围内，四阶和五阶频率变化较显著。

图4 薄壁件材料去除过程模态演变趋势

薄壁件高速铣削加工一般采用的主轴转速为3 000～12 000 r/min，若采用四齿铣刀，则对应的频率为200～800 Hz，涵盖所研究薄壁件的前二阶固有频率。因此，在实际的生产加工中应避开这些频率区域，以免发生共振，可取一阶和二阶固频中间的切削频率进行等速切削，如图4虚线①所示。为提升加工效率，可考虑采用变速切削方法，如图4虚线②所示。在加工初始阶段，基于某一转速进行铣削加工，再考虑二阶频率下降，开始提高转速，可有效远离共振频率，保证薄壁件材料去除过程的稳定性，提升加工效率。

3 结 语

基于有限元结构动力学理论，建立薄壁件加工时变结构模态分析的矩阵摄动方法，通过MATLAB建模获取薄壁件材料去除过程的整体质量矩阵和刚度矩阵，将材料去除作为当前质量矩阵和刚度矩阵的摄动项，对薄壁件模态参数变化进行解析预测，并对预测结果进行误差分析。结果表明，摄动方法分析材料去除过程工件结构模态参数变化具有较好的精度，基于该方法对薄壁件模态参数演变规律进行分析，结合模态频率变化，对薄壁件加工工艺策略进行初步规划。