飞机整体结构件加工变形的初始残余应力- 初始几何误差耦合影响与控制

2021-03-05林锋姚婉秦国华叶海潮陶江

兵工学报 2021年12期

林锋，姚婉，秦国华，叶海潮，陶江

(南昌航空大学航空制造工程学院，江西南昌 330063)

0 引言

航空航天领域的发展水平在当前已成为衡量一个国家综合国力的重要指标之一。人们对于航空材料性能的要求越来越高，而航空材料是一类特殊且要求比较严苛的材料。铝合金由于较高的比强度、低密度、良好的加工性能以及耐腐蚀等优点，广泛应用于航空领域[1]。此外，为了提高飞行机动性，许多骨架零件大量采用整体结构件[2]。在加工过程中，毛坯材料的去除率高达90%，使得零件的刚性变差，易发生加工变形且加工精度难以得到保证，每年因为加工变形而导致报废的零件数不胜数，造成巨大的经济损失，也限制了我国航空业的进一步发展。

零件在高速加工时变形的主要致因包括初始残余应力和加工残余应力[3-4]。Li等[5]介绍两种残余应力的产生原因及影响因素，并分析了这两种残余应力的预测、控制方法及其适用范围。据统计，航空零件加工过程中发生变形的约90%可以归结于毛坯内部的初始残余应力[6]。Sun等[7]根据板变形理论进行了不同初始残余应力作用下的铣削变形实验，发现初始压应力和加工残余应力耦合下板的变形会增加，而初始拉应力和加工残余应力作用下板的变形会减小。Huang等[8]分析毛坯初始残余应力和铣削应力耦合作用对零件变形的影响，对比化学铣削实验结果后认为，在高速铣削情况下，残余应力对零件变形的影响远远大于铣削应力。Wang等[9]和Lu等[10]基于航空整体薄壁件加工变形提出一种仅在残余应力作用下的加工变形预测模型，该模型可以较准确地预测出残余应力对加工变形的影响以及残余应力分布中行测的偏移对变形的影响。Li等[11]利用零件的简化模型进行有限元建模，分析了双向残余应力对变形的影响，发现初始轴向残余应力对变形的影响大于径向初始残余应力对变形的影响。Cerutti等[12]研究初始残余应力的作用下加工顺序对加工质量的影响，提出一种计算工具，利用该工具根据初始工件和加工顺序来优化加工质量。Li等[13]依据两个典型航空零件在加工后72 h内的变形规律提出三步加工工艺，以控制薄壁零件的加工变形，实验结果表明在第1次加工后将零件自由放置72 h可以有效降低约48%变形值。秦国华等[14]依据弯曲变形理论推导出零件加工变形的力学模型，将加工要求作为加工变形的约束条件，结合夹逼准则以及引入单位变量，将加工变形转化为齐次线性不等式，通过改变约束条件来求解方程组，以控制加工变形。

国内外学者对于整体结构件加工变形方面的研究较为深入，但目前对于毛坯的初始几何误差这一方面的研究较少。在实际生产加工过程中，由于材料的非均匀性，毛坯表面不可避免地存在初始几何误差。

为了进一步揭示毛坯不同数值大小初始几何误差对零件变形的影响，本文针对不同的毛坯初始几何误差建立相应的有限元模型，并选取零件在毛坯内的不同位置进行研究，得到不同初始几何误差以及零件位于毛坯内不同位置对零件加工变形的影响规律。再以零件加工后的变形最小为目标，利用变向迭代优化方法搜索零件在毛坯内的最佳位置，即当零件处于毛坯内的该位置时，零件的变形为最小，进而可以通过优化零件的加工位置来控制变形。

1 零件变形的分析方法与验证

在毛坯制造过程中，由于材料力学性能的非均匀性，使得铝合金厚板内不仅存在初始残余应力，而且存在初始几何误差(即预拉伸板的平面度)。而在高速铣削毛坯成型为飞机整体结构件的工艺过程中，毛坯内部的初始残余应力得以释放，破坏了飞机整体结构件的静力平衡状态，造成飞机整体结构件发生变形。

1.1 变形理论

由于铝合金厚板沿厚度方向的残余应力值非常小，完全可以忽略不计[15]。因此在自然状态下，铝合金厚板内部的初始残余应力在任意截面都处于自平衡状态，即内部的力和力矩处于自平衡(见图1)，可以表示如下：

图1 毛坯初始残余应力分布Fig.1 Initial residual stress distribution of blank

(1)

(2)

式中：X、Y、Z分别为轧制方向、横向方向、厚度方向；x、y、z分别为沿X轴方向、Y轴方向、Z轴方向的位置；FX(x)、FY(y)分别为X轴方向、Y轴方向的力；MX(y)、MY(x)分别为绕X轴方向、Y轴方向的力矩；σX(z)、σY(z)分别为沿X轴方向、Y轴方向的残余应力；ΩOYZ、ΩOZX分别为x、y处的截面面积。

梁类结构件变形的特点为：轧制方向Y轴上弯矩产生的弯曲变形远大于其他方向上的变形，且扭转变形极小。根据弹性变形理论，可知轧制方向上的应变为

(3)

式中：ν、E分别为零件材料的泊松比、杨氏模量。

因此，为了更加清晰地分析轧制方向上的变形，定义零件轧制方向上的等效应力如下：

σ(z)=σX(z)-νσY(z).

(4)

在残余应力分布σ(z)的作用下，梁类结构件发生了弯曲变形(即挠度)。由于一般整体结构件不同位置的截面是不相同的，在Y轴的不同位置，梁的每个弯矩也是不同的。因此不同截面绕X轴方向的力矩可以表示为

(5)

式中：Ωw(y)表示零件在任意位置y处的截面面积；z0(y)表示零件任意位置y处的截面中性轴高度，如图2所示。图2中，h为零件在毛坯中的加工位置，X′、Y′、Z′为经过零件中性层的3个坐标轴，UY(yA)、UZ(yA)、UZ(yB)分别为A点Y轴、A点Z轴B点Z轴方向的位移，UZ(y)为任意位置y处Z′轴方向的位移。

图2 受力状态Fig.2 Stress state of parts

根据材料力学可知，所选取截面处弯曲变形向下凹时，截面上的弯矩为正、反之为负。因此，对于梁结构，任意位置y处上弯矩引起的挠度UZ(y)可以表示为

(6)

式中：IX(y)为任意位置y对应的截面惯性矩。

1.2 三框件变形计算

图3所示为飞机整体三框结构件，其材料选用7075-T7451铝合金板材，弹性模量E=71 700 MPa，泊松比ν=0.33.图3中，l、w、hf为三框结构件的长、宽、高，l1、l2、l3分别为3个框的内长，hb为腹板的厚度，tf为框的厚度。

图3 三框结构件示意图Fig.3 Schematic diagram of three-frame structure

选用的毛坯外廓尺寸为1 200 mm×120 mm×60 mm，毛坯初始几何误差Δ=0.6 mm. 在采用裂纹柔度法对残余应力进行测量时，每隔1.5 mm测量一次，共测出40组应力值，最后得到毛坯内初始残余应力分布曲线，如图4所示。另外，零件在毛坯中的厚向加工位置为正向h=16.5 mm，如图2(a)所示。

图4 残余应力分布曲线图Fig.4 Distribution curves of residual stress

如图5(a)所示，将毛坯沿厚度方向划分为40层，由于零件的加工位置为16.5 mm，故处于毛坯的第12层～第31层，共计20层。

图5 有限元仿真结果Fig.5 Finite element simulated results

在有限元仿真时，按照毛坯和零件尺寸进行建模，建模过程如下：

1)在无初始几何误差的毛坯中分割出零件和去除材料两个部分。

2)将毛坯均匀划分为40层。

3)选取单元类型C3D20R进行网格划分，设置网格密度15 mm.

4)给毛坯底部施加均布载荷，使其产生一定量的初始几何误差。

5)重新导入具有初始几何误差的毛坯(不导入由均布载荷带来的附加内部应力)，利用网格编辑功能调整局部节点位置，确保零件在毛坯内处于指定位置。

6)给毛坯施加初始残余应力。

采用“单元生死”技术模拟铣削过程的材料去除效应。将实验测量得到的初始残余应力值赋到毛坯中，得到毛坯存在0.6 mm(方向为上凸)初始几何误差，有限元仿真中原始模型和简化模型的结果如图5所示。

对比图5(b)和图5(c)可以看出，简化结构与原始结构的有限元仿真结果非常相近，最大的相对误差仅有3.5%. 因此在解析法中可以利用简化结构来进行计算。

假定毛坯具有初始几何误差Δ，则在极坐标系统中，根据图4的应力值，结合(4)式计算出等效残余应力的等高曲线为

σ(ρ)=5.445 3+2.498 9ρ-2.211 7ρ2+
0.374 2ρ3-0.027 2ρ4+0.001 02ρ5-
0.000 020 7ρ6+0.000 000 215ρ7-
9×10-10ρ8，

(7)

式中：ρ为位置z处由Δ形成的圆弧半径；σ(ρ)为半径ρ处的应力值。

图6所示为应力等高线。图6中，R为毛坯底面处由Δ形成的圆弧半径，θ为半径ρ与Y″轴的夹角。根据图6的几何关系，由于应力等高曲线为以O″为圆心的圆弧，则有

图6 应力等高线Fig.6 Stress contour

(8)

(9)

(10)

尺寸l=1 100 mm代入(10)式后，再代入(8)式，经整理可得

(11)

简化后零件在任意y处都具有相同的截面，因此三框结构件的中性轴高度位置和截面惯性矩分别为

(12)

(13)

三框结构件在不同位置y处，横截面不同。因此可以计算任意位置y处截面绕X轴方向的力矩为

(14)

(14)式代入(6)式，二次积分即可得到挠度曲线。但由于在进行积分计算时无法求解出解析解，沿Y轴方向以50 mm为间距，依次从y=0 mm到y=1 100 mm进行力矩的计算，求解得到一系列力矩值，如表1所示。

由(6)式和(14)式可知，MX(y)是一个9次多项式。根据表1数据，采用多项式进行拟合即可求出弯矩MX(y)表达式，即

MX(y)=-3.996×10-35y9+2.31×10-21y8-1.016×10-17y7+1.79×10-14y6-1.6×10-11y5+7.665×10-9y4-1.91×10-8y3-1.0×10-7y2+0.013 6y-4 284.9. (15)表1 不同位置处的力矩Tab.1 Torques at different locations

由于高次项系数太小，不予考虑，进一步可得力矩表达式为

MX(y)=-1.0×10-7y2+0.013 6y-4 284.9.

(16)

(16)式代入(6)式，可以求得厚度方向上挠曲线方程为

UZ(y)=-1.66×10-6y2+αy+β，

(17)

式中：α、β为常数。

根据图2(b)中零件首尾两端的初始变形值条件，可求得零件在轧制方向上的挠曲线方程为

UZ(y)=-1.66×10-6y2+1.826×10-3y.

(18)

同样，原始结构的每一段挠度曲线也可用同样方法求得：

(19)

1.3 试验验证与分析

为了验证解析法和有限元仿真方法的正确性，下面进行试验验证。试验在意大利菲迪亚公司生产的高速铣床K211A上进行零件加工，在青岛前哨精密机械有限责任公司生产的桥式三坐标测量机ADVANTAGE 15.30.10进行变形测量。整个试验过程如下：

步骤1对毛坯进行无应力装夹。

无应力装夹是指该装夹不对毛坯产生附加的应力，一般采用压板沿四周顶住毛坯侧面，如图7所示。图7中，Q、P分别为毛坯的正面和反面。

图7 毛坯无应力装夹Fig.7 Blank clamping without stress

步骤2对毛坯P、Q表面反复粗铣。

高速铣削毛坯Q面和P面，获得零件在毛坯中的加工位置，具体工序如表2所示。

表2 毛坯表面的粗铣工艺路线Tab.2 Rough milling process of blank surface

步骤3对粗铣后的零件进行有应力装夹。

有应力装夹是指对零件施加夹紧力，将阻止材料去除过程中应力释放造成的零件变形。故采用压板在零件顶面压紧的装夹方式。

步骤4对零件进行精加工。

在意大利菲迪亚公司生产的K211A 3500×1500型铣床上加工成形为三框结构件，如图8所示。加工3个框的内腔时，采用直径为16 mm、刀头数为4、刃长为55 mm的硬质合金铣刀；加工零件外侧缘条时，则采用直径为12 mm、刀头数为4的硬质合金铣刀。主轴转速均为15 000 r/min，进给量均为0.083 mm/r. 为了尽量避免铣削对加工变形的影响，在铣削外侧缘条时，轴向切深仅为0.5 mm.

图8 加工后的三框结构件Fig.8 Three-frame structure after machining

步骤5对三框结构件进行变形测量。

零件变形的测量在青岛前哨精密机械有限责任公司生产的桥式三坐标测量机ADVANTAGE 15.30.10上进行。测量位置为零件腹板底面中线的数据，每隔50 mm进行一组数据的测量。

图9所示为零件在毛坯内的加工位置分别为正向(零件开口朝向Q面)16.5 mm和反向(零件开口朝向P面)16.5 mm时的试验值、计算值和仿真值的比较。

由图9可见：加工位置为正向16.5 mm时，最大加工变形的试验测量值为0.585 mm，毛坯存在0.6 mm上凸初始几何误差、无初始几何误差的计算值分别为0.527 mm、0.51 mm，其最大相对误差分别为9.54%和12.8%，对应的有限元仿真值分别为0.535 mm、0.525 mm，最大相对误差分别为9.4%、10.42%；当零件的加工位置为反向16.5 mm时，最大加工变形的试验测量值为-0.323 2 mm，毛坯存在0.6 mm上凸初始几何误差、毛坯无初始几何误差的零件变形计算值分别为-0.320 6 mm、-0.301 8 mm，其相对误差分别为14.5%和17.988%，相应的有限元仿真值分别为-0.325 9 mm、-0.301 mm，相对误差分别为11.95%和17.6%. 由此可知，不论是有限元仿真还是计算值，毛坯存在初始几何误差时加工变形数据会更为接近试验值。

图9 零件变形试验值、仿真值与计算值Fig.9 Experimental simulated and calculated values of part deformation

误差存在的主要原因有以下两个方面：其一，计算值中所使用的残余应力值是人为试验测量得到的，设备条件和试验方法以及测量精度的限制不可避免地带来误差，因此计算弯矩时，所使用的残余应力值是通过多项式拟合得到的曲线，会带来一定的偏差；其二，对于铝合金厚板而言残余应力在同一厚度方向上的数值并不相同，在有限元分析过程中假定每一层的残余应力值是中心平面的数值，使得仿真结果存在误差。

2 零件变形演变机制

由第1节可知，零件变形的影响因素主要有两个方面：一是毛坯的初始残余应力分布；二是零件的结构形状。

显然，如果毛坯具有初始几何误差，则必将对零件变形产生影响。其本质就是毛坯初始几何误差导致零件受到的残余应力发生了变化，如图10所示。

图10 初始几何误差对受力状态的影响Fig.10 Influence of initial geometric error on stress state

这里，分析毛坯几何误差Δ取值为0～4 mm(方向为下凹)对零件变形的影响规律。事实上，零件在毛坯中的位置不同，也会影响零件的变形，故选取4个不同的位置进行分析，如图11所示，分别为正向16.5 mm、正向9.0 mm、反向16.5 mm和反向9.0 mm.

图11 零件在毛坯内位置示意图Fig.11 Schematic diagram of the positions of parts in the blank

由于受到应力的腹板是零件变形的关键部位，表3列出了腹板所受到的应力。

利用第1节建立的有限元方法，可以获得各个位置处毛坯初始几何误差对零件变形的影响，如图12所示。

由图12(a)和表3可知，当零件正向位于毛坯内16.5 mm时，腹板受到的等效应力随着毛坯初始几何误差的增加呈现逐渐减小的趋势，而且均为拉应力。因此，零件变形向上凸，且随着初始几何误差的增加而减小。

由图12(b)和表3可知，当零件反向位于毛坯内16.5 mm时，腹板受到的等效应力值随着毛坯初始几何误差的增加呈现逐渐减小的趋势，且为压应力。因此，零件变形向上凸，且随着初始几何误差的增加而减小。

表3 腹板受到的应力值Tab.3 Stress values of web plate

由图12(c)和表3可知：当零件正向位于毛坯内9 mm，在初始几何误差为0～1.0 mm时，腹板受到的等效应力为压应力，且呈现逐渐减小的趋势，正向的零件变形为下凹，且变形逐渐减小；在初始几何误差为1.0～4.0 mm时，腹板受到的等效应力值为拉应力，且逐渐增加，则正向的零件变形也逐渐增加，变形向上凸。

由图12(d)和表3可知：当零件反向位于毛坯内9 mm，腹板受到的等效应力均为压应力，初始几何误差为0～1.0 mm时，其值呈增加趋势，故反向的零件变形上凸，且随着初始几何误差的增加而增加；初始几何误差为1.0～4.0 mm时，其值呈减小趋势，故反向的零件变形上凸，且随着初始几何误差的增加而减小。

图12 不同初始几何误差下零件变形规律Fig.12 Deformation law of parts with different initial geometrir errors

由此可知，毛坯初始几何误差对零件变形的影响规律主要有以下3个方面：

1)零件在毛坯中的加工位置不同时，无论是否存在初始几何误差，零件变形也不相同；

2)零件在毛坯中的加工位置相同时，在相同的毛坯初始几何误差情况下，若零件结构出现变化，则零件变形曲线也会出现波动或拐点；

3)毛坯初始几何误差导致零件腹板受到的应力，若为拉应力，则变形向下凹；若为压应力，则变形向上凸。若受到的应力随毛坯初始几何误差的增大而增大，则变形也增大；若受到的应力随毛坯初始几何误差的增大而减小，则变形也减小。

3 加工位置的优化方法

由变形规律可知，在毛坯初始几何误差一定时，零件在毛坯内的不同位置具有不同的变形量。因此，可通过优化加工位置来实现零件变形的控制[16-17]。

3.1 优化模型

在高速切削加工过程中，由于加工过程中产生的附加切削应力对零件变形的影响很小，优化切削工艺参数意义不大。因此，可以通过搜索零件在毛坯内的最佳加工位置，实现零件在加工过程中的变形控制。加工位置的优化模型建立如下：

findz，
minUX(y)，
s.t.zmin≤z≤zmax，

(20)

式中；zmin和zmax分别为零件在毛坯内的最低厚向位置和最高厚向位置。

3.2 步长递减算法

为了求解出(16)式中z的解，可按照下列步长递减算法来进行。

步骤1确定加工位置的初始近似值。

(21)

步骤2确定加工位置的第1个近似值。

(22)

根据加工位置的第1个近似值与初始近似值对应的变形方向，定义变形方向变化标识函数为

k1=λ1-λ0.

(23)

若k1=0，则表示变形方向未发生变化；若k1≠0，则表示变形方向发生了改变。

步骤3确定加工位置的下一个近似值。

根据相邻两个近似位置处变形标识的变化，确定当前近似位置的步长sn为

sn=(-1)kn-1η|kn-1|sn-1，

(24)

式中：η为0～1之间的任意数，也称递减系数，η∈[0，1]；kn-1=λn-1-λn-2，k0=0，n为加工位置近似值的次数，n≥1；λn-1、λn-2分别为对应于第n-1个近似值zn-1、第n-2个近似值zn-2的变形方向。

这样，加工位置的下一个近似值可表示为

zn=zn-1+sn-1.

(25)

步骤4判断当前近似值的范围。

根据(16)式可知，加工位置的约束条件为[zmin，zmax]，则当前近似值zn应该满足下列关系：

zmin≤zn≤zmax.

(26)

若znzmax，则zn=zmax+(-1)ηsn-1.

步骤5判断终止条件。

由(20)式可知，当前步长可以进一步描述为

(27)

则当且仅当

|sn|≤ε

(28)

时迭代过程终止。

此时，最佳加工位置zo即为加工位置第n次的近似值，即

zo=zn.

(29)

4 应用与分析

以典型飞机整体三框结构件的高速铣削为例，说明其在无初始几何误差的毛坯和有初始几何误差的毛坯中加工变形的控制策略。

零件结构尺寸如图3所示，所使用的毛坯结构尺寸和初始残余应力曲线如图4所示。毛坯和零件的厚度分别为60 mm和30 mm.若毛坯无初始几何误差，则零件在毛坯内的变动范围是zmin=0 mm，zmax=30 mm.给定递减系数η=0.5，阈值ε=0.000 1，初始步长为s=30，根据步长递减算法，最佳加工位置的搜索过程如表4所示。由表4可知：当毛坯不存在初始几何误差时，步长迭代算法在n=27时，当前步长为sn=0.000 1.由于|sn|<ε，优化过程达到收敛。从理论上讲，零件在毛坯内的最佳位置z=23.265 0 mm.因此，从工程实际考虑，最终选取最佳位置z=23.265 0 mm.表5所示为有初始几何误差时零件最佳位置搜索过程。

表4 无初始几何误差时零件最佳位置搜索过程Tab.4 Searching process of optimal positions of partswithout initial geometric error

由表5可知：毛坯存在初始几何误差Δ=0.6 mm时，零件在毛坯内的变动范围是zmin=0.6 mm，zmax=29.4 mm. 同样地，给定递减系数η=0.5，阈值ε=0.000 1，初始步长s=30；步长迭代算法在n=30时，当前步长sn=-0.000 1，由于|sn|<ε，优化过程达到收敛。因此，由算法可知零件在毛坯内的最佳位置z=22.754 3 mm，故实际应用中选取最佳位置z=22.754 3 mm.

表5 有初始几何误差时零件最佳位置搜索过程Tab.5 Searching process of optimal positions ofparts with initial geometric error

通过有限元仿真得到毛坯存在0.6 mm下的凹初始几何误差，零件位于毛坯内正向23.265 0 mm时零件的变形为-0.144 mm；当毛坯不存在初始几何误差，零件位于毛坯内正向22.754 3 mm时，零件的最终变形为0.043 mm. 这两组数据进一步表明毛坯的初始几何误差对零件加工变形具有重要影响。