壁板铆接变形分析及铆接顺序规划研究*

王仲奇，张津铭，常正平，景泽旺，王志敏，罗 群

（1. 西北工业大学机电学院，西安 710072；2. 航空工业西安飞机工业(集团)有限责任公司，西安 710089）

壁板结构因具有质量轻、承载效率高等优点而广泛地应用于飞机结构制造中。壁板也是保证飞机气动外形的重要组件，提高其外形装配精度对保证飞行性能具有重要作用。壁板装配主要采用铆接或螺接，由于其组成零件（蒙皮和长桁）壁薄、尺寸大、刚度小，导致壁板在铆接过程中易产生装配变形，严重影响飞机的外形装配精度和疲劳寿命。为提高壁板的铆接效率、铆接质量及稳定性，以自动钻铆系统为代表的自动化连接设备因铆接工艺参数可控，逐步取代了传统的人工铆接[1]。壁板压铆时，铆钉受到沿钉杆轴向的压铆力，使钉孔外钉杆部分发生镦粗变形，形成镦头；钉孔内钉杆部分产生不均匀膨胀变形，并挤压孔壁形成不均匀分布的干涉量。壁板在镦头挤压和钉孔不均匀膨胀的作用下，不可避免地产生局部变形。壁板装配一般涉及成百上千个铆钉的安装，大量单钉铆接产生的局部变形，随自动钻铆过程不断进行传递、耦合，影响了壁板的外形装配精度。保障飞机外形装配精度也是新一代飞行器研制的重要内容[2]，因此，开展壁板铆接变形分析对提高壁板装配精度具有一定指导意义。

近年来，国内外学者在壁板铆接方面开展了大量研究。M ü ller[3]针对铝合金及复合材料结构的连接，研究了压铆力与镦头尺寸之间的关系及压铆力对钉孔干涉状态的影响规律，并结合疲劳试验分析了压铆力对结构服役寿命的影响。Rijck等[4]研究了压铆力和镦头尺寸之间的关系，认为镦头直径和高度可用于评估压铆力对铆接质量和接头疲劳特性的影响。Cheraghi等[5]利用有限元模型，研究了压铆力、铆钉长度、铆钉直径公差、连接孔直径公差和锪窝深度等工艺参数对连接质量一致性的影响规律。Cunha等[6]提出了一种获取航空航天结构件铆接残余应变的方法，首先利用有限元方法建立单钉铆接变形模型，采集单钉铆接后钉孔孔边产生的位移，并结合线性模型推断一排钉连续铆接后的变形，研究表明预测结果与试验结果接近，可将该方法推广到更大壁板的装配过程中。Zhang等[7]将铆接过程变形划分为弹性、塑性和回弹3个阶段，并建立了相应的力学模型。牟伟强等[8]针对无头铆钉压铆过程，结合主应力法和幂指数硬化理论，构建了压铆力与镦头尺寸之间的力学模型，并通过模拟结果进行验证。刘平[9]指出压铆过程中镦头的成形遵循圆柱体镦粗过程的金属流动规律，并将压铆过程分为5个阶段，计算了压铆过程中需要的压铆力，最后运用有限元实例进行验证。冯晓旻等[10]模拟单个铆钉压铆的动态过程，研究了铆钉和母材在压铆过程中的受力和变形情况。

在铆接变形控制方面，Aman等[11]研究了单排3钉铆接结构中铆接顺序、铆钉间距和被连接件间间隙对铆接质量的影响。结果表明，随着铆钉间距增大，铆钉顺序对被连接件的膨胀变形和残余应力的影响可忽略不计，但铆钉间距和被连接件间间隙对铆接变形影响很大；铆接顺序对铆钉中残余应力和铆接后的间隙影响较小，较小的被连接件间间隙和较大的铆钉间距能够减小铆接变形。壁板铆接变形是大量铆接局部変形在时间和空间分布作用下的综合效果，对增量铆接过程进行数值模拟计算量巨大，通过合理的等效达到计算效率的提高且对精度影响有限是十分必要的[12]。Wang[13]将单钉铆接变形简化为钉孔受内压作用引起，并基于主应力法和厚壁筒受压理论构建了铆接力学等效单元，结合有限元模拟得出了水平安定面铆接变形最小的铆接顺序。在薄壁焊接结构中多采用“局部-整体”模型进行焊接变形计算的简化，如李娅娜等[14]利用“局部-整体”映射加载方法，通过仿真T形结构角焊缝局部模型，得到了局部塑性应变，并将应变映射到薄壁箱形梁的整体模型上，实现了薄壁箱形梁的焊接变形分析，仿真表明焊接顺序对焊接变形影响较大。李江飞[15]研究了同向、反向、交替等焊接顺序对薄板变形的影响，得到最优焊接顺序为反向焊接。周广涛等[16]通过建立数值计算模型，以顶板焊接变形最小为目标，对焊接顺序进行优化，并以试验验证了得到的最优焊接顺序。Masters等[17]基于“局部-整体”映射模拟方法，研究了单个自冲铆钉安装过程形成的变形对整体装配变形的影响，并将平板和顶帽装配的模拟结果与试验结果进行比较，表明该方法在早期的设计阶段可以预测连接变形。

综上所述，目前对于壁板铆接问题的研究，主要集中于铆接件的铆接质量、疲劳寿命、压铆力建模以及少钉结构局部变形分析，而对壁板整体铆接变形较少涉及。铆接顺序是影响铆接变形的主要因素之一，本文借助有限元，以单钉铆接变形分析为元模型，通过壁板铆接变形“局部-整体”映射加载，实现了壁板增量铆接变形计算，并结合智能算法，以铆接变形最小为目标，对铆接顺序进行优化。

壁板铆接变形影响因素分析

在铆接过程中，影响壁板铆接变形的因素有很多，如图1所示。壁板以不同铆接顺序进行铆接时，其变形差异很大，因此，本文主要针对铆接顺序这一变形影响因素展开研究。

图1 壁板铆接变形影响因素分析Fig.1 Analysis of influencing factors on panel riveting deformation

壁板铆接变形数值计算模型构建

针对铆接变形的仿真，通常是构建壁板、铆钉等零件的三维实体模型，设置材料属性、载荷及边界条件，然后进行全过程弹塑性模拟。由于所构建的是实体模型，单元结点数量大，计算耗时，所研究的对象也多限于单钉和几个铆钉的铆接过程。但在壁板实际装配中，往往涉及成百上千铆钉的安装，仅仅通过少量铆钉的变形仿真，难以预测出壁板铆接的整体变形。

壁板在铆接时，钉孔周围处于复杂的3向应力状态，而离钉孔轴线2.5R（R为铆钉半径）以外的部分，壁板大致处于平面应力状态[18]。采用有限元仿真时，钉孔周围受不均匀挤压，应力状态复杂，采用实体单元进行模拟。壁板其余区域出现的弯曲现象，采用壳单元可很好地进行模拟。因此，为了克服传统分析手段的缺点，本文采用“体-壳”连接方法，对三维实体模型进行简化；通过“局部-整体”映射加载，将钉孔周围复杂的应力应变状态传递到壁板整体薄壳模型上，实现壁板增量铆接过程的仿真。

1 壁板铆接体-壳连接建模

体-壳连接如图2所示，图2（a）是16节点等参实体元a，图2（b）是8节点超参壳元b。a单元的ξ=1面与b单元的ξ=-1面相互联结。定义a单元ξ=1面上的节点用1t、1b、2t、2b、3t、3b表示，它们的节点位移参数是 ui、vi、wi（i=1t，1b，2t，2b，3t，3b）；b 单元 ξ=-1 面上的节点用1、2、3表示，它们的节点位移参数为ui、vi、wi、αi、βi（i=1，2，3）。两单元中 ui、vi、wi是节点在总体坐标系内的位移分量，αi、βi是联结节点 ib、it的向量v3i绕与之相互垂直的两个向量v1i和v2i的转动。

为了得到两个单元结点位移参数的约束方程，首先将各结点位移参数转化到局部坐标系 v1i、v2i、v3i中。转化关系为：

式中，u'i、v'i、w'i是局部坐标系 v1i、v2i、v3i的3个坐标分量。

为使得两种单元在变形后联结面依旧贴合并相互协调，局部坐标系内结点位移参数间的约束方程为：

式中，ti代表壳元中节点i所对应的厚度。

当实体单元为20结点等参单元，和壳单元结点i对应的实体单元上有 3 个结点 it、im、ib，im为 it、ib的中点，同理可得两种单元的位移参数约束方程为：

2 壁板铆接“局部-整体”映射加载

壁板铆接变形是由每个铆接过程所产生的局部塑性变形累积引起的。壁板在铆接过程中，塑性变形主要集中在铆钉孔周围很小的区域内，而钉孔区域以外的区域仅发生弹性变形。采用“局部-整体”映射加载实现铆接变形预测步骤如下：

（1）从壁板整体结构中选择单钉铆接局部结构，建立单钉铆接元模型并进行仿真，提取钉孔周围节点的坐标信息和位移信息，获得局部变形场。

（2）将单钉铆接元模型的局部变形场结果映射到壁板铆接整体模型的对应节点上。

（3）通过有限元仿真重启动技术，完成整体模型铆接变形计算，流程如图3所示。

图2 体-壳连接示意图Fig.2 Schematic of solid-shell combination model

图3 壁板铆接变形数值计算流程Fig.3 Numerical calculation process of panel riveting deformation

3 壁板铆接变形数值计算模型

3.1 单钉铆接元模型

国内外学者借助有限元方法对单钉铆接变形进行分析，取得了较为丰硕的成果。如文献[1]中利用ABAQUS软件，建立了平锥头铆钉有限元模型，并采用G86型钻铆机进行铆接试验验证，结果表明所建立的有限元模型具有有效性。本文单钉铆接元模型为三维动态显式模型，包括壁板，上、下铆模，上、下压力脚衬套以及铆钉，如图4所示。单钉铆接元模型的建模过程与文献[1]相同，其计算结果验证不再赘述。通过对单钉铆接元模型进行压铆过程仿真，得到局部变形场。图5所示为压铆仿真后单钉铆接元模型的变形云图，其中圆内部分节点的变形即为向壁板整体铆接仿真模型中进行映射加载的局部变形场。

图4 单钉铆接元模型Fig.4 Single-rivet FE model

图5 单钉铆接元模型铆接变形分布Fig.5 Deformation distribution of single-rivet FE model

3.2 壁板整体铆接仿真模型

（1）几何模型。

图6所示为壁板铆接变形仿真模型。模型包含两部分：一部分是采用壳单元建模的壁板（绿色）；另一部分是采用实体单元建模的壁板上铆钉孔周围区域（黄色）。壳单元厚度设置为2mm。

图6 壁板铆接变形计算数值模型Fig.6 FE model of panel for riveting deformation calculation

（2）材料设置。

壁板材料参数如表1所示。

表1 铝合金7075-T651材料参数

（3）接触设置。

根据主从关系设置模型中各零件的接触关系，壁板零件间的摩擦系数设置为0.19。

（4）连接定义。

模型中实体单元和壳单元间连接则通过多点约束方程实现，如图7所示。

图7 体元与壳元连接定义Fig.7 Definition of solid-shell combination

（5）网格划分。

对于钉孔周围实体部分，铆接变形时钉孔受到铆钉挤压产生塑性变形，为了提高仿真精度，采用C3D8R单元，并对钉孔周围区域及铆钉进行网格细化，而对于远离钉孔的壳体部分，采用S4R单元，为提高计算效率，可适当粗化网格。

（6）边界条件。

根据实际铆接约束情况，限制壁板两侧端面的所有自由度，如图8所示。

（7）局部变形场加载。

壁板自动钻铆时，相同规格铆钉的铆接工艺参数相同，在钉孔周围区域产生的局部变形场也基本相同。所以在确定的铆钉规格下，对壁板进行成百上千个铆钉安装时，壁板增量铆接变形可从依次仿真每一个铆钉的铆接过程转化为不考虑铆钉的铆接过程，只需将此规格铆钉在单钉铆接仿真结果中所获得的局部变形场作为加载条件，依次加载到壁板上各个铆钉孔周围区域对应节点上即可。

图8 壁板边界约束条件设置Fig.8 Boundary condition of panel

4 壁板铆接变形数值计算模型验证

为验证本文所建立的“体-壳”连接壁板铆接仿真模型的可行性，将其仿真结果与对应的全三维实体模型的仿真结果进行比对，两模型的变形云图如图9所示。分别从两个模型中线上等距取了6个测点进行比较，测点位置见图9，表2所列为“体-壳”连接模型和三维实体模型中测点的法向变形值，图10为测点对比结果。

图9 “体-壳”连接模型与实体模型变形分布对比Fig.9 Comparison of riveting deformation distribution between solid-shell combination model and 3D solid model

表2 “体-壳”连接模型与三维实体模型测点法向变形值

基于遗传算法的铆接顺序优化

铆接顺序对铆接变形有着很大影响，为了保证壁板装配精度，需合理规划铆接顺序，以减小铆接变形。本文采用遗传算法，以壁板铆接变形最小为目标函数，进行铆接顺序优化。其数学模型为：

式中，F为优化目标，u为壁板铆接后的最大变形。

1 遗传算法

遗传算法是启发式算法中被广泛应用的一种算法，它是一类借鉴生物界的达尔文进化论中自然选择、适者生存、优胜劣汰进化规律，演化而来的随机化搜索最优解的方法[19]。本文基于遗传算法的铆接顺序优化流程如图11所示。

（1）编码和设置种群。

遗传算法的编码方式有很多种，对于壁板铆接过程，每一个铆钉孔只进行一次铆接操作，而铆钉数量很多，若采用二进制编码，其解码过程非常复杂，因此可对每一个钉孔标记唯一的数字编号。例如：有7个铆钉需要进行铆接，对铆钉孔编号进行排序即可代表铆接顺序，则[1 2 3 4 5 6 7 ]即为一个铆接顺序。设置种群规模为m，按钉孔编号排列随机生成铆接顺序。为更加贴合实际生产，本文规定每种顺序均从第一个铆钉编号开始。

（2）适应度函数。

在实际铆接过程中，需要使得壁板最终铆接变形越小越好，即求解最小值问题，将仿真得到的壁板变形值作为染色体的适应值，建立适应度函数为：

图10 “体-壳”连接模型与三维实体模型测点变形值对比Fig.10 Riveting deformation comparison of measuring points between solid-shell combination model and 3D solid model

图11 铆接顺序的遗传算法优化流程Fig.11 Genetic algorithm optimization process of riveting sequence

（3）选择算子。

选择过程采用轮盘赌方式实现，对于第 t代种群，A（t）={at,1，at,2，at,3，···，at,m}，m 为第 t代种群个数，其所对应适应度值为 { f（at,1），f（at,2），···，f（at,m）}，顺序优化的目标是使铆接最大变形值最小，即适应值越大，被选择的概率越大，确定选择概率为：

（4）交叉匹配。

本文采用部分匹配交叉法，它要求随机选取两个交叉点，来确定一个匹配区域，然后根据两个父个体中交叉点之间区域给出的映射关系生成两个子个体。

下面举例说明交叉操作：对于7个铆接点，生成两个铆接顺序的父个体P1、P2，随机地选择两个交叉点“|”：

P1：[ 1 2 3 | 4 5 | 7 6 ]

P2：[ 2 4 6 | 7 5 | 3 1 ]

首先，两个交叉点之间段交换，得到：

Q1：[ 1 2 3 | 7 5 | 7 6 ]

Q2：[ 2 4 6 | 4 5 | 3 1 ]

并获得中间段映射关系为：

7-4，5-5

为了解决编码重复问题，按映射关系变更父个体P1重复编码7，于是得到子个体 Q1’：[ 1 2 3 | 7 5 | 4 6 ]，同理可得到子个体Q2’：[ 2 4 6 | 7 5 |3 1 ]。

（5）变异。

变异，即交换两个随机位置上的基因。例如：对于7个铆钉的铆接，其中一种顺序为[1 3 5 7 6 4 2]，将第4和第6节点反转，得到[1 3 5 4 6 7 2]的顺序。

（6）生成新种群。

通过遗传操作生成新的铆接顺序种群，分别对每种铆接顺序进行有限元分析，得到对应的变形值，并通过适应度函数计算新种群的适应值，保存历史和全局的最优铆接顺序及其对应的适应值，对新生成的铆接顺序种群继续采用遗传操作，得到下一代种群。

（7）停止准则。

设置最大迭代次数，终止遗传计算，根据所保存的历史和全局最优铆接顺序及其对应的适应值，确定优化目标，并得到相应的最优铆接顺序。

2 实例分析

铆接顺序优化对象及铆钉孔编号如图12所示，建模过程如上文所述。

图12 铆接顺序优化实例Fig.12 Case of riveting sequence optimization

设置初始种群为8，最大遗传代数为30代，变异概率0.05，交叉概率0.9，在MATLAB中编写遗传算法的代码并运行。图13为优化过程，可以看出，经过21次迭代结果收敛。所得优化结果为：铆接最大变形为0.1769mm，对应铆接顺序为1-3-5-6-4-2-7。而按实际生产中常见的顺序铆接1-2-3-4-5-6-7进行仿真计算，得到的铆接最大变形为0.2344mm；按1-3-5-7-6-4-2的铆接顺序进行仿真，所得铆接最大变形为0.2033mm，可见采用遗传算法对铆接顺序进行优化，能有效减小壁板增量铆接时的最大变形。

本文采用的算例铆钉数目较少，对于大尺寸壁板铆接顺序规划问题，可先对壁板进行铆接区域划分，然后利用本文方法进行区域间的铆接顺序优化，从而实现铆接变形的控制。另外，壁板多钉铆接顺序规划和试验验证将是下一步研究的重点。

图13 铆接顺序优化过程Fig.13 Process of riveting sequence optimization

结论

本文针对壁板铆接过程中整体变形的求解问题，提出了以单钉铆接变形元模型为基础，“体-壳”连接为手段，通过“局部-整体”映射加载方法，实现了壁板铆接变形的快速计算。

（1）分析了壁板自动钻铆的工艺流程，给出了铆接过程中影响壁板变形的主要工艺因素。

（2）采用“体-壳”连接建模方法，提出了壁板铆接数值计算的“局部-整体”映射加载模型，并通过算例验证了该模型的可行性。

（3）以铆接变形最大值最小为目标，建立了以遗传算法为核心的铆接顺序优化模型，实现了壁板铆接顺序的优化。

[1]常正平, 王仲奇, 王斌斌, 等. 基于镦头不均匀变形的压铆力建模[J]. 航空学报,2016, 37(7): 2312-2320.

CHANG Zhengping, WANG Zhongqi,WANG Binbin, et al. Riveting force computation model based on formed head inhomogeneous deformation[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(7): 2312-2320.

[2]JEFFREYS D, LEANEY P G.Dimensional control as an integral part of nextgeneration aircraft development[J]. Proceedings of the Institution of Mechanical Engineers Part B: Journal of Engineering Manufacture, 2000,214(9): 831-835.

[3]MÜLLER R P G. An experimental and analytical investigation on the fatigue behavior of fuselage riveted lap joints[J]. Aerospace Engineering, 1995, 376(3): 201-207.

[4]RIJCK J J M D, HOMAN J J, SCHIJVE J, et al. The driven rivet head dimensions as an indication of the fatigue performance of aircraft lap joints[J]. International Journal of Fatigue,2007, 29(12): 2208-2218.

[5]CHERAGHI S H, KRISHNAN K,BAJRACHARYA B. Effect of variations in the riveting process on the quality of riveted joints[J]. The International Journal of Advanced Manufacturing Technology, 2008, 39(11/12):1144-1155.

[6]DA CUNHA F, FIGUEIRA J, DE BARROS M. Methodology to capture induced strains on riveting process of aerospace structures[C]//Proceeding of SAE Brasil 2010 Congress and Exhibit. São Paulo, 2010: 8-18.

[7]ZHANG K F, CHENG H, LI Y.Riveting process modeling and simulating for deformation analysis of aircraft’s thinwalled sheet-metal parts[J]. Chinese Journal of Aeronautics, 2011, 24(3): 369-377.

[8]牟伟强, 李原, 张开富, 等. 无头铆钉压铆力数学建模与仿真分析[J]. 西北工业大学学报, 2010, 28(5): 742-747.

MOU Weiqiang, LI Yuan, ZHANG Kaifu,et al. Mathematical modeling for and simulation analysis of flush rivet pressing force[J]. Journal of Northwestern Polytechnical University, 2010,28(5): 742-747.

[9]刘平. 铆接变形及其有限元分析[D].西安: 西北工业大学, 2007.

LIU Ping. Rivet deformation and FE simulation[D]. Xi’an: Northwestern Polytechnical University, 2007.

[10]冯晓旻, 谢兰生. 铆接变形的有限元分析[J]. 机械制造与自动化, 2009, 38(2):62-63.

FENG Xiaomin, XIE Lansheng. Rivet deformation and finite element analysis[J].Machinery Manufacturing & Automation, 2009,38(2): 62-63.

[11]AMAN F, CHERAGHI S H,KRISHNAN K K, et al. Study of the impact of riveting sequence, rivet pitch, and gap between sheets on the quality of riveted lap joints using finite element method[J]. The International Journal of Advanced Manufacturing Technology,2013, 67(1-4): 545-562.

[12]王仲奇, 常正平, 郭飞燕, 等. 大型飞机壁板无头铆钉干涉连接技术[J]. 航空制造技术, 2015(4): 34-39.

WANG Zhongqi, CHANG Zhengping,GUO Feiyan, et al. Interference fit technology of slug rivet in large aircraft panel[J]. Aeronautical Manufacturing Technology, 2015(4): 34-39.

[13]WANG H. Riveting sequence study of horizontal stabilizer assembly using finiteelement analysis and riveting equivalent unit[J].Journal of Aerospace Engineering, 2014, 27(6):04014040-1-04014040-8.

[14]李娅娜, 兆文忠, 崔晓芳. 基于“局部-整体”映射法的焊接装配变形数值仿真[J]. 焊接技术, 2008, 37(1): 43-46.

LI Ya’na, ZHAO Wenzhong, CUI Xiaofang.Welding assembly deformation numerical simulation based on "local- global" method[J].Welding Technology, 2008, 37(1): 43-46.

[15]李江飞. 薄壁多焊缝构件结构优化及焊接变形的数值模拟[D]. 石家庄: 石家庄铁道大学, 2014.

LI Jiangfei. Numerical simulation of structural optimization and welding deformation in thin-walled multi-welds component[D].Shijiazhuang: Shijiazhuang Tiedao University,2014.

[16]周广涛, 刘雪松, 闫德俊, 等. 顶板焊接顺序优化减小焊接变形的预测[J]. 焊接学报, 2009, 30(9): 109-112.

ZHOU Guangtao, LIU Xuesong, YAN Dejun, et al. Prediction for welding deformation reducing by welding sequence optimization of upper plate[J]. Transactions of the China Welding Institution, 2009, 30(9): 109-112.

[17]MASTERS I, FAN X, ROY R, et al.Modelling distortion induced in an assembly by the self-piercing rivet process[J]. Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture, 2012,226(2): 300-312.

[18]RANS C, STRAZNICKY P V,ALDERLIESTEN R. Riveting process induced residual stresses around solid rivets in mechanical joints[J]. Journal of Aircraft, 2007, 44(1): 323-329.

[19]李敏强, 寇纪淞, 林丹. 遗传算法的基本理论与应用[M]. 北京: 科学出版社,2002.

LI Minqiang, KOU Jisong, LIN Dan.The basic theory and application of genetic algorithm[M]. Beijing: Science Press, 2002.