0 引言
作为飞机气动外形的重要组成部分，壁板的装配精度直接影响飞机的使用寿命。壁板具有刚度小、尺寸大等特点，在铆接装配时钉孔周围会产生凹陷，远离钉孔中心处会发生翘曲变形[1]。壁板铆接时会涉及成千上万个铆钉连接，随着铆接的不断进行会使壁板发生整体变形，进而造成部件装配过程中的强迫装配。强迫装配附加的装配应力会增大连接处的场应力、降低结构的连接强度以及缩短结构的疲劳寿命[2]。近年来，随着飞机性能的不断提高，对机体结构的各项性能提出了更高要求，为保证飞机在高速下的性能指标和寿命，对铆接变形的研究至关重要。
目前对铆接的研究主要包括如下两个方面。一是有关单钉铆接变形和质量的研究：朱春东等[3]研究了电磁铆接时放点次数与铆钉头部变形量的关系，得出了低电压适合连接复合材料的结论；黄志超等[4]采用摆锤铆接和直压铆接对接头的应力应变性能进行了测试；FOX等[5]利用同步辐射和中子衍射测量了铆接后的三向应力状态，研究了孔周的应力分布；楼明等[6]研究了自冲铆接中模钉体积比对不同材料铆接接头的影响。二是有关多钉铆接引起的整体变形分析的研究：WANG[7]基于建立的铆接力学等效单元优化了飞机水平安定面装配过程中的铆接顺序，减小了铆接变形；吴少帅[8]对基于温度等效模型的多孔铆接变形进行了有限元分析，显著提高了计算效率；LIU等[9]基于等效法，用外部力矩引起的弯曲变形来等效铆接弯曲变形，并通过实验验证了方法的有效性；MASTERS等[10]基于局部-全局方法，通过平板装配和顶帽装配的对比，研究了自冲铆接时铆钉变形对整体装配的影响；NI等[11]基于局部变形到整体的加载方法，提出了一种天线结构的整体装配偏差分析构架；常正平等[12]提出了基于局部-整体的“体-壳”连接方法，使得多钉铆接的效率得到了显著提高，实现了对壁板变形的有效控制。
由上述研究可以看出，对多钉铆接的变形研究大多采用局部-整体的映射方法，这主要是考虑到多钉铆接分析会带来巨大计算量的限制。具体的局部-整体映射方法又分两种思路：一是采用简化的模型进行计算，该思路的优势是计算效率高，缺点是因实际结构和铆接局部变形的不准确而带来较大的计算误差，该方法更适合定性分析；二是先进行局部模型分析，再将分析的应力或位移状态映射到整体结构的相应位置，该方式的计算效率较低，但精度较高，定量分析效果更好。然而，目前有关多钉铆接使用的局部-整体映射模型中的映射区域大小以及映射区域对计算精度及效率的影响方面的研究报道较少。本文针对壁板铆接变形分析采用的局部-整体映射方法进行了深入研究，提出了基于局部位移场的分层映射方法。该方法通过对映射区域内节点进行分层，以各层节点为加载对象，有效地减小了针对单个节点为对象时的疲劳强度。通过与三维动态模型进行对比，研究了映射区域范围对计算精度与效率的影响。通过分析单排三钉预测模型与三维动态模型的对比结果验证了所提方法的有效性。
1 局部位移场分层映射方法
1.1 局部位移场数据分析
铆接变形主要是由局部应力和应变的变化造成的，局部变形的逐渐积累导致了整体模型变形[13]。铆接过程中，钉孔周围的区域处于复杂的三向应力状态，而离钉孔中心2.5R(R为钉孔半径)以外的部分处于平面应力状态[14]，塑性变形主要集中在相对小的范围内，则塑性区域外的变形可视为弹性的，此为局部映射方法的前提。基于上述思想，铆接的塑性变形和应力模拟可从整个结构变形计算中分离出来，单独地进行局部的变形分析，利用映射技术将局部模型的计算结果转换到整体模型中，最终预测整体构件的变形量[15]。
局部场数据是实现批量铆接的基础，采用映射技术使局部场数据最大程度地映射到整体模型中，首先要研究局部场数据的分布规律。本文以沉头铆钉(牌号为NAS1097AD6-6)为例，建立了单钉铆接有限元计算模型。具体过程如下。
(1)建立单钉模型。为精确地获得钉孔周围的位移场数据，本文采用三维动态模型进行研究。图1所示为依据某飞机壁板尺寸建立的单钉有限元模型，其参数如下：铆钉直径4.78 mm、长度9.78 mm，蒙皮厚度2.2 mm，长桁厚度1.6 mm，钉孔直径4.89 mm。
图1 单钉铆接示意图
Fig.1 Schematic diagram of single rivet riveting
(2)材料设置。仿真模型材料采用新型材料铝锂合金：蒙皮材料为2060-T8，长桁材料为2099-T83，铆钉材料为2117-T4，将铆模(顶铁)设为刚体，材料具体属性参数见表1。
表1 模型材料参数
Tab.1 Material parameters
材料弹性模量E(GPa)泊松比μ屈服强度σs(MPa)密度ρ(kg/m3)2060-T8720.334412 7002117-T4710.331652 7002099-T83780.285047 800
(3)接触设置。根据主从关系设置模型中各零件的接触关系，壁板与长桁间的摩擦因数为0.2。
(4)边界条件。铆接分析是复杂的非线性接触问题，为减小收敛难度，限制模型4个角的移动自由度。
(5)网格划分。为提高仿真精度，对钉孔周围的部分进行区域划分，并对钉孔部分进行细化，网格采用八节点减缩积分单元C3D8R，在远离钉孔处适当增大网格尺寸。
(6)结果与分析。考虑到被连接件的变形主要发生在沿孔的径向与轴向，因此需对径向与轴向的变形进行分析。为利于对局部场数据进行分析，对处于同一层的节点进行径向划分：以孔轴线为中心，沿径向将由处于同一距离的节点连接而成的曲线命名为曲线k(k=1,2,…,n)，在曲线k上等距选取8个节点作为测量点，分析测量点之间的位移，如图2所示，其中Dr、Da、Dt分别表示径向、轴向和切向。
图2 节点径向分层示意图
Fig.2 Schematic diagram of node radial stratification
单钉模型的应力和轴向位移云图见图3，可以看出，铆接时塑性变形影响区域发生在距钉孔中心2.5R范围内，这说明局部映射方法的前提[14](即采用局部-整体映射时，塑性变形发生在局部区域，在塑性区域外的变形可视为弹性)是正确的。图4为在钉孔中心创建的局部柱坐标系下测量点在径向与轴向的位移变化曲线，可以看出，对于处于同一条曲线上的测量点，其径向与轴向的位移基本稳定在一个定值，因此，可近似地取曲线k上所有节点的平均值作为单钉模型中的局部场数据。
(a)单钉模型应力云图
(b)单钉模型轴向位移云图
图3 单钉模型应力及轴向位移云图
Fig.3 Stress and axial displacement cloud pictures of single-rivet model
(a)径向位移曲线
(b)轴向位移曲线
图4 测量点位移曲线
Fig.4 Displacement curves of measuring points
1.2 位移场分层映射方法
位移场分层映射方法是在“局部-整体”映射思路下，将单钉模型计算后的局部场数据加载到整体模型中;将模型中待映射区域的节点进行分层操作，以各层节点为加载对象，调用局部场数据映射到相应的节点层上。考虑到单曲度蒙皮曲率较小，在钉孔局部区域可近似视为平面，所以本文采用的分层映射法的思想是通过引入一个包围盒并改变其尺寸来对铆接主要影响区域内的节点进行分层操作，以保证局部场能够准确地映射到对应的节点上，如图5所示，具体操作如下。
图5 局部位移场分层映射过程
Fig.5 Layed mapping process of local displacement field
(1)依据壁板中零件的厚度尺寸建立单钉铆接动态模型，并划分铆接影响区域。
(2)单钉模型的动态模拟仿真完成后，提取局部位移场。
(3)借助离线程序获取壁板MBD(基于三维模型的定义)的铆接孔位信息(钻铆点坐标及法矢)。以钻铆点为中心创建包围盒，获取被铆接件的映射区域，其尺寸与单钉铆接动态模型中划分的影响区域相同。考虑到壁板表面具有一定的曲率，直接对模型进行网格划分难以保证单钉模型网格与壁板网格对应，因此在有限元仿真软件ABAQUS中利用Python语言进行处理。具体过程如下：①以钻铆点坐标为起始点，从模型所有的几何边中获取距离起始点最短的几何边，并以该几何边的长度为基准获取孔的另一条边线，依据两条边线上的坐标与孔径大小来获取钻铆点处的厚度尺寸t；②以钻铆信息(坐标及法矢)为输入，根据“点-法矢”方法创建钻铆点的切平面，记为plane1；③采用偏移平面的方法，沿法矢负方向平移平面plane1，平移距离为t/2，并创建平面plane2(plane2位于模型内部)；④在平面plane2上即可创建分割圆草图，为完全划分出所有映射区域，先后沿法矢正方向和法矢负方向对映射区域进行分割；⑤循环利用不同的钻铆信息即可完成壁板的区域划分，并获取与单钉模型一致的映射区域。
(4)完成映射区域划分后，使用定比分点坐标公式[16]沿法矢负方向对铆接点坐标进行变换，则有
(1)
(2)
(3)
式中,t为被铆接件厚度；(x0,y0,z0)为铆接点坐标；a、b、c为铆接孔法矢正方向余弦；N为厚度方向上的等分数；(xp,yp,zp)为沿着铆接孔法矢负方向第p(p=1,2,…，N)次变换后的坐标。
(5) 以步骤(4)中经变换得到的坐标(xp,yp,zp)为中心创建圆柱包围盒,对铆接影响区域进行轴向分层操作，并创建面节点集合。
(6)以变换后的坐标(xp,yp,zp)为起点，将l记为该起点与步骤(5)中创建的面节点集合中的节点之间的距离。根据下式创建相应的环节点集合Set-CirNode并进行径向分层操作：
dq—d<Δ
(4)
式中，L为铆接影响区外径与铆钉孔半径的差值；M为径向上的等分数；dq为第q(q=1,2，…，M)次径向分层的等分长度；Δ为设定容差。
(7)调用局部场数据库中的位移数据，映射到创建的相应节点集合。
2 映射方法的精度与效率
2.1 映射区域范围对分析精度的影响
为验证基于分层法的“局部-整体”映射模型铆接变形的计算精度，以局部位移场数据为边界条件，将其仿真结果与三维动态模型结果进行对比。由图3可知，铆接时塑性变形主要发生在铆接孔局部区域，为避免因映射区域过大而造成的相邻钉孔映射区域重叠的情况，映射区域范围选取1.5R～4.0R(壁板钉孔间距为10R)，区域间隔为0.5R。铆接件尺寸为60 mm×60 mm×2.2 mm。在铆接模拟过程中对被连接件的4个角进行移动自由度约束。为观察铆接完成后铆接件沿孔径向的应力变化以及边缘处的轴向位移变化，沿孔径向等距选取5个节点(以圆点表示),在铆接件边缘等距选取5个节点(以交叉点表示)，测量点位置及顺序示意图见图6。
图6 铆接件测量点示意图
Fig.6 Schematic diagram of measuring points
表2所示为各测量点的应力值，可以看出，测量的应力值与三维动态模型仿真值的差值有正有负。为方便观察各测量点的变化情况，图7中的应力误差值为各测量点应力值与三维动态模型仿真值差值的绝对值。
表2 测量点应力
Tab.2 Measuring point stress MPa
三维动态模型测量点1测量点2测量点3测量点4测量点5464.22199.1690.0146.2225.13映射范围1.5R465.35177.4979.1341.5624.052.0R465.34195.2587.1044.8624.712.5R465.35200.9989.9846.1625.113.0R465.35201.9789.9146.1925.123.5R465.34201.9689.9546.2325.324.0R465.35201.9794.5046.2325.54
图7 铆接件测量点应力误差绝对值变化曲线
Fig.7 Curves of absolute value of stress difference for riveted component points
由表2和图7可知，对于不同测量点，当映射区域处于三向应力状态(1.5R～2.5R)时，随着映射区域的增大，应力误差的绝对值会逐渐减小；当映射区域由三向应力状态过渡到平面应力状态(2.5R～3.5R)时,随着映射区域的增大，应力误差的绝对值基本不变，并稳定在某个定值；当映射区域完全处于平面应力状态(3.5R～4.0R)时，应力误差绝对值的最大值出现在测量点3，为5 MPa。
表3所示为各测量点的轴向位移值，同理可以看出，测量的轴向位移值与三维动态模型仿真值的差值有正有负，因此图8中的轴向位移误差值为上述两者差值的绝对值。
表3 测量点轴向位移
Tab.3 Displacement of measuring points in axial direction mm
三维动态模型测量点1测量点2测量点3测量点4测量点5-0.020 6-0.018 9-0.016 2-0.012 5-0.008 1映射范围1.5R-0.011 3-0.010 4-0.008 9-0.006 9-0.004 52.0R-0.018 3-0.016 8-0.014 3-0.011 1-0.007 22.5R-0.020 4-0.018 7-0.015 9-0.012 3-0.007 93.0R-0.020 4-0.018 7-0.016 0-0.012 4-0.008 03.5R-0.020 3-0.018 7-0.016 1-0.012 4-0.008 14.0R-0.020 2-0.018 7-0.016 1-0.012 5-0.008 1
(a)轴向位移误差绝对值变化曲线
(b)轴向位移相对误差曲线
图8 铆接件测量点轴向位移误差变化曲线
Fig.8 Curves of axial displacement error for riveted component points
由表3和图8可知，对于不同的测量点，当映射区域处于三向应力状态(1.5R～2.5R)时，随着映射区域的增大，轴向位移误差的绝对值逐渐减小，在1.5R～2.0R映射范围内，相对误差大于10%，不满足工程应用要求；当映射区域由三向应力状态过渡到平面应力状态(2.5R～3.5R)时，随着映射区域的增大，轴向位移误差的绝对值趋于稳定，误差的绝对值约为4.7×10-4 mm，相对误差为2.2%；当映射区域完全处于平面应力状态(3.5R～4.0R)时，位移误差的绝对值虽有增加，但相对误差最大约为4.8%。
2.2 映射区域范围对分析效率的影响
为验证基于分层法的局部-整体映射模型的计算效率相对于传统三维动态模型的计算效率有所提高，将所提模型的仿真分析计算时间与三维动态模型的计算时间进行对比，并选取最优的映射区域范围。仿真完成后的计算时间见表4。本文使用的计算机配置为Intel(R) Core(TM) i5-6500CPU,3.2 GHz，内存为32 GB。
表4 计算时间对比
Tab.4 Comparison of computational time min
三维动态模型不考虑几何非线性考虑几何非线性45.46映射范围1.5R4.664.832.0R3.613.602.5R3.603.333.0R3.473.603.5R4.404.434.0R3.603.91
由表4可知：①利用映射模型分析的时间相对于利用三维动态模型的分析时间有大幅度缩短，两者之间相差9.0～13.5倍；②两种模型在映射区域范围为2R～3R时的计算时间相差较小，稳定在3.5 min左右，且无论是否考虑几何非线性，在同一映射范围内的计算时间并无明显的差异。
综上可知，考虑到映射模型的计算精度和效率，结合工程应用中的铆钉孔边距大小(最小钉孔边距为4R)，当映射区域为2.5R时，与三维动态仿真结果相比，所提模型的误差较小、计算时间更短。当映射区域大于2.5R时，沿孔径向需要进一步分层，这会增加前处理时间与求解计算成本，进而会增加总体的仿真时间，因此选取映射区域为2.5R时，在该区域范围内既能保证计算精度又能提高分析效率。
3 实例分析
为进一步验证所提方法的计算精度与效率，在不考虑铆接顺序对变形的影响下，保证相同的铆接顺序，以单排三钉模型、映射区域为2.5R为例进行分析，将其分析结果与三维动态模型的仿真结果进行对比。
3.1 几何模型
壁板模型由一块蒙皮与一根长桁组成，其具体尺寸为100 mm×60 mm×2.2 mm，铆钉类型及材料属性见表1，在铆接过程中，限制壁板两端4个角的移动自由度。为使模型映射区域中节点与单钉模型中节点相对应，采用分割技术设置与单钉模型相同的局部网格数量。铆接顺序选择为沿长桁方向依次铆接。在壁板侧边选择7个测量点，如图9所示。比较映射模型与三维动态模型的结果，两种模型铆接完成后的应力与轴向位移云图分别见图10和图11。图9中的7个测量点在两种模型下的轴向位移对比曲线见图12。两种模型下计算时间的对比见表5。
图9 测量点示意图
Fig.9 Schematic diagram of measuring points
3.2 结果讨论
(a)三维动态模型应力分布
(b)映射模型应力分布
图10 铆接件应力分布云图
Fig.10 Cloud picture of stress distribution for riveted component
(a)三维动态模型轴向位移分布
(b)映射模型轴向位移分布
图11 铆接件轴向位移分布云图
Fig.11 Cloud picture of axial displacement distribution for riveted component
(1)在应力分布方面，由图10可以看出，使用的映射模型与三维动态模型得到的铆接应力分布基本相同，最大应力相差21.5 MPa。
(2)在轴向位移分布方面，当模型发生翘曲时，由图11中轴向位移云图对比可知，使用三维动态模型得到的最大轴向位移为1.485×10-3 mm,使用映射模型得到的最大轴向位移为1.506×10-3,两者相对误差为1.4%。结合图12可知，两种模型下测量点的位移变化趋势基本一致，这表明使用映射模型可得到与直接使用三维动态模型相同的结果。
(3)由表5可知，在三钉模型中，同一计算机配置下，使用三维动态模型时每个铆钉的平均计算时间为67.1 min，而使用本文的映射模型时每个铆钉的平均计算时间为3.9 min，明显提高了计算效率。
图12 三钉模型测量点位移比较
Fig.12 Comparison of displacement of measuring point for three rivets model
表5 计算时间对比
Tab.5 Comparison of computational time
min
完成时间截至第一个钉截至第二个钉截至第三个钉平均时间三维动态模型93.1159.7201.467.1映射模型6.49.311.73.9
4 结论
(1)提出了基于局部位移场的分层映射方法，以各层节点为加载对象，有效地减小了以单个节点为研究对象的劳动强度。通过对不同参数结果进行对比分析，得出了优化的映射区域范围。
(2)映射区域范围在一定程度上影响了分析结果的精度，过大或过小的范围都会影响应力分析的精度。当映射区域范围在2.5R～3.5R时(R为钉孔半径)，与三维动态模型相比被连接件应力误差基本稳定。
(3)不同的映射区域范围影响了被连接件的变形精度，当映射区域范围在2.5R～3.5R时，被连接件轴向位移相对误差最大为2.2%。
(4)在三钉实例对比模型中，同一计算机配置下，将三维动态模型与分层映射模型的计算时间进行比较，当映射区域为2.5R时，单钉平均仿真时间从67.1 min缩短至3.9 min，明显提高了计算效率。
[1] 刘平，张开富，李原.铆接变形及其有限元分析[J].航空制造技术，2008(22):94-97.
LIU Ping, ZHANG Kaifu, LI Yuan. Riveting Defor-mation and Its Finite Element Analysis[J]. Journal of Aerospace Manufacturing Technology， 2008(22):94-97.
[2] 钱一彬.基于试验的强迫装配对飞机金属结构疲劳品质的影响[J].中国科技信息，2017(8):28-29.
QIAN Yibin. The Influence of Forced Assembly on Fatigue Quality of Aircraft Metal Structure Based on Test[J].Chinese Science and Technology Infor-mation,2017(8):28-29.
[3] 朱春东，刘卫华，黄尚宇. 低电压电磁铆接成形工艺的研究[J]. 中国机械工程，2005,16(16)：1465-1467.
ZHU Chundong，LIU Weihua，HUANG Shangyu. Research on Low Electric Voltage Electromagnetism Rivet Joint Heading Process[J].China Mechanical Engineering， 2005,16(16)：1465-1467.
[4] 黄志超，刘晓坤，占金青，等.铝合金板料摆碾铆接与直压铆接比较分析[J]. 中国机械工程，2013,24(9):1233-1238.
HUANG Zhichao， LIU Xiaokun， ZHAN Jinqing,et al. Comparative Study on Rolling Riveting and Pressing Riveting[J]. China Mechanical Engineering， 2013,24(9):1233-1238.
[5] FOX M L, WHITHERS P J.Residual Stress in and around Electro-magnetically Installed Rivets Measured Using Synchrotron and Neutron Diffraction[J]. Journal of Neutron Research,2007,15(3/4):215-223.
[6] 楼明，李永兵，黄舒彦，等. 模钉体积比对异种金属自冲铆接接头成形性能影响分析[J]. 中国机械工程，2009，20(15)：1873-1876.
LOU Ming， LI Yongbing， HUANG Shuyan，et al. Influence of Die-rivet Volume Ratio on Forming Performance of Self-piercing Riveting Joints of Dissimilar Material[J]. China Mechanical Engineering，2009,20(15)：1873-1876.
[7] WANG H. Riveting Sequence Study of Horizontal Stabilizer Assembly Using Finite-element Analysis and Riveting Equivalent Unit[J]. Journal of Aerospace Engineering,2014, 27(6):04014040.
[8] 吴少帅. 飞机壁板自动钻铆过程变形分析[D].杭州：浙江大学，2017.
WU Shaoshuai. Deformation Analysis of the Auto-matic Drilling and Riveting Process of Aircraft Wall Board[D]. Hangzhou: Zhejiang Univetsity,2017.
[9] LIU G, HUAN H, KE Y. Study on Analysis and Prediction of Riveting Assembly Variation of Aircraft Fuselage Panel[J]. International Journal of Advanced Manufacturing Technology, 2014, 75(5/8):991-1003.
[10] MASTERS I，FAN X，ROY R，et al. Modelling Distortion Induced in Assembly by the Self-piercing Rivet Process[J]. Proceedings of the Institution of Mechanical Engineers， Part B： Journal of Engineering Manufacture,2011, 226(2):300-312.
[11] NI J， TANG W， XING Y， et al. A Local-to-global Dimensional Error Calculation Framework for the Riveted Assembly Using Finite-element Analysis[J]. Journal of Manufacturing Science & Engineering, 2016, 138(3):1-13.
[12] 常正平，王仲奇，张津铭，等.采用局部-整体映射模型的壁板铆接变形预测[J]. 西安交通大学学报， 2018，52(1):115-122.
CHANG Zhengping， WANG Zhongqi， ZHANG Jinming， et al. Prediction of Riveting Deformation of Aircraft Panel Using Local-global Mapping Approach[J].Journal of Xi’an Jiaotong University, 2018,52(1):115-122.
[13] 张开富，杨殿国，李原,等.壁板自动钻铆行为及变形分析技术综述[J].航空制造技术，2010(23):67-72.
ZHANG Kaifu， YANG Dianguo， LI Yuan， et al. Deformation Analysis and Automatic Riveting for Wing Panel [J]. Journal of Aerospace Manufacturing Technology, 2010(23):67-72.
[14] RANS C, STRAZNICKY P V, ALDERLIESTEN R. Riveting Process Induced Residual Stresses around Solid Rivets in Mechanical Joints[J]. Journal of Aircraft, 2012,44(1):323-329.
[15] 张乐乐，李培，崔进. 基于整体-局部方法的焊接构架承载能力的分析[J]. 铁道学报，2011，33(8):33-38.
ZHANG Lele，LI Pei，CUI Jin. Analysis on Load-carrying Capacity of Welded Bogie Frame Based on Global-local Method[J]. Journal of the China Ra-ilway Society, 2011, 33(8):33-38.
[16] 张冠宇，赵中建. 空间线段定比分点公式的另一种表达形式[J]. 天中学刊，2006，21(2):93-94.
ZHANG Guanyu，ZHAO Zhongjian. Another Form of Expression of the Formula of the Space Line Segment[J]. Journal of Tian Zhong,2006, 21(2):93-94.
