Global Sensitivity Analysis of Chosen Harmonic Drive Parameters Affecting Its Lost Motion

The aim of the presented study was to perform a global sensitivity analysis of various design parameters affecting the lost motion of the harmonic drive. A detailed virtual model of a harmonic drive was developed, including the wave generator, the flexible ball bearing, the flexible spline and the circular spline. Finite element analyses were performed to observe which parameter from the harmonic drive geometry parameter group affects the lost motion value most. The analyses were carried out using 4% of the rated harmonic drive output torque by the locked wave generator and fixed circular spline according the requirements for the high accuracy harmonic drive units. The described approach was applied to two harmonic drive units with the same ratio, but various dimensions and rated power were used to generalize and interpret the global sensitivity analysis results properly. The most important variable was for both harmonic drives the offset from the nominal tooth shape.


Introduction
Harmonic drive gearing, which is also called strain wave gearing, is a very progressive transmission with broad spectrum of use in applications with a high demand for precise positioning, e.g., aerospace, robot arms, test devices, medical equipment, etc. Harmonic drive (HD) consists of three main components: the wave generator (WG), the flexible spline (FS) and the circular spline (CS). One of the main parameters, which influences the selection and application of the HD, is the lost motion value, also called as torsional backlash.
Backlash is the amount by which the width of a gear's tooth space exceeds the thickness of an engaging tooth measured at the pitch circle of the gears, whereas the transmission lost motion is a measurement of the angle of movement in both directions before movement at the input occurs while applying a load at the output. Backlash, also known as slop, lash, free-play, or just play, is an angular quantity caused by the circular geometry of the gear. This is known as clearance backlash, and it is required to account for manufacturing defects, providing space for lubrication, and allowing for component thermal expansion. The HDs are usually designed with zero backlash at the gear mesh.
The producers of high accuracy HD units guarantee various values of the maximum torsion ϕ during their operation life measuring the angular lost motion of the FS by low output torque values. The test is usually performed by applying the 4% of rated HD output torque T N to the FS and measuring its angular deformation ϕ 1 , ϕ 1 (ϕ 2 , ϕ 2 , respectively) in both directions, with the WG locked and the CS fixed, obtaining the hysteresis curve shown in Figure 1. φ-angular deformation by applying the +4 φ′-angular deformation by applying the − The lost motion of the HDs with various ch by simulating the lost motion test procedure usi processing the acquired data by appropriate ma An example of a dynamic model of the HD presented by Folega [1]. He discussed certain a monic drive and proposed a new original dynam transmission system. His model takes into acco also the damping.
The approach by particular HD components The main goal of this article is to present a study that is focused on the global sensitivity analysis of the variable HD component dimensions influencing the lost motion value Φ of described HD units defined as where: ϕ-angular deformation by applying the +4% of T N to the FS; ϕ -angular deformation by applying the −4% of T N to the FS. The lost motion of the HDs with various chosen component parameters is obtained by simulating the lost motion test procedure using the finite elements analysis (FEA) and processing the acquired data by appropriate mathematical methods.
An example of a dynamic model of the HD in a toothed transmission system was presented by Folega [1]. He discussed certain aspects of dynamic modeling of the harmonic drive and proposed a new original dynamic model of a harmonic drive for a power transmission system. His model takes into account the nonlinear stiffness changes and also the damping.
The approach by particular HD components modeling can vary according their complexity. The early original wave generator outer shape was an elliptoid described by Musser [2]. Further research led to the configuration curve defined by a four-term Fourier expansion equation defined by Ishikawa [3]. Gravagno et al. [4] discussed and quantitatively evaluated the influence of the wave generator shape on the pure kinematic error of HD. The Résal curve, an ellipse and two conjugate arcs curves yield the lowest kinematic error. The cam curves yield slightly higher kinematic error as a compromise including also a caterpillar effect at the poles of the WG.
There are also various teeth splines used by the common produced harmonic drives. The harmonic drive company uses S-shaped teeth with the geometry based on the patent of Ishikawa [5]. This tooth profile allows for up to 30% of the total number of teeth to be engaged at the same time. In addition, when compared to an involute tooth, the big tooth root radius boosts dental strength. This technological innovation produces high torque, high torsional rigidity, extended life, and smooth rotating as also described in [6]. The teeth shape optimization is presented, e.g., by Kayabashi [7], Kiyosawa et al. [8] and also in further patents of Ishikawa. The parametric design of double-circular-arc tooth profile and its influence on the functional backlash of harmonic drive was already described by Chen et al. [9]. They used the envelope theory to design a double circular arc tooth profile of rigid gear to achieve multi-section conjugation. The presented method is more consistent with the actual working condition of the harmonic gear compared with the method of only designing the tooth profile of the middle section conjugate. The geometry of the FS influences the torsional stiffness of the HD significantly.
The problem of the torsional stiffness affecting also the harmonic drive backlash value was analyzed in various references. Yang et al. [10] observed the backlash of the harmonic drive based on a parametric solid finite element model. They calculated the backlashes with the minimum circumferential distance between engaged involute tooth profiles along tooth height direction using the four-rollers wave generator. They observed that the validity of the backlash computation based on the FEM is higher by the contrastive analysis. Rhéaume et al. [11] presented a finite-element model to reproduce the behavior of the torsional stiffness of a harmonic drive, which allowed the evaluation of the effects of various geometrical parameters on the torsional stiffness.
Wang et al. [12] presented the measurement and analysis of backlash on harmonic drives, where the backlash of harmonic drive transmission is split into two primary components: clearance-induced lost motion and elastic deformation-induced lost motion. The tooth clearance between circular spline and flexspline, as well as the clearance from Oldham coupling in the wave generator, are the main sources of clearance, whereas the elastic deformation is primarily attributable to flexspline under torque load. Zou et al. [13] described the measurement and modeling of kinematic error and clearance in harmonic drive. They established an experimental setup to obtain the kinematic error and clearance in harmonic drive. The position-dependent kinematic error was modeled by the Fourier expression method.
The general global sensitivity analysis (GSA) approach is described in detail, e.g., by Saltelli et al. [14]. Many references deal with GSA of various gearbox components. Woll et al. [15] carried out a sensitivity analysis on the drivetrain of an offshore winch with active heave compensation. By adjusting the parameters within usual ranges and assessing the quantitative effect on the lifetime, the effects of the parameters' gear center distance and ambient temperature were explored. Li et al. [16] performed a sensitivity analysis on a synchronization mechanism for manual transmission gearbox using parameterized virtual models of the typical synchronizers, and the paper discusses the influence of the main parameters of synchronizer on the gear shift performance. A sensitivity analysis for the proper tolerance choice of gear profiles was carried out by Rajabalinejad [17]. He explained how to use a sensitivity chart to identify the acceptable and ideal tolerance fields for profile modification as a tool for assessing the influence of profile alterations on the load factor of gears with stochastic indexing faults and finding optimal solutions. Ostapski [18] analyzed the harmonic drive wave generator-flexspline system in relation to selected structural parameters and manufacturing deviations. Material research, as well as simulations of the stress state in the bearing versus manufacturing deviations, and fits between the bearing and the generator cam were used to describe the problem of failure of the elastic bearing supporting the generator in a harmonic drive.

Materials and Methods
The lost motion (torsional backlash) is a key parameter defining the accuracy of the HD. When subjected to the rated torque, HD display characteristics shown in the hysteresis curve. When a torque is applied to the output shaft of the gear with the input shaft locked, the torque-torsion relationship can be measured at the output.
The lost motion is the term used to characterize the HD torsional stiffness in the low torque region. Lost motion is evaluated by a test applying the ±4% of rated HD output torque T N to the FS measuring its angular deformation in both directions with the WG locked and the CS fixed to obtain the torsional windup of the gear.
The standard kinematic gear ratio i of the HD is defined by the WG acting as the input shaft with the angular velocity ω 1 (revolutions per minute n 1 ), the FS with the teeth Materials 2021, 14, 5057 4 of 23 number z FS (−) acting as the output shaft with the angular velocity ω 2 (revolutions per minute n 2 ) and the CS with the teeth number z CS (−) fixed according the equation: The nominal parameters of harmonic drive units, which are analyzed in the presented research and denoted as HD 42 and HD 120 according the bearing outer diameter (flexspline inner diameter), are listed in Table 1. The main goal of the presented research was to identify the most important parameters picked up from a parametric dimension set defining the HD geometry with the maximum influence to the lost motion value, and subsequently to the overall torsional backlash value by various motion states of the HD unit. The goal is to obtain the lost motion value Φ of maximum 1 arcmin alternating the identified most important HD unit parameters by next development including HD key components design changes.
The complex parametric modeling of the virtual tested HD units was carried out with detailed geometry of the WG, the ball bearing with flexible rings and the FS and CS with teeth splines. The set of the parameters influencing the backlash value was picked up. The research process was then executed in the following steps: 1.
the parametric model of the HD finite element analyses was built and the analyses of the prescribed motion state were performed to determine the lost motion values for the generated set of chosen variables; 2.
the correlation of the input parameters (flexible ball bearing inner and outer ring diameter, the radial bearing clearance and the offset of the FS and CS nominal teeth spline) and the output parameter (lost motion) was observed and the correlation coefficients were calculated; 3.
the global sensitivity analyses (GSA) of the chosen variable sets were performed based on the obtained dependence of the input and output variables calculating the GSA indices according the generated quasi-random variable sets.
The calculation routine was performed for both HD units with various nominal power value to support the generalization of the calculation results, which were interpreted in the conclusion.

Parametric Model Building
The parametric model of the HD was created in the Creo parametric software (PTC, Inc., Boston, MA, USA, 3.0) according the equations describing the variable component parameters generated in the Matlab software (Mathworks, Inc., Natick, MA, USA, 2018R2) as described by Majchrak et al. [19]. The prepared geometrical models were imported to the analysis software Ansys Workbench (Ansys, Inc., Canonsburg, PA, USA, 17.2) and processed there. There was a calculation routine built, which provided the automatic HD parametric model change according the needs of the FEA calculations. The complete parametric model of the examined HD 42 unit is presented in Figure 2. The tested HD wave generator (WG) actual outer shape radius R is defined by the cam major axis A, cam minor axis B and the exponent n for the angle 0 ≤ φ ≤ π/2 by the equation: The variables defined in the WG outer shape equation influence the circumference of the WG and hence the inner ring inner diameter of the assembled ball bearing. The model of the ball bearing with flexible rings located on the WG outer shape is created by the method described by Steininger et al. [20]. The bearings' measures are based on standard products of the PBF Krasnik factory (type 113-1145TN for the HD 42 and type 114-876TN for the HD 120) with characteristic dimensions shown in Figure 3 and the nominal dimensions listed in Table 2 [21].  The tested HD wave generator (WG) actual outer shape radius R is defined by the cam major axis A, cam minor axis B and the exponent n for the angle 0 ≤ ϕ ≤ π/2 by the equation: The variables defined in the WG outer shape equation influence the circumference of the WG and hence the inner ring inner diameter of the assembled ball bearing. The model of the ball bearing with flexible rings located on the WG outer shape is created by the method described by Steininger et al. [20]. The bearings' measures are based on standard products of the PBF Krasnik factory (type 113-1145TN for the HD 42 and type 114-876TN for the HD 120) with characteristic dimensions shown in Figure 3 and the nominal dimensions listed in Table 2 [21]. The tested HD wave generator (WG) actual outer s cam major axis A, cam minor axis B and the exponent n equation: The variables defined in the WG outer shape equatio the WG and hence the inner ring inner diameter of the as of the ball bearing with flexible rings located on the W method described by Steininger et al. [20]. The bearings' products of the PBF Krasnik factory (type 113-1145TN fo for the HD 120) with characteristic dimensions shown in sions listed in Table 2 [21].   The nominal radial internal clearance (ISO clearance class) and the ovalization of the inner ring influencing the range of diameters used by the analyses, is defined by the equation: where: d max -maximum permissible value of the inner bearing diameter d; d min -minimum permissible value of the inner bearing diameter d.
The standard bearings of the selected producer do not meet the required radial internal clearance values. The parametric model of the bearing is therefore defined as a custom bearing with standard nominal dimensions and lower value of the radial internal clearance, securing a better match of the bearing inner flexible ring to the shape of the WG without the occurrence of the undesirable blank areas in the bearing-WG surface contact area. The ball bearing in the presented research is, however, modeled without the physical separator between the rolling elements. The spacing between the balls is secured by the boundary conditions defined in the finite element analysis software.

FEA Model Definition
The material of the particular harmonic drive components are defined as steel, which is an isotropic material with defined Young's modulus E = 210,000 MPa and the Poisson's ratio µ = 0.3. Each component of the harmonic drive has to be strained just in the elastic range area, where the Hooke's law can be applied. All contact conditions between the HD components were defined by the Frictionless contact type without the friction consideration (no bonded contact type was used), hence the non-linear analyses were performed.
The parametric model of the harmonic drive unit was created using the most possible relative accuracy of the Creo/Parametric software (1 × 10 −6 ). This relative accuracy value allows virtual models to be created that are accurate to ten thousandths of millimeters depending on the maximum dimension of the particular part. This value of the relative accuracy is sufficient enough, because the real production accuracy of the harmonic drive most precise components, e.g., the gearing of the flexspline and circular spline, is manufactured in the range of hundredths of millimeters. The Ansys mesher works with even higher relative accuracy, which guarantees the required level of accuracy also by meshed model.
The gearing teeth are the most precise elements of the harmonic drive component model. The selected meshing process was performed using the rule to divide each tooth to 10 elements per tooth height and six elements per tooth thickness, independent of the tooth dimension defined in millimeters. In other words, the models of the examined harmonic drive units (HD 42 and HD 120) with the same teeth number have contained the same number of finite elements. The element size of the HD 42 model was approximately 0.035 mm and the element size of the HD 120 model was approximately 0.1 mm as shown in Figure 4. to 10 elements per tooth height and six elements per tooth thickness, independent of tooth dimension defined in millimeters. In other words, the models of the examined monic drive units (HD 42 and HD 120) with the same teeth number have contained same number of finite elements. The element size of the HD 42 model was approxima 0.035 mm and the element size of the HD 120 model was approximately 0.1 mm as sho in Figure 4. All HD component boundary conditions were defined by the FEA according to test procedure conditions described above. The lost motion was evaluated by a test plying the 4% of rated HD output torque to the FS measuring its angular deformation both directions with the WG locked and the CS fixed and the boundary conditions w adapted to these requirements.

FEA Analyses Variable Set Generation
The design of virtual experiments was carried out to obtain a representative set tested parameter combinations. The set of variables used by the analyses of the HD l motion was generated according the number of chosen input parameters influencing lost motion and subsequently the torsional backlash by the optimal Space Filling Des method. It is a Latin Hypercube Sampling (advanced form of the Monte Carlo sampl method that avoids clustering samples) optimized to better fill the parameter space. T scheme distributes the design parameters equally throughout the design space, while corners and/or mid-points are not necessarily included [22]. The method is very useful the limited computation time, which was the main reason for using it for the presen research.
The parameters influencing the HD lost motion connected with the ball bearing w flexible rings were the following: 1. the bearing inner ring inner diameter d (mm); 2. the bearing outer ring outer diameter D (mm); 3. the bearing inner radial clearance Gr (mm).
The parameters influencing the HD backlash connected with the FS and CS te All HD component boundary conditions were defined by the FEA according to the test procedure conditions described above. The lost motion was evaluated by a test applying the 4% of rated HD output torque to the FS measuring its angular deformation in both directions with the WG locked and the CS fixed and the boundary conditions were adapted to these requirements.

FEA Analyses Variable Set Generation
The design of virtual experiments was carried out to obtain a representative set of tested parameter combinations. The set of variables used by the analyses of the HD lost motion was generated according the number of chosen input parameters influencing the lost motion and subsequently the torsional backlash by the optimal Space Filling Design method. It is a Latin Hypercube Sampling (advanced form of the Monte Carlo sampling method that avoids clustering samples) optimized to better fill the parameter space. This scheme distributes the design parameters equally throughout the design space, while the corners and/or mid-points are not necessarily included [22]. The method is very useful for the limited computation time, which was the main reason for using it for the presented research.
The parameters influencing the HD lost motion connected with the ball bearing with flexible rings were the following: 1.
the bearing inner ring inner diameter d (mm); 2.
the bearing outer ring outer diameter D (mm); 3.
the bearing inner radial clearance G r (mm).
The parameters influencing the HD backlash connected with the FS and CS teeth spline were following ( The nominal and the limit dimensions of the chosen input parameters are Table 3. The Ansys workbench software generated, according to the group of five in rameters via the Latin hypercube method, 27 sets of design points for the HD 42 are listed with the calculated values of the lost motion ϕ in Table 4.  The nominal and the limit dimensions of the chosen input parameters are listed in Table 3. The Ansys workbench software generated, according to the group of five input parameters via the Latin hypercube method, 27 sets of design points for the HD 42, which are listed with the calculated values of the lost motion Φ in Table 4.  The Ansys workbench software generated, according to the group of five input parameters via the Latin hypercube method, 27 sets of design points for the HD 120, which are listed with the calculated values of the lost motion Φ in Table 5.

Input and Output Variables Correlation Analysis
The parameter correlation performs simulations based on a random sampling of the design space using Latin Hypercube sampling to identify the correlation between all chosen HD parameters. The correlation coefficient indicates whether there is a relationship between the variables and whether the relationship is a positive or a negative number. The closer the correlation coefficient is to the extremes (−1 or 1), the stronger the relationship [23].
The correlation of the input parameters (d, D, G r , r FS , r CS ) and the output parameter (Φ) was observed in the presented study by three methods:
Pearson's linear correlation (detects a linear relationship between two variables); 3.
Spearman's rank correlation (detects a monotonic relationship between two variables).

Linear Regression Model
A linear technique to modeling the relationship between a scalar response and one or more explanatory variables (dependent and independent) is known as linear regression. Simple linear regression is used when there is only one explanatory variable; multiple linear regression is used when there are more than one. Multiple linear regression is a specific example of generic linear models with only one dependent variable, and it is a generalization of simple linear regression with more than one dependent variable. The basic model for multiple linear regression is: for each observation i = 1, . . . , n.
In the formula above, there are n observations of one dependent variable and p independent variables in the calculation above. Thus, Y i denotes the i-th observation of the dependent variable, while X ij denotes the i-th observation of the j-th independent variable, with j = 1, 2, . . . , p. The values j denote the estimated parameters, while i is the i-th independently identically distributed normal error [24].

Pearson's Linear Correlation
Pearson's linear correlation product is the Pearson's correlation coefficient used to measure the strength of a linear association between two variables, where the value r = 1 means a perfect positive correlation and the value r = −1 means a perfect negative correlation [25]. Pearson's correlation coefficient, when applied to a sample, is commonly represented by r xy and defined as: where: n-sample size; x i , y i -individual sample points indexed with i.
The sample means in the Equation (6) are defined as: and:

Spearman's Linear Correlation
The Spearman's rank correlation coefficient represents a statistical dependence between the rankings of two variables. It determines how well a monotonic function can describe the relationship between two variables. The Pearson correlation between the rank values of two variables is equivalent to the Spearman correlation between those two variables [25].
If all n ranks are distinct integers, it can be commonly represented by r s using the formula: where: n-sample size (number of observations).
The difference d i between two ranks (rg) of each observation rg(X i ) and rg(Y i ) is defined as: The results of the data correlation applying Pearson and Spearman for the HD 42 are shown in Figure 6.

Spearman's Linear Correlation
The Spearman's rank correlation coefficient represents a statistical dependence tween the rankings of two variables. It determines how well a monotonic function describe the relationship between two variables. The Pearson correlation between the ra values of two variables is equivalent to the Spearman correlation between those two v iables [25].
If all n ranks are distinct integers, it can be commonly represented by rs using formula: ( ) The difference di between two ranks (rg) of each observation rg(Xi) and rg(Yi) is fined as: The results of the data correlation applying Pearson and Spearman for the HD 42 shown in Figure 6. The results of the data correlation applying Pearson and Spearman for the HD are shown in Figure 7. The results of the data correlation applying Pearson and Spearman for the HD 120 are shown in Figure 7. The results of the data correlation applying Pearson and Spearman for the HD 120 are shown in Figure 7.

Global Sensitivity Analysis
The definition of the global sensitivity analysis (GSA) stated Saltelli et al. [14] as the study of how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input. For this research, the method of conditional variances was used, including the first path and the total effects method.
The global sensitivity analysis of the chosen variables was performed in steps, which were as follows: 1. Input and output dependence polynoma definition; 2. Quasi-random variables data generation; 3. First-order and total-effect Sobol indices calculation.
The second and the third steps were repeated in a cycle until the convergence by the Sobol indices calculation was reached.

Input and Output Dependence Polynoma Definition
The GSA requires the knowledge of the functional dependence between input and output parameters. The whole harmonic drive units FEA models contained approximately 1.5 million elements. The average time demand of one FEA model variant calculation was approximately 24 h. Hence, it was not realistic to calculate thousands of input parameter combination variants. The created polynomial function defining the dependence between the variable input parameters and the harmonic lost motion value was used by the generation of thousands of inputs required for the next calculation of Sobol's indices.
A second degree polynomial was found using the least squares method, which is described for the lost motion ϕ both gearboxes HD 42 and HD 120 for the n = 5 chosen parameters by the equation: where: A, B, C, D-polynoma coefficents; x-variables (d, D, Gr, rFS, rCS); n-number of variables;

Global Sensitivity Analysis
The definition of the global sensitivity analysis (GSA) stated Saltelli et al. [14] as the study of how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input. For this research, the method of conditional variances was used, including the first path and the total effects method.
The global sensitivity analysis of the chosen variables was performed in steps, which were as follows:
First-order and total-effect Sobol indices calculation.
The second and the third steps were repeated in a cycle until the convergence by the Sobol indices calculation was reached.

Input and Output Dependence Polynoma Definition
The GSA requires the knowledge of the functional dependence between input and output parameters. The whole harmonic drive units FEA models contained approximately 1.5 million elements. The average time demand of one FEA model variant calculation was approximately 24 h. Hence, it was not realistic to calculate thousands of input parameter combination variants. The created polynomial function defining the dependence between the variable input parameters and the harmonic lost motion value was used by the generation of thousands of inputs required for the next calculation of Sobol's indices.
A second degree polynomial was found using the least squares method, which is described for the lost motion Φ both gearboxes HD 42 and HD 120 for the n = 5 chosen parameters by the equation: where: A, B, C, D-polynoma coefficents; x-variables (d, D, G r , r FS , r CS ); n-number of variables; i = 1 . . . n-first variable index; j = i + 1 . . . n-second variable index.
According the control of the presented considerations, the input data from the FEanalysis data sets were put to the obtained equation and the output results were compared. The quadratic determination correlation coefficient for the presented equation reached the value of R 2 = 0.98 for HD 42 and R 2 = 0.99 for HD 120.
The FEA results and the fitted curve function results comparison is shown for the HD 42 in Figure 8 for the 27 data sets described in Table 4. According the control of the presented considerations, the input data from the FEanalysis data sets were put to the obtained equation and the output results were compared. The quadratic determination correlation coefficient for the presented equation reached the value of R 2 = 0.98 for HD 42 and R 2 = 0.99 for HD 120.
The FEA results and the fitted curve function results comparison is shown for the HD 42 in Figure 8 for the 27 data sets described in Table 4. The FEA results and fitted curve function results comparison is shown for the HD 120 in Figure 9 for the 27 data sets described in Table 5. The FEA results and fitted curve function results comparison is shown for the HD 120 in Figure 9 for the 27 data sets described in Table 5.

Quasi-Random Variables Data Generation
Quasi-random variables are not fully random in the sense that they cannot be predicted. An algorithm that generates low-discrepancy sequences must somehow bias the selection of new points to keep them away from the points already existing in order to

Quasi-Random Variables Data Generation
Quasi-random variables are not fully random in the sense that they cannot be predicted. An algorithm that generates low-discrepancy sequences must somehow bias the selection of new points to keep them away from the points already existing in order to maintain an even spread of points. They are, nonetheless, similar to random points in that they are uniformly dispersed across the whole variable space as described in [14]. The quasi-random variables can be calculated by various methods generating various data sets. The methods of quasi-random variables data generation used in the presented research contained: 1.
The Latin Hyper Cube method has a wide range of uses, but is mostly used in computer engineering. It was invented by McKay et al. from Los Alamos Laboratory in 1979. It has been proven that the Latin Hyper Cube sampling method can increase efficiency compared to Monte Carlo, but this only applies to certain classes of functions [26].
Some quasi Monte Carlo methods work best with specially chosen sample sizes such as large prime numbers or powers of small prime numbers. Halton sequences can be used with any desired sample size. There may be no practical reason to require a richer set of sample sizes than, for example, powers of two. However, first time users may prefer powers of 10, or simply the ability to select any sample size they choose [27].
Sobol sequences are an example of quasi-random sequences. The Russian mathematician Ilya M. Sobol introduced these in 1967. These sequences use a base of two to form successively finer uniform partitions of the unit interval and then reorder the coordinates in each dimension. These sequences are digital (t, s) sequences in base 2, based on Niederreiter's former work. The best uniformity of distribution as n→∞, good distribution for fairly small initial sets and a very fast computational algorithm are the three main requirements, which are the aim of Sobol sequences [28].
The scatter plots of the generated quasi-random variables by the Latin Hyper Cube, Halton and Sobol methods are shown for both gearboxes HD 42 and HD120 in Figure 10.
ls 2021, 14, x FOR PEER REVIEW 16 sets. The methods of quasi-random variables data generation used in the presente search contained: 1. Latin Hyper Cube set; 2. Halton set; 3. Sobol set.
The Latin Hyper Cube method has a wide range of uses, but is mostly used in puter engineering. It was invented by McKay et al. from Los Alamos Laboratory in It has been proven that the Latin Hyper Cube sampling method can increase effici compared to Monte Carlo, but this only applies to certain classes of functions [26].
Some quasi Monte Carlo methods work best with specially chosen sample sizes as large prime numbers or powers of small prime numbers. Halton sequences can be with any desired sample size. There may be no practical reason to require a richer s sample sizes than, for example, powers of two. However, first time users may prefer ers of 10, or simply the ability to select any sample size they choose [27].
Sobol sequences are an example of quasi-random sequences. The Russian math tician Ilya M. Sobol introduced these in 1967. These sequences use a base of two to successively finer uniform partitions of the unit interval and then reorder the coordi in each dimension. These sequences are digital (t, s) sequences in base 2, based on Ni reiter's former work. The best uniformity of distribution as n→∞, good distributio fairly small initial sets and a very fast computational algorithm are the three main req ments, which are the aim of Sobol sequences [28].
The scatter plots of the generated quasi-random variables by the Latin Hyper C Halton and Sobol methods are shown for both gearboxes HD 42 and HD120 in Figur

First-Order and Total-Effect Sobol's Indices Calculation
Saltelli et al. [14] define the first-order sensitivity index Si of Xi on Y accordin equation:

First-Order and Total-Effect Sobol's Indices Calculation
Saltelli et al. [14] define the first-order sensitivity index S i of X i on Y according the equation: where the conditional variance: is called the first-order effect of X i on Y. S i is a number always between 0 and 1. A high value signals an important variable. Total effects are a direct consequence of Sobol's variance decomposition approach and estimation procedure. The total effect index accounts for the total contribution to the output variation due to factor X i , i.e., its first-order effect plus all higher-order effects due to interactions [14].
The lost motion values of the fitted function calculated for the Sobol data set are shown for the HD 42 in Figure 11. The lost motion values of the fitted function calculated for the Sobol data set are shown for the HD 120 in Figure 12. The lost motion values of the fitted function calculated for the Sobol data set are shown for the HD 120 in Figure 12. The data presented in Figures 11 and 12 serve as the base for the calculation of the B and C matrices necessary to the Sobol's indices calculations. The cyclic procedure The data presented in Figures 11 and 12 serve as the base for the calculation of the A, B and C matrices necessary to the Sobol's indices calculations. The cyclic procedure of their calculation started with a small number of the quasi-random variables, which was gradually increased until the convergence was reached.
The progress of the cyclic calculation of the first order and total effect Sobol's indices for the HD 42 are shown in Figure 13. The progress of the cyclic calculation of the first order and total effect Sobol's ind for the HD 120 are shown in Figure 14.

Discussion
This research investigated the influence of a change in a five-input parameter se the value of the harmonic drive (HD) lost motion (torsional backlash) by the prescri test procedure. The flexspline (FS) of the harmonic drive was loaded by the ±4% of the rated output torque by the locked wave generator (WG) and fixed circular spline (CS) the lost motion value of the flexspline was measured.
The overall results of the HD 42 global sensitivity analysis calculations are presen in Figure 15. The progress of the cyclic calculation of the first order and total effect Sobol's indices for the HD 120 are shown in Figure 14. The progress of the cyclic calculation of the first order and total effect Sobol's ind for the HD 120 are shown in Figure 14.

Discussion
This research investigated the influence of a change in a five-input parameter se the value of the harmonic drive (HD) lost motion (torsional backlash) by the prescri test procedure. The flexspline (FS) of the harmonic drive was loaded by the ±4% of the rated output torque by the locked wave generator (WG) and fixed circular spline (CS) the lost motion value of the flexspline was measured.
The overall results of the HD 42 global sensitivity analysis calculations are presen in Figure 15.

Discussion
This research investigated the influence of a change in a five-input parameter set to the value of the harmonic drive (HD) lost motion (torsional backlash) by the prescribed test procedure. The flexspline (FS) of the harmonic drive was loaded by the ±4% of the HD rated output torque by the locked wave generator (WG) and fixed circular spline (CS) and the lost motion value of the flexspline was measured.
The overall results of the HD 42 global sensitivity analysis calculations are presented in Figure 15. The overall results of the HD 120 global sensitivity analysis calculations are pre in Figure 16. Global sensitivity analysis results show the minor impact of the bearing dime variations on the value of the HD lost motion, or the torsional backlash, respective The bearing inner diameter d was considered to have positive deviations fro nominal dimension simulating the clearance contact between the WG outer shape a bearing inner ring. The bearing outer diameter D was considered with the negative ations given to the nominal dimension, simulating the clearance between the FS shape and the bearing outer ring. The bearing internal diameter d has an influence by the correlation of about 17.95% to 17.09% by the standard data regression metho the HD 42, or 20.62% to 20.33% for the HD 120, but this is a minor influence consi the Sobol's indices values.
The lower the bearing internal radial clearance, the better the deformation of ner bearing ring on the wave generator outer shape, and the outer bearing rin matches the FS tube inner surface better. This fact causes a better equidistant match WG outer shape and the FS internal surface. The radial clearance of the bearing wit The overall results of the HD 120 global sensitivity analysis calculations are presented in Figure 16. The overall results of the HD 120 global sensitivity analysis calculations are pres in Figure 16. Global sensitivity analysis results show the minor impact of the bearing dime variations on the value of the HD lost motion, or the torsional backlash, respectivel The bearing inner diameter d was considered to have positive deviations fro nominal dimension simulating the clearance contact between the WG outer shape a bearing inner ring. The bearing outer diameter D was considered with the negative ations given to the nominal dimension, simulating the clearance between the FS shape and the bearing outer ring. The bearing internal diameter d has an influence by the correlation of about 17.95% to 17.09% by the standard data regression metho the HD 42, or 20.62% to 20.33% for the HD 120, but this is a minor influence consid the Sobol's indices values.
The lower the bearing internal radial clearance, the better the deformation of t ner bearing ring on the wave generator outer shape, and the outer bearing rin matches the FS tube inner surface better. This fact causes a better equidistant match WG outer shape and the FS internal surface. The radial clearance of the bearing wit ible rings Gr has an influence rated by the correlation of about 20.23% to 32.17% Global sensitivity analysis results show the minor impact of the bearing dimensions variations on the value of the HD lost motion, or the torsional backlash, respectively.
The bearing inner diameter d was considered to have positive deviations from the nominal dimension simulating the clearance contact between the WG outer shape and the bearing inner ring. The bearing outer diameter D was considered with the negative deviations given to the nominal dimension, simulating the clearance between the FS inner shape and the bearing outer ring. The bearing internal diameter d has an influence rated by the correlation of about 17.95% to 17.09% by the standard data regression methods for the HD 42, or 20.62% to 20.33% for the HD 120, but this is a minor influence considering the Sobol's indices values.
The lower the bearing internal radial clearance, the better the deformation of the inner bearing ring on the wave generator outer shape, and the outer bearing ring also matches the FS tube inner surface better. This fact causes a better equidistant match of the WG outer shape and the FS internal surface. The radial clearance of the bearing with flexible rings G r has an influence rated by the correlation of about 20.23% to 32.17% by the standard data regression methods for the HD 42, or 10.98% to 23.87% for the HD 120, but this is a minor influence considering the Sobol's indices values.
The higher value of the bearing radial internal clearance causes the imperfection of the component group WG-bearing-FS surfaces in contact with undesirable blank areas without surface contact, especially on the WG minor axis affecting the incorrect FS and CS teeth position accompanied by the teeth interference. The interference, then, occurs not just in the WG minor axis region, but also in the region of the teeth contact path start and end inclined from the WG minor axis at the angle of approximately 60 • in both directions.
The bearing parameter group interconnectedness can cause the accumulation of their effect on the WG-bearing-FS surface match accuracy, and also their effect on the lost motion according to the parameter combination. The same value of the bearing internal clearance by various values of the bearing diameters causes significantly different impacts to the WG-bearing-FS surface match and also the lost motion. The various bearing parameter sets cause various influence to the HD lost motion and subsequently the torsional backlash value, which has an impact on the Sobol's indices values. The parameters of the bearing with flexible rings (d, D and G r ) are intertwined, which influences the lower coefficients by the outer bearing ring diameter D for the HD 42 and similar values to the inner bearing diameter d for the HD 120.
The total-effect Sobol's indices of the bearing parameter set for the HD 120 are similar (0.05999; 0.05448; 0.05198), because all of them affect the bearing radial clearance and subsequently the HD component surface match and also the lost motion value in the same way. The total-effect Sobol's indices of the bearing parameter set for the HD 42 are different from the ones calculated for the HD 120. This is caused by the variation in the chosen ISO tolerance grade IT6, defined for various HD nominal dimensions (0.016 mm for the HD 42 versus 0.022 mm for the HD 120) which generates various bearing clearance values (0.02 mm for the HD 42 versus 0.036 mm for the HD 120).

Conclusions
Global sensitivity analysis by both HD variants let to the same parameter affecting the lost motion, and subsequently the torsional backlash the most; the offsets from the nominal FS and CS teeth shape r FS and r CS . They have the largest size among the observed parameters group and thus their criteria weight is highest too. Their Sobol's indices values are close to each other. The change in r FS and r CS affects the dimension over pin for the FS and between pin for the CS, and also the lost motion in the same way; equal changes in the offsets generate equal change in the output parameter, that is, the HD lost motion.
The highest correlation is observed for the HD 42 by the Spearman's rank correlation (67.58%), but the Sobol's indices total effect method is also pointing out a high level of correlation between 51.52% (FS offset), or the 55.13% (CS offset), respectively. The first order Sobol's indices show just the correlation of 29.02% (FS offset) to 27.93% (CS offset), respectively. The highest correlation is observed for the HD 120 by the Spearman's rank correlation (70.27%), but the total effect Sobol's indices also show a high level of correlation of 42.34% (FS offset), or 41.68% (CS offset), respectively. The first order Sobol's indices show just the correlation of 39.62% (FS offset) to 38.38% (CS offset), respectively.
The different values of the Sobol's indices found with the different offsets are the consequence of the various diameters of the examined harmonic drive units. They have the same teeth number and gear ratio, but the modules are quite different: the tooth height of the HD 42 is only 0.337 mm, while the tooth height of the HD 120 is 1 mm. The flexspline thickness affecting the torsional stiffness is also different for both HD units and thus, the values of the lost motion performed by the FEA of the HD 42 and HD 120 are not equal, which has also impact to the different values of the Sobol's indices.
The HD teeth spline offset from the nominal shape has the highest impact to the HD lost motion value by both virtual tested units. The proper selection of the r FS and r CS values will lead to the required low value of the HD backlash during its operation in the desired motion state.