Int J Adv Manuf Technol (2011) 53:47–61 DOI 10.1007/s00170-010-2839-4 ORIGINAL ARTICLE Simulation of variational compliant assemblies with shape errors based on morphing mesh approach Pasquale Franciosa & Salvatore Gerbino & Stanislao Patalano Received: 22 October 2009 / Accepted: 8 July 2010 / Published online: 23 July 2010 # Springer-Verlag London Limited 2010 Abstract Variation analysis of assemblies is a strategic task in many industrial applications. Parts manufactured through plastic deformation processes exhibit appreciable shape deviations from the nominal geometry due mainly to spring-back phenomena. When these parts are assembled, initial shape deviations at part level highly influence the final assembly shape. This work focuses on the modeling and simulation of shape errors in order to perform variation analysis of compliant assemblies. The aim is to simulate variational shape of parts according to a small number of control points chosen on the part geometry through a morphing mesh procedure. These points are typically related to measurement or inspection points of manufactured parts. From the mesh model of parts, mesh nodes are moved by applying the morphing procedure. In particular, in order to assure control points belong to the “perturbed” shape, a linear-constrained approach is adopted. The so-morphed parts are used to accomplish the variational assembly analysis following the classical place, clamp, fasten, and release cycle. In order to achieve statistical results, Monte Carlo simulation is performed: a set of control points driving the perturbed parts is generated at each iteration; these parts are then assembled and results are stored. Numerical results are compared with ones coming from commercial software that uses a linear approach based on the sensitivity matrix. P. Franciosa : S. Patalano University of Naples Federico II, Naples, Italy P. Franciosa e-mail: pasquale.franciosa@unina.it S. Patalano e-mail: stanislao.patalano@unina.it S. Gerbino (*) University of Molise, School of Engineering, Termoli (CB), Italy e-mail: salvatore.gerbino@unimol.it Keywords Monte Carlo FEA . Morphing mesh . Shape errors . Geometric covariance . Variational assemblies 1 Introduction Predicting the final assembly shape of compliant non-ideal parts is a strategic topic for many industrial applications. Three sources of variation are typically identified in compliant assembly process: part shape variations, fixture variations, and fastening tool variations. Several numerical methods have been developed over the years to simulate and to take into account all these sources of variation. Most of these methods are based on linear variational models to provide a fast solution to the problem [1, 2]. Starting from deviations related to fastening or clamping points, a constant “Sensitivity Matrix,” linking deviations at part level to variations at assembly level, is built. Deviations data may be obtained from statistical analysis of the manufactured parts in the assembly line. The Sensitivity Matrix may be easily calculated by using any commercial FEA solver, running two consecutive FEA calculations and using the method of influence coefficients (MIC). Originally, the Sensitivity Matrix was defined in process involving single-station assemblies. Then, in Camelio et al. [3], this concept was extended also to multi-station assembly processes to consider the effects of variation accumulation when parts are moved from one station to another. Camelio et al. [4] proposed a more computationally efficient method than Liu et al. [1], also taking into account geometric covariance by means of the statistical principal component analysis method. However, linear approaches do not provide adequate results when large deformations occur, or part-to-part contact conditions have to be taken into account. Liao et al. [5] and Xie et al. [6] showed how contacts among parts being assembled highly influence final assembly shape. They proposed to use a non-linear FEA approach to solve 48 the contact problem. Their methods are more accurate than a linear one, but they are very time consuming especially if combined with Monte Carlo simulations. To overcome this lack, an interesting linear contact algorithm was proposed both in Dalhlström et al. [7] and Ungemach et al. [8], by combining a linear contact search and a contact equilibrium criterion. Their methods may be easily integrated with the MIC. In Gerbino et al. [9], a linear methodology, called statistical variation analysis and finite element analysis (SVA-FEA), which allows to statistically simulate both single- and multi-station assemblies, starting from the Sensitivity Matrix, was proposed. Also, in SVA-FEA, a linear contact algorithm was implemented by using multipoint constraints (MPC) elements of MSC.Nastran®. In Franciosa et al. [10], a comparison study between SVAFEA and a commercial CAT software, TAA® (module of CATIA® CAD system, by Dassault Systèmes), with reference to a multi-station assembly, was presented. The comparison showed high numerical correlation both for mean and standard deviations. Many methods and numerical algorithms have been developed over the years, but there are very few works on validation of the methods compared with real production inspection data. Hu et al. [11] presented a simulation method, which was verified for one instrument panel consisting of two parts, assembled with different weld sequences, with part and process variation. A validation of a more complex assembly was presented by Sellem et al. [12]. In Sellem et al. [13], authors concluded also that modeling aspects such as type and size of elements, incompatible and compatible meshes, and the use of contacts have minor impact on the simulation result. The present work focuses on the modeling of part shape deviations and provides a methodology to simulate variational compliant assemblies when these deviations occur. The aim is to generate variational shapes according to a small set of local deviations defined on nominal parts. Typically, local deviations, assigned at point level, cannot describe by itself the whole geometry. In the literature, this issue is analyzed by introducing the “geometric covariance” concept. The geometric covariance states the geometrical relation among the neighboring points on the same surface. Obviously, geometric covariance assures surface continuity and smoothness. Merkley [14] proposed to use random Bezier curves. Shape deviations were parameterized by constraining the displacement of the control points of Bezier curves. This method may be also extended to rectangular Bezier patches. However, for complex shapes, the parameterization of the patch becomes not a trivial task, so this method may be inadequate. Bihlmaier [15] extended Merkley’s work and proposed to model the variations of a surface as a finite summation of Int J Adv Manuf Technol (2011) 53:47–61 sinusoidal waves, each with a different amplitude and wavelength. Then, any surface profile is modeled as a summation of sinusoids, having different wavelengths and amplitudes, represented in the frequency domain using the Fourier transform. A more general method was proposed in Tonks [16]. To take account of surface variations, a hybrid method was used to model the surface covariance. Legendre polynomials were used to model the long wavelengths, and the frequency spectrum was used to model the shorter wavelengths. The hybrid method for geometric covariance was validated by means of experimental data. Huang et al. [17] suggested decomposing the shape deviation by using the direct cosine transform technique. The field of variation was divided into a set of independent defects. However, this method well works only for rectangular-based surface. An interesting approach was proposed in Samper et al. [18] by using a modal decomposition analysis. Firstly, the nominal CAD geometry is meshed and the orthogonal modal matrix is calculated. This matrix is then used to extract the principal modal shapes. To reach more accuracy, during the modal decomposition, a large amount of measurement data is needed. Moreover, a modal solver is required to calculate the modal matrix. A similar method was proposed also in Ungemach et al. [8], where the first eigenvector mode, resulting from a preliminary buckling analysis, was adopted to generate the initial variational geometry to be used in the assembly process simulation. Both the modal decomposition and the buckling eigenvector analyses have no physical significance. However, they give a good approximation of the geometrical covariance. Starting from modal decomposition, statistically modal analysis was proposed in Huang et al. [19] to statistically generate geometrical errors. This method decomposes the error field into few principal modes. In the present work, a morphing mesh-based approach is used to generate variational shapes, according to deviations occurring in a small set of points defined on the nominal geometry. The proposed method is based on a constrained deformation approach. Each deviation point (in the following named control point) determines a local deformation of the surface. A weighted function is generated in order to take into account the influence of all control points. Once parts are “perturbed,” according to the morphing approach, the place, clamp, fasten, and release (PCFR) cycle [20] is accomplished. This method does not require measurement data, and it can be adopted in the early steps of the design process in “what-if” scenarios. Anyway, measurement data can be used to validate an existing design. This paper is arranged as follows: Section 2 summarizes the proposed methodology; the morphing mesh procedure is illustrated in Section 3 and applied to a case study in Int J Adv Manuf Technol (2011) 53:47–61 Section 4; Section 5 shows the proposed methodology to simulate the PCFR cycle; Section 6 briefly describes the graphical user interface (GUI) developed to perform the variation analysis; a complete case study is analyzed in Section 7, in which also a comparison with a commercial software, able to analyze variational assemblies with a linear approach, is shown; finally, Section 8 draws the conclusions. 2 Methodology overview Figure 1 shows the general workflow proposed to perform the variational analysis of compliant assemblies with shape errors. The software architecture is based on MatLAB® environment [21], which controls the whole analysis, by running all FEA simulations in background mode through Comsol Multiphysics® [22] routines. Comsol Multiphysics® allows to solve PDE equations by means of a finite element approach. Expressions and functions can be easily handled as this environment allows to numerically evaluate integrals and derivatives defined on geometrical boundaries or subdomains. FE model is created from the nominal CAD geometry. Boundary conditions, contacts and identity pairs, and material properties are defined, accordingly. Identity pairs are used to simulate two parts that may move together as if they are perfectly glued. Fig. 1 General workflow of the methodology 49 Once FE model is defined, it is imported into MatLAB® where the classical PCFR cycle is simulated. Typically, during assembly processes, parts are positioned on fixtures. Then, they are clamped and fastened (for example welded). Finally, some fixtures are removed, and the assembled parts are released reaching, under the spring-back effect, their final configuration. In order to get statistical results, Monte Carlo approach is adopted (NMC is the total number of Monte Carlo simulations in Fig. 1). For each step of the simulation, the nominal geometry is perturbed according to deviations occurring at some control points (morphing mesh procedure). In this way, a deformed geometry is generated to be within a specified range (tolerance). Then, on this specific configuration (randomly generated) for each part, the PCRF cycle is performed at assembly level and the final assembly configuration is stored. Statistical results are calculated by considering all assembly configurations. Supposing to fasten parts with weld spots, how weld guns interact with (deformed) parts is described in the following, and a simplification to make faster the analysis is also proposed. 3 Morphing mesh procedure Morphing mesh is a well-known technique used in computer graphic applications as a powerful tool for shape modeling and designing [23, 24]. The user defines a set of control points, by giving a desired displacement and the relative radius of influence for each point. Control points directly influence the final shape of the “morphed” objects, and this shape can be fine tuned by adjusting the influence radius or the position of each control point. Variational parts, to be used in the PCFR assembly simulation, are generated here by means of a morphing mesh procedure (MMP). MMP works on a large amount of points describing the whole geometry. These points are directly derived from the FE model, used to perform computer simulation of the assembly process. Generally speaking, a FE model is made of nodes connected by elements. The number of nodes depends on the shape function of each element [26]. Let [V] be the matrix of nodes to be moved and [F] be the matrix of elements; [V] is an Nn ×3 rectangular matrix, where Nn is the total number of nodes; whereas [F] is an Ne ×3 rectangular matrix, where Ne is the total number of elements. Moreover, let [Pc] be the matrix of control points; [Pc] is an Ncp ×3 rectangular matrix, where Ncp is the total number of control points (see Fig. 2). Ri is the influence radius related to the ith control point. The region embedded within the influence radius is called here influence hull. When production data are known, control points can be directly derived from measurement data. Otherwise, during 50 Int J Adv Manuf Technol (2011) 53:47–61 In order to evaluate morphing matrix [M], Eq. 1 is specified with respect to all control points [Pc]. Then, one can write: 8 ΔPc1 ¼ fW ðPc1 Þg
½ M Š > > > > < ΔPc2 ¼ fW ðPc2 Þg
½ M Š ::: ! ½ΔPc Š ¼ ½W ðPc ފ
½ M Š    > > ΔPcNcp ¼ W PcNcp
½ M Š > > : where [ΔPc] is the Ncp ×3 matrix of displacements related to control points, whereas [W(Pc)] is a square matrix of Ncp order. Then, morphing matrix [M] can be calculated by solving the following linear system: Fig. 2 Definition of control points and relative influence hulls ½ M Š ¼ ½W ðPc ފ the early design stage, control points can be randomly generated by using a Monte Carlo simulation. The displacement, ΔVj, of any mesh node Vj, is calculated as a weighted mean of all control points [Pc], as stated in Eq. 1. 8 XNcp   > W Vj
Mi;x ΔVj;x ¼ > i¼1 i > < XNcp      W Vj
Mi;y ! ΔVj ¼ W Vj
½ M Š ΔVj;y ¼ i¼1 i > > XNcp   > : W Vj
Mi;z ΔVj;z ¼ i¼1 i ð1Þ where {W(Vj)} is the weighted (1×Ncp) row vector, and [M] is the morphing matrix (Ncp ×3). The ith weighted element depends on (I) position of the point Vj, (II) position of the control point Pci, and (III) its influence radius Ri, as stated in Eq. 2. The influence radius defines the 3D region within which any node is influenced by the related control point.     Wi Vj ¼ f Vj ; Pci ; Ri ; 8i ¼ 1; 2; ::; Ncp ð2Þ f is the basic function and may be assumed as a piecewise Bezier curve or a B-spline-based function as proposed in Raffin et al. [27]. However, in this work, in order to easily handle this function, a third degree polynomial function is adopted. Equation 2 can then be re-written as:  Vj   W i Vj ¼ f Ri  Pci  ¼ f ðdÞ ¼ 1 ð4Þ 1
½ΔPc Š It should be noted that matrix [W(Pc)] is singular only if two control points are coincident. When this happens, the system 5 has no unique solution and morphing matrix [M] is un-determined. To solve this issue, the linear system in Eq. 5 may be solved with a least-squares approach in which the pseudo-inverse of the matrix [W(Pc)] is evaluated [25]. Commercial routines are based on the singular value decomposition. In this paper, the MatLAB®’s backslash (\) operator is adopted. Once morphing matrix [M] is known, displacements of any mesh node Vj can be calculated by applying relationship 1. As it is, MMP allows to move nodes of the mesh taking into account the displacement of a set of control points. However, no information is available for points belonging to mesh elements. To overcome this limitation, an interpolation approach is used. The interpolation function allows to describe any point of the geometry starting from mesh nodes (whose displacements are calculated by applying Eq. 1). The postinterp function, embedded in MatLAB® via Comsol Multiphysics®, is adopted here. In order to get a more powerful control on the deformation of the nominal geometry, different influence hulls can be adopted. If the radius of influence is constant, 3
d 2 þ 2
d 3 ; 8i ¼ 1; 2; ::; Ncp ð3Þ Relation 3 states that function f is equal to 1 when Vj and Pci are coincident and tends toward zero for points Vj whose distance from Pci is greater than zero. To assure smoothing shape, function f has zero slopes at the two end points (Fig. 3). In Eq. 3, d is the a-dimensional distance from point Vj to control point Pci. ð5Þ Fig. 3 Definition of the basic function, f Int J Adv Manuf Technol (2011) 53:47–61 51 the influence hull becomes a sphere. Instead, an ellipsoid can be derived by defining three principal radii (defined along the principal directions of the ellipsoid). More general influence hulls can be defined as proposed in Raffin et al. [27]. However, only spheres and ellipsoids are considered here since they offer a more intuitively control. For a sphere hull, the a-dimensional distance, ds, is defined as: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  2  2ffi Vj;x Pci;x þ Vj;y Pci;y þ Vj;z Pci;z ds ¼ Rs ð6Þ whereas for an ellipsoid hull, it becomes: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u   2  2 u Vj;x Pci;x 2 Vj;y Pci;y Vj;z Pci;z de ¼ t þ þ R2e;1 R2e;2 R2e;3 Fig. 4 MMP case study: nominal geometry and control points ð7Þ where Rs is the radius of the sphere, while Re,1, Re,2, and Re,3 are the principal radii of the ellipsoid. Having assigned the above parameters, MMP may now generate variational parts. The user has to set the position of control points, their influence radius and relative displacements, and the interpolation function. 8 > < Um ¼ GðP; Pci ; Ri Þ ΔP ¼ GðP; Pci ; Ri Þ ! Vm ¼ GðP; Pci ; Ri Þ ð8Þ > : Wm ¼ GðP; Pci ; Ri Þ Then, with respect to the global coordinate frame (GCF), directly derived from the CAD model and here called (X0, Y0, Z0), any geometrical point P with (x, y, z) components is transformed by means of the MMP as in Eq. 8. G depends on both the basic and the interpolation function, whereas (Um, Vm, Wm) is the displacement of the point P with respect to GCF. Below, a case study clarifies how the proposed procedure works. 4 Case study: morphing mesh The proposed MMP is applied to the sheet metal part depicted in Fig. 4. FE model is made of about 3,000 nodes with tetrahedral mesh elements. 3D elements are required here, even for a model usually described with shell elements, as in the current release of Comsol Multiphysics® contact pairs cannot be defined on shell models. That is why all FEA simulations are performed with tetrahedral elements. Node displacements are interpolated by means of linear shape functions. Table 1 shows the statistical values assigned at control points, Pc, along the global Y axis (see Fig. 4). Two kinds of MMPs are performed. In the first one (MMPI), each control point influences the whole geometry (global influence). This means that each influence radius is equal to the maximum distance between any control point and the points of the part. A local influence, instead, is assigned into the second simulation (MMPII). In both cases, 1,000 simulations are performed using Monte Carlo method. Figures 5 and 6 show contour plots related to mean values and standard deviations, respectively (in the following, the Euclidean norm is adopted). Looking at Fig. 6, a twist effect due to deviations of control points Pc1 and Pc2, which are opposite each other, can be pointed out. Figure 7 depicts the influence hulls adopted in the MMPII simulation. Results in terms of mean values and standard deviations are, instead, shown in Figs. 8 and 9, respectively, for MMPII. Looking at Fig. 8, the reader can see that a local deformation of the geometry occurs according to the adopted influence hulls. 5 PCFR simulation A typical PCFR cycle starts with placing parts on fixture/s, followed by clamping parts on this/these fixture/s, then joining parts together, and finally releasing fixture/s and clamps. Table 1 Statistical values assigned at control points (Gaussian distribution assumed) Pc1 Pc2 Pc3 Pc4 Pc5 Mean (mm) Std dev (mm) 4.0 −3.0 −2.5 2.0 1.0 1.0 0.8 0.5 0.4 0.2 52 Int J Adv Manuf Technol (2011) 53:47–61 Fig. 5 MMPI—mean contour plot A critical aspect is related to the joining and the releasing phases. In fact, fastening tools force mating surfaces together at the joining area. Compliant parts tend to deform locally during these operations, and this changes contact conditions: Surfaces initially in contact may move apart, producing gaps or interferences. Moreover, depending on the joining method, multi-physics phenomena may occur in the contact zones, such as material plastic deformation, thermostructure interaction, thermo-electrical-structure interaction, and friction phenomena. Over the years, linear approaches, mainly based on the evaluation of the Sensitivity Matrix to do linear PCFR simulation, have been proposed. However, these methods do not consider the real behavior at the contact interface. In the present work, a non-linear methodology aiming to numerically simulate the whole PCFR cycle, taking into account contact area among parts being assembled, is proposed. To do it, a FEM approach is adopted. Figure 10 shows the general workflow used to perform the PCFR cycle. Generally speaking, three consecutive FEA runs are required to perform the simulation. The first FEA run is related to the positioning and clamping phase (phase I), considering initial variational parts, derived from MMP. Fastening phase (phase II) is achieved into the second FEA simulation, taking into account the deformed geometry derived from the first FEA run. Finally, the releasing phase (phase III) is resolved in the third FEA run, by applying the residual stresses (initial stresses) calculated during the previous phases. After the ith phase (i=I, II, III), the geometrical domain is updated according to the displacement field, (Ui, Vi, Wi). In this way, an updated geometry frame, here called deformed Fig. 6 MMPI—standard deviation contour plot Fig. 8 MMPII—mean contour plot Fig. 7 MMPII—influence hulls assigned at control points for the local influence of the geometry Int J Adv Manuf Technol (2011) 53:47–61 53 5.1 Positioning/clamping phase Fig. 9 MMPII—standard deviation contour plot geometry frame (DGF) is calculated. For the ith phase, DGFi is defined as: 8 > < Xi ¼ X0 þ ::: þ Ui DGFi : Yi ¼ Y0 þ ::: þ Vi > : Zi ¼ Z0 þ ::: þ Wi First, FEA run aims to simulate the location of parts onto fixtures. Generally speaking, fixtures are designed to determine location (mate constraint) of a part and to keep the effect of the location (contact constraint) [28]. In the context of the present paper, the positioning phase determines location while the clamping phase keeps such constraints. Changing the order in the positioning and clamping phases influences final results. However, in this paper, positioning and clamping on fixtures are assumed acting simultaneously. Then, the positioning/clamping phase will be simply named positioning phase briefly later on. Let DGF0 and DGFm be the nominal geometry frame and the morphing geometry frame, respectively. Looking at Fig. 11, the reader can see that parts, after applying MMP, may not fit boundaries corresponding to fixtures (defined on DGF0). First, FEA run is performed on the DGFm. Contact pairs are introduced to avoid part-to-part penetration. ð9Þ For example, the “fastening phase” is performed with respect to the “positioning geometry frame.” This frame is defined as: DGFpos 8 > < Xpos ¼ X0 þ Um þ Upos : Ypos ¼ Y0 þ Vm þ Vpos > : Zpos ¼ Z0 þ Wm þ Wpos ð10Þ In particular, DGF0 corresponds to the nominal CAD geometry. PCFR cycle calculation has to be iterated for each Monte Carlo step. Solving non-linear equations is required for each step. Convergence problems may additionally arise during the iterative solution of non-linear tasks. Furthermore, all these may become very difficult to automate. Due to the above mentioned reasons, an original methodology to solve for the contact problem between parts being joined will be showed in the following. In addition, the proposed methodology assumes the following hypotheses: & & & friction among parts is neglected; only spot weld joints are considered; and, local thermal effects due to welding process are not accounted. Fig. 10 Workflow of the non-linear methodology proposed for PCFR cycle 54 Int J Adv Manuf Technol (2011) 53:47–61 Fig. 11 Positioning modeling Furthermore, boundary constraints are added in order to fit the actual boundary on the fixture. To do this, boundary condition (BC) is defined as follows: 8 > < BCXm ¼ Um;bnd;fixture þ "x;pos BCYm ¼ Vm;bnd;fixture þ "y;pos > : BCZm ¼ Wm;bnd;fixture þ "z;pos ð11Þ where (Um,bdn,fixture, Vm,bdn,fixture, Wm,bdn,fixture) are the displacement fields calculated with MMP and related to boundaries to be fit on the fixtures. Variations on fixtures may be accounted, when necessary, as (εx,pos, εy,pos, εz,pos). 5.2 Fastening phase Fastening phase is performed with respect to DGFpos. When two parts are fastened with weld joints, firstly, weld guns force mating surfaces to close their gap, and then parts are joined together. Fastening phase is made in two consecutive steps: closing gap sub-phase (I) and joining sub-phase (II). In this phase mating surfaces do not penetrate each other; in fact, during positioning, penetrations among parts have been already avoided, by assigning contact pairs. Pairs (Ps,pos, Ps,0) and (Pd,pos, Pd,0) depend on the initial deformation field related to DGFpos. Point Ps,pos is calculated as intersection of the line, having direction Nweld and passing through point Pweld, with the mesh belonging to the boundary bnds,pos. The same procedure may be adopted for the opposite point. In this work, only triangular meshes are considered, but the proposed methodology may be extended to quadrilateral or mixed meshes. As described in Section 3, the mesh is characterized by means of nodes and elements. Here, the matrix of nodes is updated according to the DGFpos. Point Ps,pos is calculated by evaluating, on the soupdated mesh, the point Pint, which is the closest to point Pweld (see Fig. 13). To do it, for the ith triangular element, the intersection point Pint,i between the line (Nweld, Pweld) and the plane πi , whose element Fi belongs to, is calculated. Then, if point Pint,i belongs to the planar face Fi (inpolygon function available in MatLAB® is used to check this condition), it can be written:   Pint : mini Pint;i Pweld  ; 8i ¼ 1; ::; Ne ð13Þ Once gaps gs and gd have been evaluated, BCs may be assigned with respect to the local coordinate frame (Xweld, Yweld, Zweld). Then, it can be written: ( ð gs þ "fast Þif Δs ¼ Rweld BCs;Zweld ¼ free displacement; otherwise ( ð gd þ "fast Þif Δd ¼ Rweld ð14Þ BCd;Zweld ¼ free displacement; otherwise  Δi ¼ Ppos  Pi;def ; 8i ¼ s; d 5.2.1 Closing gap phase A weld spot is here characterized by its radius, Rweld, the nominal position of its center, Pweld, and the direction, Nweld, along which fastening tools act (see Fig. 12). Moreover, the source boundary (bnds) and the destination boundary (bndd) have to be introduced. In addition, a local coordinate frame (Xweld, Yweld, Zweld) is assigned having its local Zweld axis directed along Nweld vector. Looking at Fig. 12, the aim is to close gaps, defined as: ( gs ¼ Ps;pos gd ¼ Pd;pos Ps;0 Pd;0 : ð12Þ Fig. 12 Fastening modeling Int J Adv Manuf Technol (2011) 53:47–61 55 is the volume of weld spot sphere (see Annex A for additional mathematical details). mus ¼ 1 Ωs Ωs ¼ Ωs ΔΩs
us
dΩs mus if Δd <¼ Rweld free displacement; otherwise BCd ¼ ΔΩs ¼ Z ( 1 if Δs <¼ Rweld ð15Þ 0 otherwise Z ΔΩs
dΩs Ω  s  Δi ¼ Ppos Pi;0 ; 8i ¼ s; d Fig. 13 Line–mesh intersection procedure Relationship 15 states that any point of bndd, belonging to the sphere of radius Rweld and center Pd,0 (see Fig. 12), is constrained with a displacement equal to the mean displacement, μus, related to bnds. 5.3 Releasing phase where Ppos is any point defined into the DGFpos. The factor εfast is introduced in Eq. 14 to account also a variation, if any, in the nominal closure point of the weld guns. Readers have to consider that in the local coordinate frame, (Xweld, Yweld, Zweld), BCs are not assigned along local Xweld and Yweld directions: displacements are free along those directions. Finally, in the second FEA run, a constrained displacement, as in Eq. 4, is assigned at any point Ppos, belonging to the sphere of radius Rweld and center Pi,def and belonging to boundary bndi. Equations stated in Eq. 14 may be easily implemented in Comsol Multiphysics® by using boundary expressions and boundary constraints. Third, FEA run aims to simulate the final elastic springback, taking into account MPC conditions as defined in Section 5.2.2. In this phase, the geometry frame is updated with respect to the displacement field derived from fastening phase (DGFfas—note that DGFfas corresponds to DGFclose). Final displacements, (Urel, Vrel, Wrel), are calculated by applying to any part of the assembly the initial stress field calculated during previous steps (positioning and fastening phases). Then, it can be written: s rel;i ¼ s pos;i þ s fast;i ; 8i ¼ 1; :::; Np ð16Þ 5.2.2 Joining phase During the joining phase, bnds and bndd are constrained to move together (coupling constraint condition). Now, the geometry frame has to be updated taking into account the displacement field derived from the closing gap phase (DGFclose). An MPC condition is defined in order to couple bnds and bndd. Let us be the displacement field of bnds, available during joining phase. Coupling constraint conditions may be written as in relationship 15, where μus is the mean displacement related to points belonging to the weld spot sphere and belonging to the boundary bnds; Ωs where Np is the total number of assembled parts, while σpos and σfast are the stress fields related to the positioning and fastening phases, respectively. 6 Implementation and graphical user interface The proposed methodology was implemented into a GUI developed in MatLAB® environment, linked to Comsol Multiphysics® working in background mode (Fig. 14). FE model can be imported as structured array data directly from Comsol Multiphysics®. 56 Int J Adv Manuf Technol (2011) 53:47–61 Fig. 14 Graphical user interface MMP may be controlled from the VARIABILITY menu, where the user can set control points and their variability, in terms of mean value and standard deviation (statistical normal distribution is here assumed). Control points may be easily set by picking any point from the VIEW MODEL interface. JOINT menu allows the user to define fastening joints, setting radius, direction and nominal position. Boundary constraints and their related variability are set from the BOUNDARY menu. Final assembly shape, in terms of mean value and standard deviation displacements, can be shown directly in the Comsol Multiphysics® GUI, where the user can setup post-processing options. To understand what happens at point level through the whole assembly simulation, from INSPECTION menu, inspection points can be selected and their displacements may be exported into EXCEL® format. Comsol Multiphysics® solves contact problems by using the augmented Lagrangian method. This means that the software solves for equations in a segregated way (augmentation components are introduced for the contact pressure). A good alternative to augmented Lagrangian method is the penalty method. However, the latter is more stable but less accurate than the augmented Lagrangian method [29]. To implement penalty method within Comsol Multiphysics®, a contact normal penalty factor expression has to be defined as k/h, where k is a stiffness constant in the order of the stiffness of the material, while h is the mean dimension of the mesh. The proposed methodology requires that simulations of different assembly configurations (depending on the number of Monte Carlo simulations) are solved in automatic way (the user sets only initial solver parameters). In order to achieve this, a stable non-linear solution is required. Therefore, in the present work, the penalty Fig. 15 PCFR case study - nominal geometry, control points (Pci,j), spot welds (SPi), and inspection points (MSi) Int J Adv Manuf Technol (2011) 53:47–61 57 Table 2 Statistical values assigned at control points (Gaussian distribution assumed) Pc Part A Part B Pc1 Pc2 Pc3 Pc1 Pc2 Pc2 (Y) (Y) (Y) (X) (X) (Z) Mean (mm) Std Dev (mm) 0.3 2.0 1.0 2.0 −1.0 −1.0 0.1 0.3 0.6 0.5 0.1 0.1 method is adopted, despite that it is generally less accurate than augmented Lagrangian method. 7 Case study: PCFR cycle The proposed PCFR assembly methodology is applied to a two-part assembly (Young’s modulus E=210.000 N/mm2; Poisson’s ratio ν=0.30; uniform thickness T=2 mm), shown in Fig. 15. Overall dimensions of the assembly are about 800×600×500 mm. FE model is made of about 3,000 nodes with tetrahedral mesh elements (about 44,000 degrees of freedom). Two spot welds are assigned (SP1 and SP2 in Fig. 15), each with Fig. 16 a MMP phase. Generating variational parts. b Cumulative contour plot. Displacements at positioning phase. c Cumulative contour plot. Displacements at fastening phase. d Cumulative contour plot. Displacements at releasing phase 58 a radius of 20 mm. A fixed fixture is assigned for both parts (“fixed boundary” in Fig. 15). Contact pairs are defined at boundary interface between “part A” and “part B.” Moreover, linear shape functions are adopted to interpolate node displacements. Table 2 shows the statistical values assigned at control points. Five hundred Monte Carlo simulations were performed in less than 2 h on a notebook with Core Duo 1.83 GHz, 2 GB RAM, Win XP 32 bit. Figure 16 depicts mean contour plots for each phase of the PCFR cycle. The initial gap existing between part A and part B gets closed in the fastening phase. To understand what happens when contact pairs are not assigned, a new simulation is performed (run-time about 20 min). Results are shown into Fig. 17a, b with respect to inspection points (MS in Fig. 15), in terms of mean and standard deviations, respectively. In order to quantify the difference between the simulations with and without contact pairs, root sum square (RSS) index is calculated. Results are shown in Table 3. RSS index is quite low for MMP and positioning phases Fig. 17 a Mean deviations at inspection points. b Standard deviations at inspection points Int J Adv Manuf Technol (2011) 53:47–61 Table 3 Contact pair vs no-contact pair—RSS indices Phase MMP Positioning Fastening Releasing RSS (mean) RSS (std dev) 0.1460 0.1232 0.2445 0.2844 0.0619 0.0819 0.1091 0.1132 (when the number of Monte Carlo simulations increases, RSS index toward zero is expected). Moreover, RSS index becomes higher for fastening and releasing phases. This is due to the high influence that contact pairs assume both in the closing gap phase and in the final releasing phase. The above results are related to the non-linear approach described in the present paper. In the following, a comparison is made with the linear approach behind the commercial TAA module integrated into CATIA® CAD system, which implements the model described in Sellem et al. [2], mainly based on the sensitivity matrix. The case study depicted into Fig. 15 is Int J Adv Manuf Technol (2011) 53:47–61 analyzed in TAA environment, where contact points are also activated. Here, one should note that TAA allows creating shell mesh elements, whereas the proposed method, as described in Section 4, works only with 3D (tetrahedral) mesh elements for the limitations within Comsol Multiphysics on applying contacts on shell elements. Figure 18a, b shows mean values and standard deviations coming out from TAA simulations related to the same inspection points depicted into Fig. 15, respectively. Looking at mean displacements, readers can see that there is a good agreement among deviations related to points MS1 to MS6 and MS11 to MS15. The strong disagreement for points MS7 to MS10, instead, can be explained by the different approaches used to simulate geometric covariance. In TAA implementation, input deviations are assumed statistically independent. This means that no geometric covariance effect is accounted, but only “material covariance” [2] is simulated. The proposed approach, instead, allows to integrate both the geometric covariance, through the morphing mesh approach, and the material Fig. 18 a Proposed approach vs TAA method—mean values. b Proposed approach vs TAA method—standard deviations 59 covariance, by deforming parts during the PCFR cycle. Moreover, the control point Pc1,B strongly influences inspection points MS7 to MS10, through its influence hull. Looking at Fig. 18b, it is clear that, with respect to inspection points MS1 and MS15 (belonging to fixture boundaries—see also Fig. 15), standard deviations are zero for TAA results. This means that no cumulative effect can be simulated by using TAA approach. In other words, initial deviations that may occur at fixture boundaries do not affect final assembly deviations. In the described approach, instead, assembly deviations are cumulatively accounted and are considered at the final assembly level. 8 Final remarks In this paper, a methodology to perform statistical variational analysis of compliant parts, mainly based on a nonlinear FEA approach, is proposed. Statistical evaluation is 60 performed via Monte Carlo simulations to generate variational parts to be assembled and then released. The variational shape of parts being assembled is generated by using a morphing mesh-based approach borrowed by Computer Graphics. The user can manage and control the input shape by handling a small set of control points, assigning their position, the influence radius, and the relative displacement. In this way, no measurement data are directly necessary to perform a preliminary assembly simulation during the early design stage. Variational parts are, then, assembled following the classical PCFR cycle. The geometry domain is updated for each step of the PCFR cycle. The attention is focused on weld spot joints: how to simulate weld guns clamping parts is described. To improve the solution convergence, an efficient method, aiming to simulate fastening and releasing phases, is proposed. To do it, scalar and integral expressions, defined at boundary level, are defined and implemented. The proposed methodology was implemented into a friendly MatLAB®’s GUI, linked to Comsol Multiphysics® into background mode. The GUI drives user loading the initial FE model, assigning the input parameters, running Monte Carlo simulations, and finally, exporting output data. The proposed morphing mesh-based procedure is applied to a single sheet metal part, highlighting how, by managing control points and their influence hulls, different shape can be generated. Finally, a two-part assembly is simulated. Results point out that contacts may influence final assembly variations. Contacts allow to reach more realistic results but, on the other hand, a very time-consuming simulation is required in a contest of Monte Carlo simulations where hundreds of FEA calculations are usually run. Selecting the right balance between simulation time and goodness of numerical results is not a trivial task. Solving this issue requires more investigation. Furthermore, a comparative analysis is performed with the commercial TAA module, integrated into CATIA® CAD system, based on a linear approach. The analysis points out that TAA does not allow to manage “cumulative” effects when initial deviations are assigned at fixture boundaries, and no geometric covariance can be introduced. In the present approach, instead, deviations are accounted during the whole assembly process. Moreover, geometric covariance can be introduced through the morphing mesh-based model. More attention should be addressed also to different fastening joints, such as rivets and bolts, and to multistation assembly processes. ANNEX A: Integral definition on a boundary Let Ω be a domain. The aim is to define an integral on a sub-domain, Ωs. Supposing Ωs is a spherical sub-domain, Int J Adv Manuf Technol (2011) 53:47–61 and P a point belonging to Ω, the following condition is satisfied: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx x0 Þ2 þ ðy y0 Þ2 þ ðz z0 Þ2  R ! kP P0 k  R 8P 2 Ω ð17Þ where P0 and R are the center and the radius of the sphere, respectively. A step-wise function, ΔΩs , may be defined to account all points belonging to domain Ωs. Then, it has: ( 1 if kP P0 k  R ΔΩ s ¼ ð18Þ 0 otherwise Now, let l be a physical variable defined into the domain Ω, aiming to be integrated on the sub-domain Ωs, one can write: Z lint ¼ ΔΩs
l
dΩ ð19Þ Ω Taking in account Eq. 18, lint becomes (where Ωm =Ω−Ωs): lint ¼ Z Ωs ΔΩs
l
dΩs þ Z Ωm ΔΩs
l
dΩm ¼ Z Ωs ΔΩs
l
dΩs ð20Þ The main value of the physical variable l may be achieved by applying the definition of mean for a continuous function. Finally, one can write: Z 1 lm ¼ ΔΩs
l
dΩs ð21Þ Ωs Ωs References 1. Liu CS, Hu JS (1997) Variation simulation for deformable sheet metal assemblies using finite element methods. ASME J Manuf Sci Eng 119(3):368–374 2. Sellem E, Rivière A (1998) Tolerance Analysis of Deformable Assemblies, Proc. of DETC98, ASME Design Engineering Technical Conference, Atlanta, GA. (Paper Number: DETC98DAC5571) 3. Camelio J, Hu SJ, Ceglarek D (2003) Modeling variation propagation in multi-station assembly systems with compliant parts. ASME J Mech Des 125(4):673–681 4. Camelio J, Hu SJ, Marin SP (2004) Compliant assembly variation analysis using component geometric covariance. ASME J Manuf Sci Eng 126(2):355–360 5. Liao X, Wang GG (2007) Non-linear dimensional variation analysis for sheet metal assemblies by contact modeling. Finite Elem Anal Des 44:34–44 6. Xie K, Wells L, Camelio JA, Youn BD (2007) Variation propagation analysis on compliant assemblies considering contact interaction. ASME J Manuf Sci Eng 129(5):934–942 7. Dahlström S, Lindkvist L (2007) Variation simulation of sheet metal assemblies using the method of influence coefficients with contact modelling. ASME J Manuf Sci Eng 129(3):615–622 Int J Adv Manuf Technol (2011) 53:47–61 8. Ungemach G, Mantwill F (2009) Efficient consideration of contact in compliant assembly variation analysis. ASME J Manuf Sci Eng 131(1):1–9 9. Gerbino S, Patalano S, Franciosa P (2008) Statistical variation analysis of multi-station compliant assemblies based on sensitivity matrix, J. Comput Appl Technol 33(1):12–23 10. Franciosa P, Gerbino S, Patalano S (2010) Variation analysis of compliant assemblies: a comparative study of a multi-station assembly. Anales de Ingenieria Grafica 21:45–52 (in print) 11. Hu M, Lin Z, Lai X, Ni J (2001) Simulation and analysis of assembly processes considering compliant, non-ideal parts and tooling variations. J of Machine Tools and Manufacture 41:2233–2243 12. Sellem E, Riviére A, Hillerin CAD, Clement A (1999) Validation of the Tolerance Analysis of Compliant Assemblies, ASME Design Engineering Technical Conference, Las Vegas, Nevada, Sept. 12–15 13. Sellem E, Sellakh R, Riviére A (2001) Testing of Tolerance Analysis Module for Industrial Interest, 7th CIRP CAT Seminar, ENS de Cachan, France, April 24–25 14. Merkley K (1998) Tolerance analysis of complaint assemblies. Ph. D. Dissertation, Brigham Young University, Utah 15. Bihlmaier BF (1999) Tolerance analysis of flexible assemblies using finite element and spectral analysis. M.S. Thesis, Brigham Young University, Utah 16. Tonks M (2002) A robust geometric covariance method for flexible assembly tolerance analysis. M.S. Thesis, Brigham Young University, Utah 17. Huang W, Ceglarek D (2002) Decomposition mode-based off share form error by discrete cosine transformation with imple- 61 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. mentation to assembly and stamping systems with compliant shares. Annals of CIRP 51:21–26 Samper S, Formosa F (2007) Form defect tolerancing by natural modes analysis. ASME J Comput Inf Sci Eng 7(4):44–51 Huang W, Kong Z (2009) Simulation and integration of geometric and rigid body kinematics errors for assembly variation analysis. J of Manufacturing Systems 27:36–44 Chang M, Gossard DC (1997) Modeling the assembly of compliant, non-ideal parts. Comput-Aided Des 29(10):701–708 MatLab® (2009) User’s manual. University of Waterloo, Mathworks Comsol Multiphysics® (2009) User’s manual. COMSOL AB Borrel P, Rappoport A (1994) Simple constrained deformations for geometric modelling and interactive design. ACM Trans Graph 13:137–155 Raffin R, Neveu M, Jaar F (2000) Curvilinear displacement of free-form-based deformation. Vis Comput 16:38–46 Strang G (2009) An introduction to linear algebra, 4th edn. Wellesley Cambridge Press, Wellesley. ISBN 978-09802327–14 Zienkiewicz OC, Taylor RL, Zhu JZ (2005) The finite element method, its basis and fundamentals, 6th edn. Elsevier, Amsterdam Raffin R, Neveu M, Derdouri B (1998) Constrained Deformation for Geometric Modelling and Object Reconstruction, Proc. Winter School Computer Graphics, pp. 299–306 Whitney DE (2004) Mechanical assemblies: their design, manufacture, and role in product develpoment. Oxford University Press, New York Wriggers P (2002) Computational contact mechanics. Wiley, New York