薄壁件在航空航天工业中扮演着重要角色，其中，飞机蒙皮是典型的薄壁零件。传统的蒙皮加工常采用化学铣切加工，但由于加工精度低、污染严重等弊端，逐渐被镜像铣加工所替代。镜像铣是一种比较绿色环保的加工方式，加工过程中，双五轴机床在薄壁件的两侧同步移动，一侧为铣削侧，另一侧采用浮动支撑结构对蒙皮进行支撑以增加飞机蒙皮的局部刚度。然而由于蒙皮尺寸大、刚度小等特点，易在切削力的作用下产生弯曲变形，导致轴向切削深度产生偏离，形成加工误差，降低零件的加工精度。因此，有必要在精确预测变形的基础上实现对加工误差的有效补偿，对变形的控制成为飞机蒙皮高加工精度的关键。

为了达到所需的加工精度，在工艺规划过程中必须考虑加工变形，诸多研究人员对薄壁件的变形预测做了大量研究，常用的预测变形方法主要有数值分析法。Kang等^{[ [1]  KANG Y G, WANG Z Q. Two efficient iterative algorithms for error prediction in peripheral milling of thin-walled workpieces considering the in-cutting chip[J]. International Journal of Machine Tools and Manufacture, 2013, 73: 55-61.
1 ]}提出了一种柔性迭代算法和双迭代算法来预测薄壁件的表面形状误差。Shi等^{[ [2]  SHI J H, GAO J, SONG Q H, et al. Dynamic deformation of thin-walled plate with variable thickness under moving milling force[J]. Procedia CIRP, 2017, 58: 311-316.
2 ]}基于第一剪切变形理论和拉格朗日方程讨论了变厚度工件的动态响应，并提出了一种模拟薄壁零件无颤振变形的数值方法。除以上所用方法外，有限元法也是数值分析法中常用的一种。研究人员常使用迭代模型研究薄壁件变形^{[ [3]  SUTHERLAND J W, DEVOR R E. An improved method for cutting force and surface error prediction in flexible end milling systems[J]. Journal of Engineering for Industry, 1986, 108(4): 269-279.
[4]  TSAI J S, LIAO C L. Finite-element modeling of static surface errors in the peripheral milling of thin-walled workpieces[J]. Journal of Materials Processing Technology, 1999, 94(2–3): 235-246.
[5]  WAN M, ZHANG W H, TAN G, et al. Systematic simulation procedure of peripheral milling process of thin-walled workpiece[J]. Journal of Materials Processing Technology, 2008, 197(1–3): 122-131.
3-5 ]}。Ratchev等^{[ [6]  RATCHEV S, LIU S L, HUANG W, et al. Machining simulation and system integration combining FE analysis and cutting mechanics modelling[J]. The International Journal of Advanced Manufacturing Technology, 2007, 35(1): 55-65.
6 ]}提出了一种基于商业ABAQUS的模拟环境，可同时计算切削力和变形。Wang等^{[ [7]  WANG L P, SI H. Machining deformation prediction of thin-walled workpieces in five-axis flank milling[J]. The International Journal of Advanced Manufacturing Technology, 2018, 97(9): 4179-4193.
7 ]}使用商用ANSYS计算了五轴侧铣薄壁零件的变形。吴子腾等^{[ [8]  吴子腾, 张立强, 杨青平, 等. 飞机蒙皮镜像加工变形迭代预测方法[J]. 机床与液压, 2022, 50(20): 25-31.WU Ziteng, ZHANG Liqiang, YANG Qingping, et al. Iterative prediction method for deformation of aircraft skin in mirror machining[J]. Machine Tool & Hydraulics, 2022, 50(20): 25-31.
8 ]}针对蒙皮镜像铣提出了一种加工变形迭代预测方法，然而，由于在迭代时需要反复调用有限元分析，上述方法预测变形的效率很低。为了避免重复进行有限元分析，Li等^{[ [9]  LI Z L, TUYSUZ O, ZHU L M, et al. Surface form error prediction in five-axis flank milling of thin-walled parts[J]. International Journal of Machine Tools and Manufacture, 2018, 128: 21-32.
9 ]}提出了结构刚度修正方法，只需对薄壁零件进行1次有限元分析。Huang等^{[ [10]  HUANG W W, ZHANG Y, ZHANG X Q, et al. Wall thickness error prediction and compensation in end milling of thin-plate parts[J]. Precision Engineering, 2020, 66: 550-563.
10 ]}将结构刚度修正方法和特殊的有限元网格方法相结合，进行端铣过程中的变形预测。Ge等^{[ [11]  GE G Y, XIAO Y K, FENG X B, et al. An efficient prediction method for the dynamic deformation of thin-walled parts in flank milling[J]. Computer-Aided Design, 2022, 152: 103401.
11 ]}通过对工件刚度矩阵进行了简化和动态修改，以避免重新网格化，首次提出了一种有效的变形预测方法，降低了求解过程的计算规模。Li等^{[ [12]  LI W T, WANG L P, YU G. Force-induced deformation prediction and flexible error compensation strategy in flank milling of thin-walled parts[J]. Journal of Materials Processing Technology, 2021, 297: 117258.
[13]  LI W T, WANG L P, YU G. Chatter prediction in flank milling of thin-walled parts considering force-induced deformation[J]. Mechanical Systems and Signal Processing, 2022, 165: 108314.
12-13 ]}提出了一种基于静态子结构方法和结构刚度修正方法的力诱导变形预测模型，通过耦合两个子结构的变形模型，显著降低了整个薄壁件变形模型的自由度，该方法实现了小尺寸薄壁零件的高效、准确加工仿真。Xi等^{[ [14]  XI X L, CAI Y L, WANG H T, et al. A prediction model of the cutting force-induced deformation while considering the removed material impact[J]. The International Journal of Advanced Manufacturing Technology, 2022, 119(3): 1579-1594.
14 ]}通过考虑移除材料对整体刚度矩阵的影响，采用超单元法提高了变形预测的效率。

虽然以上研究所提出的变形预测方法提高了加工效率，但大部分方法是针对传统机床研究的。在蒙皮镜像双五轴加工中，针对蒙皮此类大型薄壁件，使用有限元求解的计算难度大。由于现有研究中的有限元模型大多针对小尺寸薄壁件，且很少有研究考虑加工变形沿刀路的连续性，因此为了减小工件变形的计算负担和揭示实际的变形过程，仍需要一种有效的大型薄壁零件加工变形预测方法。为了解决上述问题，本文提出了一种基于子结构的蒙皮镜像铣削变形预测柔性方法，将薄壁件分为主、子结构两个部分，根据两个部分的刚度方程，建立了整个薄壁件的降阶变形模型，可以有效减少有限元求解的计算时间。通过考虑变形沿刀路的连续性，提出了一种变形预测的柔性模型，将前一个刀具位置的收敛轴向切削深度作为初始迭代值，计算下一个刀具切削位置处零件的变形，所提出的变形预测柔性方法提高了变形预测的效率。

镜像铣系统由两台同步运动的五轴机床组成，一侧为铣削侧，另一侧为支撑侧，如图1所示，支撑侧的6个浮动支撑装置提供支撑力以增强蒙皮的局部刚度，减小蒙皮的变形。在支撑侧配有超声波测厚仪，能通过喷涂射流式耦合剂对蒙皮的壁厚进行实时测量^{[ [15]  吴子腾, 张立强, 杨青平, 等. 蒙皮镜像加工误差实时补偿优化方法研究[J]. 机械科学与技术, 2023,42(4): 1-7.WU Ziteng, ZHANG Liqiang, YANG Qingping, et al. Research on real-time compensation optimization method for skin mirror processing error[J]. Mechanical Science and Technology for Aerospace Engineering, 2023,42(4): 1-7.
15 ]}。蒙皮的四周为柔性夹爪夹持工装，蒙皮加工过程中，当切削力作用时，虽然浮动支撑提供了支撑力，但是在无法支撑的悬空区域内仍会产生加工变形e（x），如图2所示，因此有必要预测该变形的大小以进一步抑制加工变形，减小加工误差并提高加工效率。

在飞机蒙皮镜像铣削过程中，切削力、工件变形和切削深度之间相互耦合，实际轴向切削深度和薄壁件变形，以及初始名义轴向切削深度之间存在如下关系。

式中，y为实际轴向切削深度；H为初始名义轴向切削深度；e（x）为工件变形量。

在加工过程中，由于镜像铣的铣刀刚度远大于蒙皮刚度，故不考虑刀具的轴向变形和径向变形，只考虑工件的变形。

当飞机蒙皮在切削力的作用下产生变形时，实际轴向切削深度减小，如表1第1行所示。切削力与实际轴向切削深度之间为正相关，因此切削力随实际轴向切削深度的减小而减小，由于变形量也与切削力为正相关，因此当切削力减小时，零件的变形量也减小。由式（1）可知，零件的变形量减小时，实际轴向切削深度增大，如表1第2行所示。三者之间的关系如此往复变化，需要不断迭代直至达到平衡状态。切削力、工件变形量和实际轴向切削深度三者之间的关系变化见表1。

由于实际轴向切削深度和蒙皮变形量之间的耦合关系，迭代公式可表示为

式中，y_k为第k次迭代计算的切削深度；e（x_k）为第k次迭代计算的零件变形量。轴向切削深度的迭代更新过程如图3所示，当达到平衡时，收敛条件限制如下。

式中，ε为所允许的加工精度，达到该精度要求时表明变形已经收敛。

飞机蒙皮刚度低，加工时容易产生误差，传统的有限元分析通过求解大型稀疏矩阵方程计算各刀具位置上每个网格节点的变形，即

式中，K为薄壁件有限元模型的刚度矩阵；F、e分别为工件坐标系中各节点沿x、y、z 3个方向的载荷向量和位移向量。

传统有限元方法通过选择适当的边界条件可以获得满意的变形预测效果，但对于飞机蒙皮此类大型薄壁件而言，由于划分的网格单元数量过多，在进行有限元分析时耗时非常长。因此引入子结构技术对薄壁件的有限元模型进行缩减，通过减少变形计算过程中涉及的网格节点数量，避免求解式（4）类的大型方程组。值得注意的是，使用此方法并不会损失结构的复杂性。

将薄壁件划分为两部分，分别为未加工区域和加工区域，如图4所示，薄壁件的四周为夹紧区域，两个部件之间接触区域上的节点和夹紧区域上的节点为边界节点，其他区域为内部节点。将整个未加工区域视为子结构，将加工区域视为主结构（即非子结构），每个结构均包含内部节点和边界节点。

为了建立整个薄壁件的刚度方程，首先建立每个部件的刚度方程，然后利用边界节点的位移兼容性将两个部件的刚度方程进行耦合。主结构和子结构两个部分的刚度方程各分为两种情况（图5），一种情况是边界节点固定，并施加内部载荷，如图5（b）所示，此时边界节点位移为0，该条件下的刚度方程为

式中，下标i和b分别为内部和边界节点；K_jj（j=i，b）为对应的刚度矩阵；e_i，1为结构件内部节点的位移；F_i为等效节点载荷；F_b，1为结构件的反作用力。

另一种情况是结构内部没有载荷作用，实际位移发生在边界节点上，如图5（c）所示，此条件下的刚度方程为

式中，e_i，2为结构件内部节点的位移；F_b，2为边界节点的载荷；e_b为边界节点的位移。

将上述两种情况相加，得到内部节点的实际位移和作用在边界节点上的载荷，可表示为

式中，e_i为结构件内部节点的位移之和；F_b为边界节点的载荷之和。

由式（5）~（8）得到每个结构零件的刚度方程为

将子结构和主结构（即未加工区域和加工区域）两部分的刚度方程进行耦合，由于不同结构之间不存在刚度矩阵元素的交换，因此不同结构之间的刚度矩阵元素为0，进而得到整个薄壁件的变形方程如下。

上标1和2分别为加工区域和未加工区域零件，由于加工时只有加工区域受到铣削力的作用，未加工区域并未受到铣削力作用，所以$F_{\text{i}}^2=0，F_{\text{b}}^2=0$，边界节点刚度和位移之间的关系为

因此，式（10）可简化为

薄壁件变形模型简化后，避免了求解整个薄壁件全局刚度矩阵的逆矩阵，只需要求得内部节点和边界节点的部分刚度矩阵逆解。由于部分刚度矩阵的阶数远小于整体刚度矩阵，因此变形计算效率得到大幅提高。切削力模型使用文献[ [8]  吴子腾, 张立强, 杨青平, 等. 飞机蒙皮镜像加工变形迭代预测方法[J]. 机床与液压, 2022, 50(20): 25-31.WU Ziteng, ZHANG Liqiang, YANG Qingping, et al. Iterative prediction method for deformation of aircraft skin in mirror machining[J]. Machine Tool & Hydraulics, 2022, 50(20): 25-31.
8 ]中所提出的模型，将式（12）代入式（2），进行切削参数的迭代。

传统的变形预测模型在计算除初始刀具位置外某个刀具位置的变形时，都将名义轴向切削深度H作为迭代初始值，即

由图3可知，薄壁件的表面和设计表面分别为迭代空间的上限和下限，传统的预测方法将名义轴向切削深度作为每个刀具位置的初始迭代值，即迭代空间的下限是相同的，没有考虑在相邻刀具位置之间发生的变形之间的关系，两个相邻的刀具位置如图6所示。可知，传统的预测方法忽视了“沿刀具路径的变形也是连续的”这一问题，从而增加了不必要的迭代空间。由于铣削蒙皮时的刀路具有连续轨迹，铣削在进给方向上也是一个连续过程，因此沿刀具轨迹的变形是连续的。根据任意两个相邻刀具位置处的变形连续性原理可知，相比于名义轴向切削深度H，前一个刀具位置m处的收敛轴向切削深度更接近后一个刀具位置m+1处的收敛变形。因此，提出变形预测柔性模型，将前一个刀具位置m处k次迭代后的收敛轴向切削深度y_m，k作为当前刀具位置m+1处的迭代初始值，如图7所示。在刀具位置m+1处，将y_m，k作为迭代空间的下限以预测刀具位置m+1处的变形，相比刀具位置m处，显然缩小了整个迭代空间，加快了变形的收敛速度。所提出的变形预测柔性模型可表示如下。

式中，y_m，k为第m个刀具位置处，k次迭代后收敛的轴向切削深度；y_m+1，k′为第m+1个刀具位置处收敛的轴向切削深度；k′为m+1位置处的迭代次数。与式（13）相比，式（14）的迭代空间更小，表明变形预测柔性模型的收敛速度更快。

使用ABAQUS建立零件的几何模型，蒙皮尺寸为1000 mm×1200 mm×6 mm，直径为1800 mm，材料为7050铝合金，杨氏模量和泊松比分别为70 GPa和0.33。将刀具和浮动支撑装置设置为刚体，浮动支撑装置共提供75 N的支撑力，蒙皮四周设置为完全固定约束，在对模型进行网格划分时，使用四面体网格，长方形槽的网格密度相对更密集，整个模型共划分了125787个单元，255624个节点。对蒙皮进行槽铣试验，机床转速和进给速度分别为8000 r/min和1500 mm/min，轴向切削深度为0.5 mm，长方形槽长度为300 mm，在进行加工变形切削试验时，需要关闭镜像铣机床的误差实时补偿功能。长方形槽上面100 mm为预留下刀长度，下面200 mm为正常切削和仿真位置，并在这段距离间等距测量15个采样点（即刀具位置），镜像铣示意如图8所示，采样点位置如图9所示。

文献[ [8]  吴子腾, 张立强, 杨青平, 等. 飞机蒙皮镜像加工变形迭代预测方法[J]. 机床与液压, 2022, 50(20): 25-31.WU Ziteng, ZHANG Liqiang, YANG Qingping, et al. Iterative prediction method for deformation of aircraft skin in mirror machining[J]. Machine Tool & Hydraulics, 2022, 50(20): 25-31.
8 ]中使用传统的有限元模型和迭代预测方法，将本文所提简化变形模型和柔性方法与文献[ [8]  吴子腾, 张立强, 杨青平, 等. 飞机蒙皮镜像加工变形迭代预测方法[J]. 机床与液压, 2022, 50(20): 25-31.WU Ziteng, ZHANG Liqiang, YANG Qingping, et al. Iterative prediction method for deformation of aircraft skin in mirror machining[J]. Machine Tool & Hydraulics, 2022, 50(20): 25-31.
8 ]中的方法进行比较。首先进行简化变形模型的验证，在验证简化变形模型用于有限元分析时，迭代方式均采用传统的迭代方法，得到15个采样点达到收敛平衡时的仿真数据如图10所示。15个采样点中，以点14为例，其加工变形有限元仿真结果如图11所示，在仿真中，使用传统有限元模型计算完所有切削位置大约需要23400 s，而使用所提简化模型则只需要约12600 s，计算时间减少了46.15%。

随后，对所提变形预测柔性方法进行验证，有限元模型均使用简化变形模型。分别采用传统变形迭代预测方法和柔性方法对长方形槽进行迭代仿真，第一次切削计算时的轴向切削深度为H，铣削力系数使用文献[ [8]  吴子腾, 张立强, 杨青平, 等. 飞机蒙皮镜像加工变形迭代预测方法[J]. 机床与液压, 2022, 50(20): 25-31.WU Ziteng, ZHANG Liqiang, YANG Qingping, et al. Iterative prediction method for deformation of aircraft skin in mirror machining[J]. Machine Tool & Hydraulics, 2022, 50(20): 25-31.
8 ]中的校定数据。在初始刀具位置处使用变形预测传统方法，在其余刀具位置处使用所提柔性方法，蒙皮的加工变形仿真流程如图12所示，在不同刀具位置使用两种方法的迭代次数有限元仿真结果如图13所示，在进行有限元仿真时，在满足式（3）的条件下，加工精度为0.001 mm。从图13可以看出，除初始位置点1外，在其他刀具位置，变形预测柔性方法的迭代次数均小于传统方法。两种方法在所有刀具位置的迭代次数之和分别为36次（变形预测柔性方法）和62次（传统方法），柔性方法的整体迭代次数减少了41.94%，这表明柔性方法具有更快的收敛速度，提高了加工效率。

选取前2个切削位置的加工变形迭代过程，进一步对比两种方法的迭代次数，结果如图14所示。图14（a）是在初始切削位置点1处产生的加工变形迭代过程，在第2个切削位置点2处，分别使用传统变形预测方法和变形预测柔性方法，得到的加工变形迭代过程分别如图14（b）和（c）所示。在使用文献[ [8]  吴子腾, 张立强, 杨青平, 等. 飞机蒙皮镜像加工变形迭代预测方法[J]. 机床与液压, 2022, 50(20): 25-31.WU Ziteng, ZHANG Liqiang, YANG Qingping, et al. Iterative prediction method for deformation of aircraft skin in mirror machining[J]. Machine Tool & Hydraulics, 2022, 50(20): 25-31.
8 ]中的方法进行第1次计算时，轴向切削深度H为0.5 mm，在计算出加工变形量后根据式（13）计算出第1次迭代时的轴向切削深度，即图14（a）和（b）轴向切削深度的初始点数值。在使用变形预测柔性方法进行第1次计算（图14（c））时，轴向切削深度以第1个切削位置处的收敛切削深度为迭代初始值，经过3次迭代后达到收敛平衡。而文献[ [8]  吴子腾, 张立强, 杨青平, 等. 飞机蒙皮镜像加工变形迭代预测方法[J]. 机床与液压, 2022, 50(20): 25-31.WU Ziteng, ZHANG Liqiang, YANG Qingping, et al. Iterative prediction method for deformation of aircraft skin in mirror machining[J]. Machine Tool & Hydraulics, 2022, 50(20): 25-31.
8 ]中的方法由于仍以初始名义轴向切削深度作为迭代初始值，因此经过4次迭代后才达到收敛平衡（图14（b）），进一步说明变形预测柔性方法可以减少迭代次数，提高加工效率。

在5000 mm级飞机蒙皮双五轴镜像铣机床中进行切削试验，采用SINUMERIC 840D SL数控系统，如图15所示。使用激光扫描仪进行蒙皮型面的扫描，将扫描所得点云拟合成实际的工件状态，得到实际工件的数学模型后，通过计算机辅助软件生成刀轨。应注意在加工过程中关闭镜像铣机床的误差实时补偿功能，以保证变形数据的真实性。使用支撑侧末端的超声波测厚仪对长方形槽进行壁厚测量，通过所测壁厚数据计算加工变形量。

取图9所示的15个切削位置进行切削试验，采用变形预测柔性方法预测蒙皮变形，所得变形预测仿真结果和试验结果如图16所示，可以看出，变形预测结果与试验结果一致性较高，在长方形槽中间位置产生的加工变形较大而两端变形较小，这是因为槽的中间处蒙皮刚度较小而两侧靠近夹爪处的刚度较大。此外，仿真试验为理想条件，变形量相对较小，在实际加工中，受夹紧力、双五轴机床同轴度及切削热等因素的影响，蒙皮产生的实际变形大于仿真产生的变形，进一步证实了所提变形预测柔性方法的有效性。

（1）传统的蒙皮加工变形有限元模型计算效率低，在对整个大型飞机蒙皮零件使用有限元方法进行分析时，网格数量多将导致计算耗时较长，本文所提子结构分析方法简化了飞机蒙皮的变形模型，降低了薄壁件刚度方程的阶数，在不损失结构复杂性的前提下提高了计算效率，使整体计算时间减少了46.15%。

（2）所提变形预测柔性方法考虑了沿刀具路径的变形连续性，通过有限元仿真比较所提方法和传统迭代预测方法，结果显示，本文所提柔性方法的迭代次数更少，整体迭代次数减少了41.94%，表明本文所提方法比传统迭代预测方法拥有更快的收敛速度，更高的加工效率。