A Second-Order Method For Assembly Tolerance Analysis

of Monte Carlo simulation increases with larger sample sizes.Obviously, the computational effort of large sample sizes can besignificant, but Monte Carlo simulation offers many advantagesbecause of its flexibility. Monte Carlo simulation allows anycomponent distribution to be specified and will calculate theresulting assembly distribution.1.Randomly Change AllModel VariablesFor each sample2.Evaluate Model Function3.Count RejectsFigure 3: Steps of Monte Carlo SimulationMonte Carlo simulation and the Linearized Method providedifferent capabilities. The Linearized Method can perform ananalysis and a tolerance allocation quickly, so it is suitable fordesign iteration. The Linearized Method is limited in that itcannot output non-normal distributions or handle non-normalcomponent distributions. Also, the Linearized Method will notbe accurate for highly nonlinear assemblies. Monte Carlosimulation allows non-normal input distributions and anonlinear analysis. However, Monte Carlo simulation iscomputationally expensive and does not accommodate rapiddesign iteration. For example, if a single input parameter ismodified, the entire Monte Carlo simulation must be re-run.Table 1: Comparison of Method FeaturesFeaturesSpeedTolerance allocationClosed-loop constraintsNonlinear approximationNon-normal inputdistributionsNon-normal outputdistributionsLinearizedMethod Monte CarloSimulation SOTAMethod Table 1 summarizes the features of the Linearized Method,Monte Carlo simulation and the Second-Order ToleranceAnalysis (SOTA) method proposed in this paper. The SOTAmethod attempts to combine the features of the LinearizedMethod and Monte Carlo simulation.The next section of this paper, Section 2, discusses researchrelated to the SOTA method. Section 3 presents the SOTAmethod. Section 4 compares the results of the SOTA methodwith the Linearized Method and Monte Carlo simulation for asample problem.2. RESEARCH REVIEW2.1 Linearized MethodThe Linearized Method, explained in Section 1.2, providesa quick way to perform nonlinear tolerance analysis for bothexplicit and implicit assembly dimensions of a vector-looptolerance model. Because of its speed, the Linearized Methodis ideal for design iteration and tolerance allocation. Multipleresearch studies have continued to refine the LinearizedMethod, making it more general and accurate.The Direct Linearization Method (DLM) [Marler 1988]prescribed a systematic approach to vector-loop modeltolerance analysis. DLM has enabled the Linearized Method tobe applied to a broad range of tolerance problems. Mostimportantly, DLM has allowed a general tolerance analysismethodology to be incorporated into a computer programsuitable for integration with a CAD system.More recently, the Global Coordinate Method [Gao 1993]for determining the partial derivatives of the loop equations wasdeveloped. This method simplified the calculations of thesederivatives.In the same paper, Gao benchmarked theLinearized Method against a comparable Monte Carlosimulation system. The benchmark results showed that theaccuracy of the Linearized Method corresponds to Monte Carlosimulation with a sample size of 30,000 for quality levels ofthree sigma.The Linearized Method has demonstrated its usefulness asa design tool. However, the method is inadequate for highlynonlinear tolerance problems and non-normal inputdistributions.2.2 Monte Carlo SimulationGenerally, Monte Carlo simulation is applied to an explicitfunction of random variables. However, the variables ofinterest in the equations of a vector-loop tolerance model areinherently implicit. McCATS, a Monte Carlo based toleranceanalysis method developed recently [Gao 1995], is able toadequately handle the implicit equations of a vector-looptolerance model.The McCATS system starts by generating random variatesfor the assembly variables. These random variates are sent toan assembly function that then solves the nonlinear system ofloop equations iteratively for the dependent assemblydimensions. The assembly dimensions are stored, new randomvariates are generated, and the assembly function is calledagain. This procedure is continued until the desired number ofassemblies has been simulated. Solving the loop equationsiteratively for each assembly simulation is critical to theaccuracy of the tolerance model.Including the capability for kinematic constraints in MonteCarlo simulation enabled Monte Carlo methods to be applied toa much broader range of design problems. However, therequired iterative assembly function does add more calculationsto an already computationally intense method.3Copyright 1999 by ASME

[]m2 (u i ) b 2j µ 2 (x j ) 2b j b jj µ 3 (x j ) b 2jj µ 4 (x j )2.3 Method of System MomentsThe Method of System Moments (MSM) [Cox 1979,Shapiro 1981] is a technique for estimating system output basedupon the relationship between input and output variables andinformation about the distribution of the inputs. MSM is alsoknown as nonlinear propagation of error and propagation ofmoments. MSM estimates the first four moments of a functionof random variables. The first four statistical moments are showin Figure 4.nj 1 [2bn 1 j 1 k j 1(9)]b b 2jk µ 2 (x j )µ 2 (xk )njj kkm3 (ui ) b 3j µ 3 (x j )n(10)j 1nn 1j 1j 1 k j 1m4 (u i ) b 4j µ 4 (x j ) n 6b b2 2j kµ 2 ( x j ) µ 2 (x k )(11)Where,First Moment:mean - measure of locationbj xThe termSecond Moment:variance - measure of spread-sFourth Moment:kurtosis - measure of peakednessFigure 4: First Four Statistical Momentsn 1 2nij 1 k j 1jjk(7)The number of terms in the expressions for the highermoments increases dramatically. For instance, the second, thirdand fourth moments require that the expected value be found forEquation 6 raised to the second, third and fourth power. Theexpressions for the higher moments are simplified if the originis shifted to the mean values. In terms of the new notation, thefirst four moments of the assembly dimensions areapproximated by the following four equations:m1 (u i ) b jj µ 2 (x j )nj 1(12)E (u i2 ) m2 (u i ) [m1 (ui )]E (u ) m3 (u i ) 3m2 (u i )m1 (u i ) 2[m1 (u i )]33iE (u ) m4 (u i ) 4m3 (ui )m1 (ui )4i 6m2 (u i )[m1 (ui )] 3[m1 (u i )](13)(14)(15)4(6) x j )(x k x k )1 n 2ui var(x j )2 j 1 x 2j 2ui x j xkµi (x j ) represents the ith distribution moment of the2The approximate system mean, the first distributionmoment, is calculated by taking the expected value of the aboveexpression, which gives:E (u i ) u i b jk 22 u x x (x1 2ui2 x 2jE (ui ) m1 (ui ) u iMSM is formulated by expanding the function of interest inTaylor series about its mean values. Retaining second-orderterms the expansion yields: u1 n u2u i u i i (x j x j ) 2i (x j x j )2 j 1 x jj 1 x jb jj jth component dimension.Equations 10 and 11 have been truncated significantly inorder to simplify the expressions. The complete third momentequation is lengthy, and the complete fourth moment equation isformidable. The complete equations for the third and fourthmoments may be found in [Cox 1979].After calculating Equations 8 through 11, the four momentsabout the mean may be found from: sThird Moment:skewness - measure of symmetryn ui x j(8)To estimate the four moments of an assembly distributionusing the full quadratic model requires the first eight momentsof the component dimension distributions and the partial22derivatives ui , ui and ui . x j x 2j x j xkIn a comparison of advanced tolerance analysis methods[Greenwood 1987], the Method of System Moments wasrecommended as the best method.3. THE SOTA METHODThe second-order tolerance analysis (SOTA) method isproposed as a general analysis method for vector-loop tolerancemodels. The SOTA method is comprised of a nonlinear systemsolver, finite difference approximations for the first and secondorder partial derivatives, the Method of System Moments(MSM), and a Generalized Lambda Distribution (GLD)empirical fit. The difference equations and nonlinear solver areused together to supply MSM with the required relationshipsbetween the component dimensions and the resultant assemblydimensions. MSM is then used to calculate the first fourmoments of the assembly dimensions. Finally, GLD is used tofit the calculated moments and approximate the distribution of4Copyright 1999 by ASME

the assembly dimensions. The SOTA method process is shownin Figure 5.Change One (or Two)Model Variable(s)1.For each variable2.Evaluate Model Function3.Calculate Sensitivity4.Calculate MomentsTwo of the function evaluations that appear in this threepoint difference formula also appear in the central differenceformula, Equation 16. Therefore, if there are n componentdimensions, 2n function evaluations can be avoided if the samefunction evaluations are used for both difference Equations 16and 17. If this is done, the quadratic partial derivatives willrequire 2n function evaluations plus one evaluation at thenominal and, without any further evaluations, the linear partialscan also be obtained.The approximation of the cross-derivatives is morecomplicated. These partial derivatives are found by using thecentral difference of a central difference.A B 5.Fit Distribution6.Calculate Rejects3.1 Difference FormulasIn order to approximate the three sets of partial22derivatives, ui , ui and ui , three separate difference x 2j x j xkformulas are required. The linear partial derivatives areapproximated by a central difference formula: ui ui (x j x j , ui ) ui (x j x j , ui ) 2 x j x j(16)For this notation, ui (x j x j , ui ) represents a functionevaluation for the implicit assembly dimensions ui , where thex j , are at their nominal value exceptfor the jth dimension, which is perturbed by a value x j . So,component dimensions,Equation 16 indicates two function evaluations for eachcomponent dimension x j . Note that each function evaluationrequires an iterative solution of a system of nonlinear equations.A three-point difference formula is used for theapproximation of the quadratic partial derivatives [Burden1993]. 2 ui ui (x j x j , ui ) 2ui (x j , ui ) u i (x j x j , ui ) x 2j x 2jui (x j x j , xk xk , ui ) ui (x j x j , xk xk , ui ) 2ui(A B ) 2 xk x j xkFigure 5: Steps of the SOTA Method x jui (x j x j , xk xk , ui ) ui (x j x j , xk xk , ui )2 x j(17)2 x j(18)For n component dimensions there are (n2-n)/2 uniquecross-derivatives. Four new function evaluations must beperformed for each derivative. Therefore, 2n2-2n evaluationsare required to obtain the cross-derivatives. Together with the2n 1 function evaluations for the linear and quadratic partialderivatives, the total becomes 2n2 1 function evaluations.Thus, with 2n2 1 function evaluations, the required partialderivatives for the SOTA method are obtained.3.2 Distribution FitThe Generalized Lambda Distribution was chosen as thebest empirical model for fitting the statistical moments based onease of implementation. The single form of the GLD make themethod of matching of moments easily applied to the momentscalculated by MSM, whereas, the Johnson and Pearson systemsrequire multiple distribution forms to cover a full range ofmoments. In addition, a GLD table, indexed by skewness andkurtosis values, was readily available for use in a computerprogram because of earlier research [Gao 1995]. The GLD'srange of coverage is smaller than the Johnson and Pearsonsystems, however, it does cover most practical distributionshapes likely to be encountered in mechanical assemblies.4. EXAMPLEThe following One-Way Clutch example problemillustrates the performance of the SOTA method compared toMonte Carlo Simulation and the Linearized Method. Theexample problem was analyzed using the SOTA method, theLinearized Method and Monte Carlo simulation at four differentsample sizes. The four sample sizes were 30,000 samples,100,000 samples, 106 samples and 109 samples.5Copyright 1999 by ASME

A one-way clutch transmits torque in a single direction.The clutch assembly consists of the following components: ahub, an outer ring, four rollers, and four springs. When the hubrotates in a counter-clockwise direction, the roller wedgesbetween the hub and the ring, locking these two parts together.When the hub turns in a clockwise direction, the spring iscompressed by the roller, the roller slips, and the hub is allowedto rotate freely. The one-way clutch assembly and the singlevector loop used to model this assembly are shown in Figure 6.φ1ySpringRingDHubBHubARollerC RollerRingETable 5 contains the predicted parts-per-million (PPM)assemblies that fall outside of the specification limits of thepressure angle. All the Total Rejects results were within 1000PPM of the Monte Carlo simulation of 109 samples.Table 4: Statistical Moments ResultsAnalysisMC 1e97.014953StandardDeviation0.219668MC 6-0.094773.027695MC 100k7.0154530.220172-0.101683.021511MC e 5: PPM Rejects Resultsφ2Figure 6: One-Way Clutch AssemblyThe function of the one-way clutch mechanism is governedby the pressure angle φ1. There are three manufactureddimensions that control the pressure angle. The mean values,standard deviations and distribution types of these threedimensions are shown in Table 2. The vectors representing theroller radius, vectors C and D, were treated as the samevariation source. As the pressure angle, φ1, is critical to thefunction of the clutch, it was given specification limits as shownin Table 3.MeanAnalysisLower RejectsUpper RejectsMC 1e944062166Total Rejects6572MC 1e6446722066673MC 100k458020806660MC n order to compare the results of the six analyses, arelative error measure was calculated for the estimates of thestatistical moments and the estimate of PPM rejects. MonteCarlo Simulation with 109 samples was assumed to be the mostaccurate analysis and was, therefore, used as the baseline for therelative error comparison.Table 2: Input VariablesNameAC, DEMean27.645 mm11.430 mm50.800 mmStandard Deviation0.01666 mm0.00333 mm0.00416 mmDistributionNormalNormalNormalPercent Error of the Mean0.050%0.045%0.040%Table 3: Assembly SpecificationNameφ1Nominal7.0184 Lower Limit6.4184 0.035%Upper Limit7.6184 0.030%0.025%0.020%4.1 Analysis ResultsThe One-Way Clutch analysis results are displayed in Table4 and Table 5. Table 4 shows the calculated values for the firstfour statistical moments of the pressure angle. The One-WayClutch assembly was a good test problem since the pressureangle exhibits nonlinear behavior. The results show the pressureangle to be negatively skewed and slightly more peaked than aNormal distribution. Because all the input distributions weresymmetric, this skewness indicates that the pressure angle is aninherently nonlinear function. Of course, the skewness valuecalculated by the linear analysis was zero since the linearanalysis cannot estimate this non-linearity.0.015%0.010%0.005%0.000%MC 1e6MC 100kMC 30kSOTALinearFigure 6: Error of the MeanFigure 6 compares the error of the pressure angle meanwith respect to the Monte Carlo 109 analysis. All five analysesestimated the mean very accurately to within 0.05% error. TheSOTA method was the most accurate with only 0.0002% error.6Copyright 1999 by ASME

Figure 7 displays the error of the standard deviation. Thestandard deviation values were also very accurate with the errorranging from 0.10% for the Monte Carlo 106 to 0.40% for theMonte Carlo 30k. Both Figure 6 and Figure 7 clearly illustratehow Monte Carlo Simulation should increase in accuracy as thesample size increases.Distribution was fit to the four moments in all six analyses. Thereject results are a composite result of all four statisticalmoment estimates and, therefore, provide a good overallmeasure of accuracy for an analysis.Percent Error of the Kurtosis2.00%Percent Error of the Standard C 1e60.05%0.00%MC 1e6MC 100kMC 30kSOTAThe truncation of nonlinear terms of the linear analysis isevident in the skewness results. With symmetric inputdistributions, a linear analysis will always predict a skewnessvalue of zero. Figure 8 shows the absolute error of theskewness. With the exception of the linear result all skewnessvalues are relatively accurate.MC 30kSOTALinearFigure 9: Error of the KurtosisLinearFigure 7: Error of the Standard DeviationMC 100kFigure 10 shows the error of the upper, lower and totalPPM rejects relative to the Monte Carlo 109 analysis. The linearanalysis predicted symmetric rejects: 3109 ppm for the lowerlimit and 3109 ppm for the upper limit. The linearappoximation of the nonlinear pressure angle function resultedin an underestimate for the lower rejects and an overestimate forthe upper rejects.Error of the PPM RejectsError of the SkewnessLower Rejects0.10Upper RejectsTotal 04-2000.03-4000.02-600-8000.01-10000.00MC 1e6MC 100kMC 30kSOTA-1200Linear-1400MC 1e6Figure 8: Error of the SkewnessThe percent error of the kurtosis values is shown in Figure9. With the exception of the Monte Carlo 30k result, all thekurtosis values had errors under 1%.The process of calculating rejects for an assemblyspecification involves the four statistical moments and fitting adistribution to these moments. The Generalized LambdaMC 100kMC 30kSOTALinearFigure 10: Error of the RejectsThe second-order approximation of the SOTA methoddramatically improved the estimate of rejects over the linearapproximation. While the SOTA method still slightlyunderestimated the lower rejects and slightly overestimated theupper rejects, the total rejects estimate was within 54 ppm of7Copyright 1999 by ASME

the Monte Carlo 109 analysis, roughly equivalent to the MonteCarlo 100k results.4.2 Computational EffortFor the three methods, Monte Carlo simulation, the SOTAmethod and the Linearized Method, a relative measure of effortis easily formulated. The common operation of these threetolerance analysis methods is that each must perform a linearsolution of the loop equations. For example, the SOTA methodand Monte Carlo simulation require a linear solution of the loopequations for each iteration of Newton's method. Of course theLinearized Method only requires a single linear solution. So, ifthe Linearized Method is given an effort value of 1, the relativeeffort of Monte Carlo simulation and the SOTA method may beevaluated by the following expressions:MC Effort (sample size) x (average Newton iterations)SOTA Effort (2n2 1) x (average Newton iterations)For the SOTA Effort expression the variable n is thenumber of component dimensions.It would be expected that the number of iterations ofNewton’s method be greater for Monte Carlo simulation thanfor the SOTA method. For each nonlinear solution, MonteCarlo simulation changes the nominal value of all thecomponent dimensions, whereas the SOTA method onlychanges one or two component dimensions for each solution.Furthermore, the step size used by the SOTA method willgenerally be very small compared to the variations required byMonte Carlo. The average number of Newton iterations alongwith the effort metrics for the six analyses is shown in Table 6.Table 6: Relative Computational EffortAnalysisMC 1e9MC 1e6MC 100kMC 1Effort3,400,000,0003,400,000340,000102,3004115. CONCLUSIONSThe SOTA method is a general, nonlinear toleranceanalysis method for vector loop tolerance models. The SOTAmethod provides the benefits of speed, tolerance allocation,closed-loop constraints, a nonlinear approximation and thecapability for non-normal input and output distributions.For the One-Way Clutch example problem the SOTAmethod shows a dramatic improvement in accuracy over thelinear approximation for the estimates of the four statisticalmoments of the pressure angle. The estimate of total rejectsresult for the SOTA method was comparable to the Monte Carloresult using 106 samples. This accuracy level is significant sincethe computational effort of the SOTA method was five orders ofmagnitude less than the Monte Carlo simulation with 106samples.Seven additional example problems were analyzed [Glancy1994] and demonstrated similar results.REFERENCES[Burden 1993] Burden, Richard L. and J. Douglas Faires,Numerical Analysis Fifth Edition, PWS Publishing Co.,1993[Chase 1995] Chase, Kenneth W., J. Gao and S. P. Magelby.“General 2-D Tolerance Analysis of MechanicalAssemblies with Small Kinematic Adjustments.” J. ofDesign and Manufacturing, 5 (1995): 263—274[Cox 1979] Cox, N. D., “Tolerance Analysis by Computer”, J.of Quality Technology, Vol. 11, No. 2, April 1979[Cvetko 1997] Cvetko, Robert. “Characterization of AssemblyVariation Analysis Methods”, MS Thesis, Brigham YoungUniversity, Dec. 1997[Dudewicz 1974] Dudewicz, E. J., J. S. Ramberg and P. R.Tadikamalla, “A Distribution for Data Fitting andSimulation”, An. Tech. Conf. Trans. A.S.Q.C., 28 (1974):407—418[Gao 1995] Gao, J., K. W. Chase and S. P. Magleby,“Comparison of Assembly Tolerance Analysis by theDirect Linearization Method and Modified Monte CarloSimulation Methods”, Proc. of the ASME DesignEngineering Tech. Conf., 1995, 353—360[Glancy 1994] Glancy, Charles G., “A Second-Order Methodfor Assembly Tolerance Analysis”, MS Thesis, BrighamYoung University, Dec. 1994[Greenwood 1987] Greenwood, W. H., “A New ToleranceAnalysis Method for Designers and Manufacturers”,Dissertation, Brigham Young University, 1987[Marler 1988] Marler, Jaren D., “Nonlinear Tolerance AnalysisUsing the Direct Linearization Method”, MS Thesis,Brigham Young University, 1988[Shapiro 1981] Shapiro, Samuel S., Alan J. Gross, StatisticalModeling Techniques, Marcel Dekker, Inc., New York,1981[Ramberg 1979] Ramberg, J. S., P. R. Tadikamalla, E. J.Dudewicz and E. F. Mykytha, “A probability distributionand its uses in fitting data”, Technometrics, 21 (1979):201—2148Copyright 1999 by ASME

a quick way to perform nonlinear tolerance analysis for both explicit and implicit assembly dimensions of a vector-loop tolerance model. Because of its speed, the Linearized Method is ideal for design iteration and tolerance allocation. Multiple research studies have continued to refine the Linearized Method, making it more general and accurate.