基于拉丁超立方抽样的薄板装配误差分析

邓乾旺 文 文

湖南大学汽车车身先进设计制造国家重点实验室，长沙，410082

0 引言

轿车白车身装配的尺寸偏差对整车的性能有着极大的影响，误差分析是装配尺寸偏差控制的有效手段。

基于刚性体假设的薄板装配偏差分析技术已经得到深入研究，许多学者在装配偏差建模方法上开展了卓有成效的工作。文献［1－2］采用尺寸链模型直接线性化方法对薄板类零件进行2D、3D系统的装配公差分析。文献［3］对零件进行确定性定位分析，利用装配约束函数的雅可比矩阵，求解零件偏差对零件定位偏差的敏感度矩阵。

由于薄板零件的特殊性，基于刚体假设的误差分析理论不适用于薄板零件。近年来，利用有限元分析工具对薄板装配误差进行仿真已成为国际上研究的前沿课题。文献［4］提出的影响系数法将有限元方法与统计方法结合，研究零件偏差与装配偏差的关系，并建立了两者之间的影响系数矩阵，大大减小了有限元仿真次数，但由于影响系数法的理论局限性，该方法未考虑焊接时的零件间接触等非线性问题。文献［5］对影响系数法进行改进，在使用影响系数法的同时考虑零件间的接触问题，提高了仿真精度，但因此也降低了计算效率。

目前，各种统计工具的运用研究以及模型的简化处理使得有限元法在薄板零件装配误差仿真中的应用有了很大的进步［6－7］，但车体零件众多，各组件的层次关系及装配过程复杂，从而使得精确的误差分析存在较大困难。本文将拉丁超立方抽样（Latin hypercube sampling）代替以往仿真时的蒙特卡罗直接抽样（direct Monte Carlo sampling），并在装配模型中考虑装配时可能发生的零件与零件、零件与焊枪间的接触作用，构造出焊接完成后零件的变形对初始偏差的响应面模型，并比较两种抽样方法的收敛速度与精度。

1 薄板零件装配误差分析过程

典型车体薄板零件焊接一般包括以下4个过程：

（1）定位。将零件放置在夹具的定位元件上，由于零件存在误差，故零件无法与夹具完全贴合。薄板零件常采用N－2－1（N＞3）定位原则进行定位；进行有限元仿真时，采用在定位点对应位置的节点施加位移约束的方法进行模拟。

（2）夹紧。夹紧力和焊枪力将零件固定在夹具上，此时零件受力发生变形，零件与零件、零件与焊枪发生接触。

（3）焊接。焊枪通电，两薄板零件连接在一起，整体刚度发生变化。

（4）释放。撤去夹具，零件在自由状态下，仅进行简单支撑。板件在重力与内应力的作用下产生回弹，从而引起装配变形，这是装配件产生偏差的直接原因。

薄板零件装配误差分析的常用方法是直接蒙特卡罗法。通过对零件进行简单随机抽样、装配、测量得到装配后零件的变形，然后求出该变形零件某些统计量的无偏估计值，根据公差要求进行零件的公差分析。薄板零件装配误差分析的整个过程如图1所示。在公差分析的过程中，用有限元模拟代替实际的零件装配过程能够很好地缩短开发周期、降低开发成本。此时，抽样样本点的个数直接决定了有限元分析的次数。传统的简单随机抽样需要大量的抽样样本，会导致仿真消耗大量的计算时间，而拉丁超立方抽样能够以极少的抽样使计算收敛，而较少的计算使提高单次模拟的精度变得可能。

2 拉丁超立方抽样

拉丁超立方抽样是一种特殊的多维分层抽样方法，广泛应用于工程实际中。该方法中，抽样点的分布具有均匀性，能够很快速地达到收敛。其基本理论如下：

假设K维随机变量X的各个元素的概率分布函数为Fi（i＝1，2，…，K）。向量X的各元素相互独立，每个元素进行N次抽样，xjk为第k（k＝1，2，…，K）个元素的第j（j＝1，2，…，N）次抽样的值，定义N×K维矩阵P。P的每一列由数列｛1，2，…，N｝中各元素的随机排列组成。令随机变量ξjk服从区间［0，1］上的均匀分布，则抽样后得到的结果为［8］

式中，pjk为N×K维矩阵P的j行k列元素。

设存在函数h（x），定义函数h（x）的均值E（h（x））的无偏估计为

则简单随机抽样时的无偏估计的方差为

拉丁超立方的无偏估计的方差为

可以 证 明 （N－1）cov（h（x1n），h（x2n））/N概率逼近负值。因此拉丁超立方抽样比简单随机抽样更容易收敛。

3 接触分析

有限元仿真时，两零件在不考虑接触时很可能发生穿透或干涉，但是这些情况在现实的装配过程中是不可能发生的。实际焊接过程中，零件与零件、零件与焊枪、零件与夹具间必然会发生接触，而且接触会发生在从夹紧到回弹的整个装配过程中，因此在虚拟装配仿真中要充分考虑接触。接触分析在冲压成形和碰撞分析中非常关键。柔性薄板装配分析常用的影响系数法与直接蒙特卡罗法都忽略了零件间的接触作用。影响系数法只能构造出偏差源与变形的线性关系，因此不考虑接触的作用。直接蒙特卡罗法要进行大量的仿真计算，如果在单次仿真时考虑的模型过于复杂会使计算效率极低。

图2 零件接触单元

本文采用面面接触形式对接触过程进行模拟。零件与零件接触的接触类型为柔性体与柔性体接触，如图2所示。

考虑到零件的弹性，零件与焊枪（夹具）的接触属于柔性体与刚性体的接触，该接触类型的接触算法采用非线性罚函数法，在可能发生接触的地方添加接触单元。添加接触的模型在进行计算前先检查是否有穿透发生，如果发生，便在节点与表面间引入一个界面接触力（大小与穿透深度、接触刚度成正比），根据穿透深度构建惩罚函数。这个算法的意义相当于在两零件间放置了一系列的法向弹簧来限制穿透。图3、图4表示了零件是否考虑接触后得到的不同结果，可以看出，是否考虑接触的两种分析方法得出的结果差别比较大，而且考虑接触会更加接近真实模型。

图3 未考虑接触时零件发生穿透

图4 考虑接触时零件未发生穿透

本文通过采用细化网格的方法来减小单一节点穿透的可能，以提高接触算法的收敛速度。当零件在接触厚度内仅发生极小的干涉时，可以忽略这些极微小的接触来达到快速收敛的目的。

4 响应面模型

由于焊接误差分析过程中抽样得到的统计结果是概率收敛的，因此仅仅比较单次拉丁超立方抽样与简单随机抽样的结果是没有意义的，而需要进行多次比较。由于单次简单随机抽样的蒙特卡罗法就需要2000～5000次才能收敛，那么直接采用逐次仿真方法计算出的结果进行比较显然十分困难。为了实现两种方法的结果比较，本文构造了零件最终变形与偏差源的响应面模型，直接将抽样结果代入到高精度的响应面函数中计算出零件的变形。

响应面法采用实验设计理论对指定的设计点集合进行实验，得到目标函数和约束函数的响应面模型，从而预测非实验点的响应值［9］。

空间设计变量与响应Y的2阶多项式响应面模型如下［10］：

同理，可以构造出3阶、4阶甚至更高阶的多项式响应面近似模型，多项式系数A＝［a0a1…aL］T采用最小二乘法确定，即

式中，X为基函数矩阵；Y为P（P＞L）个实验点的响应。

响应面的精度常用响应面精度系数R2来验证，R2的表达式为

5 数值算例

将图5所示的两块薄板零件焊接在一起，零件存在偏差源x1和x2，薄板零件的厚度为1mm，弹性模量E＝20.7GPa，泊松比μ＝0.3。在零件可能发生接触的部分添加接触单元，摩擦因数为0.2。

图5 焊接零件及偏差源

通过样本采样后仿真出零件的变形来构造响应面函数，为减小试验次数，样本采用两变量五因素正交试验表构造，正交试验结果如表1所示。

表1 仿真样本 mm

构造出的焊接变形y的响应面函数如下：

通过试验设计点求得的响应面的误差决定系数R2＝98.25%。因此构造得到的响应面模型具有较高的精度，基本可以代替实际模型进行装配结果的比较，响应面模型如图6所示。对是否考虑接触这两种情况进行仿真得到的结果如图7所示。

图6 响应面图

图7 是否考虑接触时的两种仿真结果比较

假设两个偏差源x1、x2分别服从N1～（4，0.16）和N2～ （－1，0.16）的正态分布。为比较两种方法收敛的精度及效率，对零件使用拉丁超立方抽样法进行1000组误差分析，得到零件的1000个均值，同样计算出直接蒙特卡罗法的1000个均值，将两种结果落在指定区间的次数做对比，实现两种方法收敛精度的比较。通过这两组均值落在实际平均值附近一个较小区间的次数来实现这两种方法的比较，比较结果如表2所示。

表2 两种抽样方法结果比较

从表2可以看出，2000个实验样本的简单随机抽样计算出的零件变形收敛在区间［3.910，3.925］的频数是748，而80样本数的拉丁超立方收敛频数是954。 当收敛区间减小到［3.912，3.922］时，80样本数的拉丁超立方的收敛效率依然能保持较高的水平，达到872次，但是简单随机抽样的结果有很大的变化，变为593次，而且相比于1000样本数的简单随机抽样增加的幅度不大。这表明，当计算结果精度要求比较高时，简单随机抽样通过增加抽样次数提高精度的效果有限，但是拉丁超立方抽样效果较好。图8和图9分别是拉丁超立方抽样与简单随机抽样得到的均值的直方图，可以看出拉丁超立方抽样直方图均值分布比较集中，收敛效果明显好于简单随机抽样。

6 结语

图8 拉丁超立方抽样直方图

图9 简单随机抽样直方图

本文在薄板装配误差分析中，对直接蒙特卡罗法进行改进，采用拉丁超立方抽样代替以往的简单随机抽样。这种抽样方法能够大大减小抽样次数，提高收敛效率，进而能够对更加精确复杂的模型进行计算。本文在误差分析模型中考虑了零件与零件的接触以及零件与焊枪的接触，仿真结果表明，考虑接触的模型能够有效提高仿真的精度。为比较两种不同方法概率收敛的结果，构造了变形的响应平面，从表2可以看出，相同精度下拉丁超立方抽样的收敛次数远小于简单随机抽样，对于精度要求越高的分析，其收敛效果更好。

［1］Chase K W，Gao J，Magleby S P.General 2－D Tolerance Analysis of Mechanical Assemblies with Small Kinematic Adjustments［J］.Journal of Design and Manufacturing，1995，5（4）：263－274.

［2］Gao J，Chase K W，Magleby S P.General 3－D Tolerance Analysis of Mechanical Assemblies with Small Kinematic Adjustments［J］.IIE Transactions，1998，30（4）：367－377.

［3］Cai W，Hu S J，Yuan J X.A Variational Method of Robust Fixture Configuration Design for 3－D Work Pieces［J］.ASME Journal of Manufacturing Science and Engineering，1997，119（5）：593－601.

［4］Liu S C，Hu S J.Variation Simulation for Deformable Sheet Metal Assemblies Using Finite Element Methods［J］.ASME Journal of Manufacturing Science and Engineering，1997，119：369－374.

［5］Dahlström S.Variation Simulation of Sheet Metal Assemblies Using the Method of Influence Coeffcients with Contact Modeling［J］.Journal of Manu－facturing Science and Engineering，2007，129（3）：615－622.

［6］Dahlström S，Camelio J A.Fixture Design Methodology for Sheet Metal Assembly Using Computer Simulations［C］//Proceedings 2003ASME International Mechanical Engineering Congress and Exposition.Washington D C，2003：15－21.

［7］Camelio J A，Hu J S，Ceglarek D.Modeling Variation Propagation of Multi－Station Assembly Systems with Compliant Parts［C］//Proceedings of ASME 2001Design Engineering Technical Conference and Computers ＆Information in Engineering Conference.Pittsburgh，PA，2001：9－12.

［8］Stein M.Large Sample Properties of Simulations Using Latin Hypercube Sampling［J］.Technometrics，1987，29（2）：143－151.

［9］熊俊涛，乔志德，韩忠华.基于响应面法的跨声速机翼气动优化设计［J］.航空学报，2006，27（3）：399－402.

［10］孙光永，李光耀，张勇，等.基于鲁棒性的概率优化设计在薄壁构件耐撞性中的应用［J］.中国机械工程，2007，18（4）：479－483.