DOI: 10.12677/mos.2025.141121, PDF, HTML, XML,  被引量   
作者: 孙 宁：上海理工大学机械工程学院，上海
关键词: 航空薄壁件；镜像铣削；充液拉深；去应力退火；多工序仿真；Aerospace Thin-Walled Parts； Mirror Milling； Liquid Filled Deep Drawing； Stress Relief Annealing； Multi Process Simulation

文章引用：孙宁. 航空薄壁件多工序连续工艺链式仿真方法研究[J]. 建模与仿真, 2025, 14(1): 1344-1358. https://doi.org/10.12677/mos.2025.141121

1. 引言

随着航空航天领域对轻量化的需求，钢、铝和钛合金等薄壁构件得到了广泛地应用。由于薄壁零件本身的特性，其切削量大，加工工艺复杂，容易产生变形。铝合金是航空工业中使用最多的薄壁件，具有轻质、轻量化和薄壁化等优点。如何对薄壁零件铣削过程中的变形进行精确控制，是目前航空航天领域的一个难点问题[1]。目前国内外学者普遍认同的观点是：铝合金材料自身的残余应力和铣削后残余应力的再分布是导致铝合金薄壁件铣削变形的根源[2]。目前，国内外已开发出多种成形模拟软件，将工艺设计与成形分析结合起来，并将其用于薄壁零件的设计与制造[3]，通过有限元仿真，可以省去繁琐的试错工作，并在此基础上对零件进行优化设计，从而减少试验和生产成本[4]。

在多工艺条件下，各道次切削都会产生残余应力，这必然会对随后的成形工艺和工件的性能等带来不利的影响[5]，而退火处理则能有效地消除残余应力，提高其随后的成形性[6]。孔祥景[7]利用ABAQUS软件，对航空X形薄壁构件在单道和多道工序铣削时的残余应力进行了研究，提出了一种消除应力退火的方法，并对该方法进行了实验验证。詹梅等[8]以ABAQUS为基础，构建薄壁零件的旋压、回弹和退火的全过程数值模拟模型，并对零件二道次成形时的壁厚分布和工艺参数的影响进行了研究。周建君等[9]在 ABAQUS软件中实现了前、后两个功能模块的二次开发，使复杂的建模程序得以简化，提高了计算的效率。Sugiant等[10]针对复杂铝合金结构开展了淬火过程的数值模拟研究，采用钻芯法对零件内应力进行了检验，证实了数值模拟的正确性。蒋涛[11]采用生死单元方法模拟了柴油机的切削过程，得出了机体在加工过程中的变形规律。

本文拟以火箭发动机箱底薄壁件为对象，开展包括充液拉深、热处理及镜像铣削等多道工艺的数值模拟与试验研究，获得不同工艺参数下的残余应力分布规律，并利用ABAQUS软件对各个工艺参数进行数值模拟，并通过工艺测试来验证模拟结果的正确性，从而为火箭发动机箱底结构的设计、制造和工艺优化提供理论依据。

2. 多工序仿真方法

2.1. 研究对象及工艺选择

本文以直径4000 mm，总厚度6.6 mm的2219铝合金的燃料贮箱底为例展开讨论研究，其各部分特征如图1所示，该贮箱底为半椭球形，表面分布若干个凸台特征，其基体部分厚6.6 mm。2219铝合金的物理力学性能如表1所示。

贮箱底的加工步骤主要分为：充液拉深、去应力退火、镜像铣削内外表面与表面特征加工，各阶段加工参数如表2、表3、表4所示。经过多工序加工后，贮箱底内部往往会产生复杂的残余应力分布，这些残余应力如果未得到合理控制，可能会导致材料的局部失效，进而影响整个燃料贮箱的性能和安全性。因此针对本文的贮箱底案例，将使用ABAQUS软件进行加工仿真，探究多工序加工后贮箱底的残余应力分布情况。

Figure 1. Schematic diagram of fuel storage tank bottom

图1. 燃料贮箱底示意图

Table 1. Physical and mechanical properties of 2219 aluminum alloy

表1. 2219铝合金物理力学性能表

传导率/(W/mK)	质量密度 $\rho$ /(kg∙m−3)	杨氏模量MPa	泊松比	比热[mJ/(ton∙℃)]
159	2.84e−9	71,000	0.33	834,000,000

Table 2. Parameters of liquid filled deep drawing

表2. 充液拉深参数表

拉深位移/(mm)	压边力/(T)	液室压力/(MPa)
950	1400	8.3

Table 3. Stress relief annealing parameters table

表3. 去应力退火参数表

保温温度/(℃)	保温时间/(min)	冷却方式
535	60	随炉冷却

Table 4. Mirror milling parameter table

表4. 镜像铣削参数表

主轴转速/(r/min)	进给量/(mm/min)	切深/(mm)
8400	5000	2

2.2. 充液拉深仿真

利用Abaqus有限元软件对2219铝合金的充液拉深试验进行模拟。有限元模型包括凸模、压边圈和板坯等三个部分，由于凹模已被液体介质完全代替，故不必建模，如图2所示。凸模和压边圈建为刚体模型，板坯建为可变形体，直径为4米，厚度为6.6 mm。随后在Mesh模块中划分网格，单元类型为线性六面体单元C3D8R，采用八节点缩减积分，并在应力集中处进行网格细化，在上下表层按照一层0.1 mm的厚度划分10层，如图3所示，以便后续曲面铣削的有限元仿真，并赋予板坯材料属性，如表5、表6所示。

Figure 2. Key components and assembly relationships

图2. 关键部件及装配关系

Figure 3. Grid division of slab

图3. 板坯网格划分

Table 5. Aluminum alloy 2219 J-C damage

表5. 铝合金2219 J-C损伤

d1	d2	d3	d4	d5	熔点	转变温度	参考应变率	损伤演化
0.13	0.13	1.5	0.011	0	590	350	1E−06	0.04

Table 6. Temperature dependent plasticity of aluminum alloy 2219

表6. 铝合金2219塑性随温度变化

温度(℃)	传导率[mJ/(s∙mm∙℃)]	比热[mJ/(ton∙℃)]	屈服应力(MPa)	塑性应变	杨氏模量(MPa)	泊松比
25	159	834,000,000	107.79	0	71,000	0.33
100	168	838,000,000	103.42		65,193
200	176	880,000,000	100.98		56,262
300	180	964,000,000	71.67		37,980
400	180	1,090,000,000	25.36		31,200
500	180	1,337,000,000	13.1		25,000

在Assembly模块中将各零件进行装配，在Step模块中建立两个动态分析步，如表7所示，设置场输出应力、等效塑性应变和坐标等输出量，并在Interaction模块中设置零件的相互作用属性，建立各接触面之间的相互作用关系，如表8所示。设定凸模在y轴负方向上的位移为950 mm，按照Smooth类型的幅值曲线进行加载，其它方向上的自由度均为零，在板坯下表面施加垂直于表面的均布压力，压力设置为8.3 MPa，模拟液体介质对底面的液压。最后，提交作业。

Table 7. Function of dynamic analysis steps

表7. 动态分析步作用

Step-1	凸模下移、压边圈固定、下表面均布压力
Step-2	凸模和压边圈上移，移除均布压力

Table 8. Interaction contact properties

表8. 相互作用接触属性

接触位置	摩擦类型	摩擦系数
凸模与板坯	罚函数	0.1
压边圈与板坯	罚函数	0.05

对于一些具有旋转轴的圆环类零件，建立柱坐标系可以更好的描述其力学特性，因此在后处理界面建立柱坐标系，其中S11代表径向，S22代表周向，S33代表轴向。仿真结果如图4所示。

Figure 4. Treatment diagram after liquid filled deep drawing

图4. 充液拉深后处理图

2.3. 去应力退火仿真

通过数据传递方法，读取充液拉深后的变形、应力结果到新的求解模型，进行退火仿真，参数如表3所示。

在Abaqus中去应力退火仿真采用顺序耦合的方法，指先进行温度场分析，然后将温度场的结果作为热载荷施加到结构上进行分析，具体流程如图5所示。

在有限元分析中，可以通过设置热传递分析步的加载时间来实现保温时间一个小时，通过设置膜层散热系数与环境温度来实现随炉冷却，如图6所示。其仿真结果如图7所示。

Figure 5. Schematic diagram of sequential coupling process

图5. 顺序耦合流程示意图

Figure 6. Boundary conditions for heat transfer analysis step

图6. 热传递分析步边界条件

Figure 7. Post stress annealing treatment diagram

图7. 去应力退火后处理图

2.4. 镜像铣削仿真

通过数据传递方法，将充液拉深后的变形和应力结果读取到新的求解模型中，以进行镜像铣削仿真，相关参数见表4。在加工过程中，由于曲面复杂，刀具移动导致工件材料不断被去除。为了解决这一问题，本文使用了“动态载荷 + 生死单元”技术，如图8所示。该技术通过离散化连续加工过程，建立若干分析步来模拟实际加工；根据分析顺序，将三维动态铣削仿真得到的铣削力施加到待切除材料单元的节点上，利用“生死单元”功能“杀死”需要去除的材料单元，从而实现铣削力切除工件材料的效果；最后，通过建立任务Job对铣削仿真模型进行求解。

Figure 8. Schematic diagram of dynamic load + birth-death element

图8. 动态载荷 + 生死单元示意图

现有技术通常仅允许在全局坐标系中对零件节点施加载荷，导致施加的铣削力方向无法根据曲面零件的形状进行调整，从而无法有效预测其残余应力变形。本节提出了一种铣削力方向变换的方法，可以有效预测曲面零件的残余应力变形。该方法利用Abaqus软件的二次开发功能，首先提取网格节点坐标，并计算每个节点的变换矩阵。然后，将铣削力进行相应变换并施加到模型中，解决了铣削力方向无法根据零件曲率调整的问题。通过此方法，可以准确施加曲面零件的铣削力场，实现复杂结构件的加工残余应力变形预测。

假设目标节点A的坐标为 $\left(x_1,y_1,z_1\right)$ ，任意相邻节点B的坐标为 $\left(x_2,y_2,z_2\right)$ ，根据式(1)判断出与目标节点距离最近的2个相邻节点。

$d_{AB}=\sqrt{{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2+{\left(z_2-z_1\right)}^2}$ (1)

式中， $d_{AB}$ 是目标点A与任意相邻点B的距离。

接着，建立目标节点A与相邻单元B、C (假设B、C为最近的2个相邻点，C坐标为 $\left(x_3,y_3,z_3\right)$ )的向量，如式(2)和(3)所示。

$AB=\left(x_2-x_1,y_2-y_1,z_2-z_1\right)$ (2)

$AC=\left(x_3-x_1,y_3-y_1,z_3-z_1\right)$ (3)

而目标单元的法向量 $h$ 可以由 $AB$ 与 $AC$ 的叉积得到，如式(4)所示。

$h=AB\times AC=\left[\begin{array}{c}\left(y_2-y_1\right)\cdot \left(z_3-z_1\right)-\left(z_2-z_1\right)\cdot \left(y_3-y_1\right)\\ \left(z_2-z_1\right)\cdot \left(x_3-x_1\right)-\left(x_2-x_1\right)\cdot \left(z_3-z_1\right)\\ \left(x_2-x_1\right)\cdot \left(y_3-y_1\right)-\left(y_2-y_1\right)\cdot \left(x_3-x_1\right)\end{array}\right]=\left[\begin{array}{c}{h}_x\\ h_y\\ h_z\end{array}\right]$ (4)

为了获得变换矩阵，还需对法向量 $h$ 进行标准化，得到等标准法向量 $h_{unit}$ ，如式(5)和(6)所示。

$\Vert h\Vert =\sqrt{h_x^2+h_y^2+h_z^2}$ (5)

$h_{unit}=\left(\frac{h_x}{\Vert h\Vert }，\frac{h_y}{\Vert h\Vert }，\frac{h_z}{\Vert h\Vert }\right)$ (6)

式中， $\Vert h\Vert$ 为法向量 $h$ 的模长。

选择向量 $\alpha =\left(0,0,1\right)$ 作为初始向量，由此计算另外两个正交方向的标准化向量 $u$ 和 $v$ ，如式(7)和(8)所示。

$u_{unit}=\frac{u}{\Vert u\Vert }=\frac{h_{unit}\times \alpha }{\Vert h_{unit}\times \alpha \Vert }$ (7)

$v_{unit}=\frac{v}{\Vert v\Vert }=\frac{h_{unit}\times u_{unit}}{\Vert h_{unit}\times u_{unit}\Vert }$ (8)

由此，可以获得变化矩阵M，如式(9)所示。

$M=\left[\begin{array}{ccc}{u}_{unit}& v_{unit}& h_{unit}\end{array}\right]$ (9)

仿真与实验获得的铣削力数据可以通过式(9)中的变换矩阵进行方向变换并保存于文件中。在Abaqus前处理中，将转换后的铣削力数据输入集中力场，实现曲面铣削力的准确施加。

通过以上步骤，完成Abaqus的二次开发，能够精确对曲面薄壁件施加铣削力场，为研究镜像铣削对变形的影响提供有效方法。由于曲面铣削力难以测量，本节将曲面化为平面，在Abaqus中对小型长方体使用2刃PCD铣刀进行三维铣削仿真，铣削模型和结果如图9和图10所示。测得X方向(刀具进给方向)的铣削力为60 N，Y方向(垂直刀具进给方向)的铣削力为180 N，Z方向(刀具轴线方向)的铣削力为25 N。

将上述铣削力通过式(9)进行方向变换，同时在Abaqus中创建多个分析步以模拟镜像铣削时的刀轨，如图11所示。结合生死单元完成仿真前处理的设置，提交作业后结果如图12所示。

Figure 9. 2 PCD milling cutter model with cutting edge

图9. 2刃PCD铣刀模型

Figure 10. Plane milling model and results

图10. 平面铣削模型及结果

Figure 11. Multi step analysis loading

图11. 多分析步加载

Figure 12. Mirror milling processing diagram

图12. 镜像铣削处理图

3. 实验验证及数据对比

3.1. 实验设备

充液拉深实验采用航天某院研制的4米级整体箱底生产线，通过数控板材流体高压成型机制造出直径4米的火箭整体箱底，液室压力为10 MPa，下压行程900 mm。热处理阶段利用航天某院的热处理自动化生产线对箱底内外表面进行去应力退火处理。镜像铣削阶段则由上海拓璞数控科技有限公司研发的大型飞机曲面蒙皮/火箭箱底双五轴镜像铣削技术与装备完成。最后，采用山东某公司生产的便携式X射线残余应力检测仪对各加工阶段的箱底表面进行残余应力检测。

3.2. 测试数据及结果分析

充液拉深的测试点和现场图如图13所示，其测试结果如表9所示，在Origin软件中绘制相应的柱状图如图14所示。热处理的测试点和现场图如图15所示，其测试结果如表10所示，相应的柱状图如图16所示。镜像铣削的结果如图17、表11以及图18所示。

Figure 13. Distribution diagram of liquid filled deep drawing test points (left) and on-site test diagram (right)

图13. 充液拉深测试点分布图(左)和现场测试图(右)

Table 9. Residual stress values at liquid filled deep drawing test points

表9. 充液拉深测试点残余应力值

测试编号	外壁残余应力值(MPa)
1	−224.6
2	−78.7
3	−75.7
4	52.0
5	27.8
6	15.1
7	8.0
8	−38.4
9	−165.9

Figure 14. Bar chart of residual stress at liquid filled deep drawing test points

图14. 充液拉深测试点残余应力柱状图

在充液拉深时，由于液体在壳体内部施加的压力是均匀的，但由于受几何形状的影响，成型件边缘部分由于受到的约束和压缩作用，导致其边缘部分产生压应力；而中间部分的材料在拉伸过程中会经历较大的变形，导致拉应力的产生。

Figure 15. Distribution diagram of heat treatment test points (left) and on-site test diagram (right)

图15. 热处理测试点分布图(左)和现场测试图(右)

Table 10. Residual stress values at heat treatment test points

表10. 热处理测试点残余应力值

测试编号	外壁残余应力值(MPa)
1	−52.8
2	−53.5
3	−31.2
4	−27.5
5	−37.4
6	−20.9
7	26.1
8	21.5
9	24.0

Figure 16. Bar chart of residual stress at heat treatment test points

图16. 热处理测试点残余应力柱状图

在热处理阶段，材料的晶格结构会产生变化，进而致使内部应力重新分布，残余应力的幅值便会有所减小。此外，过程中还有可能发生相变，相变会使得材料的物理性质以及应力状态随之改变。这种相变现象会造成局部区域出现应力集中或者应力释放的情况，进而形成正应力与负应力。

另外，在热处理过程中，温度分布往往难以达到完全均匀，这就使得材料不同区域的热膨胀与收缩情况不一致。这种温度上的差异会引发应力产生，最终可能造成部分区域呈现拉应力，而其余区域则呈现压应力。

Figure 17. Distribution diagram of test points for mirror milling outer wall (left) and on-site test diagram (right)

图17. 镜像铣削外壁测试点分布图(左)和现场测试图(右)

Table 11. Residual stress values at mirror milling test points

表11. 镜像铣削测试点残余应力值

测试编号	外壁残余应力值(MPa)	测试编号	外壁残余应力值(MPa)
1	−209.6	9	−165.4
2	−104.3	10	−129.1
3	−224.2	11	−119.1
4	−217.3	12	−64.2
5	−136.0	13	−186.8
6	−169.0	14	−198.9
7	−164.8	15	−215.2
8	−230.2

Figure 18. Bar chart of residual stress at mirror milling test points

图18. 镜像铣削测试点残余应力柱状图

在采用2刃铣刀进行铣削加工时，由于切削过程中刀具的切削角度和切削方式，通常会在工件表面产生压应力。铣削过程中，切削产生的热量会导致工件表面温度升高。切削液的使用可以有效降低切削温度，减少热应力的产生。冷却的同时，材料表面由于温度变化而产生的热应力可能会转变为压应力。

3.3. 仿真与实验结果对比

最后，将仿真和实验的结果进行对比，由于实验测的是沿着加工方向，也就是周向的残余力，因此，选择S22作为对照，结果对比如图19、图20以及图21所示。

通过对比，发现二者存在一定的偏差，误差的原因有可能是实验存在误差或者仿真条件无法完全复刻现场工况，但是总的来说，大体趋势是相同的，可以认为本仿真方法具有可靠性。

4. 结论

本文通过实验验证了基于ABAQUS的航空薄壁件多工序连续工艺链式仿真方法的可行性和可靠性，并详细描述了不同工序之间的数据传递技术。采用“动态载荷 + 生死单元”技术、顺序耦合方法和铣削

Figure 19. Comparison of liquid filled deep drawing

图19. 充液拉深对比图

Figure 20. Comparison of stress relief annealing

图20. 去应力退火对比图

Figure 21. Comparison of mirror milling

图21. 镜像铣削对比图

力方向变换方法等关键技术，确保了各工序间的无缝衔接和仿真的准确性。主要得出了多道工序中残余应力分布特点如下。

充液拉深阶段：箱底顶端呈现拉应力，底端呈现压应力。这是由于液体在壳体内施加均匀压力时，受几何形状影响，成型件边缘受约束和压缩产生压应力，中间部分材料拉伸变形产生拉应力。

热处理阶段：经过去应力退火处理后，箱底表面的拉压残余应力分布趋于平衡，且数值相对变小。材料晶格结构变化使内部应力重新分布，同时可能发生的相变以及温度不均匀导致的热膨胀与收缩不一致等因素共同作用，促使残余应力幅值减小。

镜像铣削阶段：采用“动态载荷 + 生死单元”方法进行仿真，结果表明箱底表面大多呈现为压应力。在铣削加工中，刀具切削角度和方式以及切削热等因素影响下，工件表面产生压应力，这对于后续箱底的使用有积极影响，能够减小变形和裂纹产生的可能性，增加其使用寿命。

参考文献

[1]	Ning, H., Zhigang, W., Chengyu, J. and Bing, Z. (2003) Finite Element Method Analysis and Control Stratagem for Machining Deformation of Thin-Walled Components. Journal of Materials Processing Technology, 139, 332-336. https://doi.org/10.1016/s0924-0136(03)00550-8
[2]	金志刚. 大型叶片数控加工精度提升技术分析[J]. 现代制造技术与装备, 2017(12): 131, 134.
[3]	Yao, C., Zhang, J., Cui, M., Tan, L. and Shen, X. (2020) Machining Deformation Prediction of Large Fan Blades Based on Loading Uneven Residual Stress. The International Journal of Advanced Manufacturing Technology, 107, 4345-4356. https://doi.org/10.1007/s00170-020-05316-8
[4]	Yanagimoto, J., Banabic, D., Banu, M. and Madej, L. (2022) Simulation of Metal Forming—Visualization of Invisible Phenomena in the Digital Era. CIRP Annals, 71, 599-622. https://doi.org/10.1016/j.cirp.2022.05.007
[5]	LI, Y., GAN, W., ZHOU, W. and LI, D. (2023) Review on Residual Stress and Its Effects on Manufacturing of Aluminium Alloy Structural Panels with Typical Multi-Processes. Chinese Journal of Aeronautics, 36, 96-124. https://doi.org/10.1016/j.cja.2022.07.020
[6]	Li, J., Carsley, J.E., Stoughton, T.B., Hector, L.G. and Hu, S.J. (2013) Forming Limit Analysis for Two-Stage Forming of 5182-O Aluminum Sheet with Intermediate Annealing. International Journal of Plasticity, 45, 21-43. https://doi.org/10.1016/j.ijplas.2012.10.004
[7]	孔祥景. 多工序航天薄壁件铣削加工残余应力的控制方法研究[D]: [硕士学位论文]. 上海: 上海理工大学, 2021.
[8]	詹梅, 李虎, 杨合, 等. 大型复杂薄壁壳体多道次旋压过程中的壁厚变化[J]. 塑性工程学报, 2008(2): 115-121.
[9]	周建君, 许俊海, 范青山, 等. ABAQUS二次开发在自冲铆接模拟中的研究[J]. 机械制造与自动化, 2021, 50(5): 146-148.
[10]	Sugianto, A., Narazaki, M., Kogawara, M. and Shirayori, A. (2009) A Comparative Study on Determination Method of Heat Transfer Coefficient Using Inverse Heat Transfer and Iterative Modification. Journal of Materials Processing Technology, 209, 4627-4632. https://doi.org/10.1016/j.jmatprotec.2008.10.016
[11]	蒋涛. 柴油机机体残余应力模拟与加工变形研究[D]: [硕士学位论文]. 镇江: 江苏科技大学, 2017.