Copula method of the kinematic accuracy of cylindrical roller bearing

: The kinematic accuracy of cylindrical rollers bearing isn’t only influenced by the dimensional and form precision of its parts, but also more influenced by the cooperation among them. For cylindrical rollers bearing, the main indices of the kinematic accuracy are the ending beat and the radial run out of bearings. There must be dependent relationship between the cooperative action among the parts of bearings and the ending beat or the radial run out of bearings. This relationship is hardly expressed by mathematical formula. However, because the parts dimension deviation and the run out of bearings follow statistic law and have their distribution characteristics, while the copula joint distribution function can connect the run out of bearings with multiple parts dimension accuracy, the copula function is introduced to analyze the dependent relationship between bearing parts cooperative action and its kinematic accuracy. Based on the copula function of statistic theory, the kinematic accuracy of bearings can be forecast by parts geometric accuracy.


Introduction
The kinematic accuracy of bearings isn't only restricted by its parts dimension precision but also influenced by the cooperation among them. In recent years, many researchers have been studying how to improving the kinematic accuracy from various aspects, example as force, structure, assembly, lubrication and etc. There were many works about the dependent connection between the parts precision accuracy and the bearings kinematic accuracy reported. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. But, very few scholars studied the bearing kinematic accuracy from statistical standpoint. For cylindrical rollers bearing, the radial run out and the ending beat of bearings is the most important indices influencing on the kinematic accuracy of bearings which often directly affect machine performance, especially high-speed precision instrument [16]. The radial run out data and the ending beat data of bearings are statistics1 from measurement and obey statistical law. their distribution characters are decided synthetically by bearing parts precision accuracy which owe to many conditions such as worker technique, environment temperature and humidity, the precision grade of machining tool and so on. On account of various 1 This work is supported by the National Key R&D Program of China (No. 2019YFB2005004) causation, every part dimension and form couldn't be standard value and has more or less deviation. While, the bearings is a assembly which consist of multiple parts such as rollers, cage, outer ring , inner ring and etc. When assembling parts together into the bearings, the parts were taken randomly and these parts precision synthetically decide the bearings kinematic accuracy. So, the kinematic accuracy of the bearings is integral effect. It needn't only to study the parts precision, but also to study the cooperation among them. While, the relationship between the kinematic accuracy of bearings and the components precision is quite complex and difficult to be expressed by mathematical model. However, because of their distribution characteristics in dimension, the Copula function of joint distribution can connect multiple variables distributions to construct the transmission and map relationship model between the multiple parts precision distribution and the bearing kinematic accuracy. So, the bearing kinematic accuracy can be known with the dimensional precision of parts before assembly.
In 1959, Sklar propounded the coupla theory firstly which can be applied to construct Copula function to join or "couple" multivariate distribution functions to their one-dimensional marginal distribution function [17]. In the end of the 20th century, the theory and methods of copula were rapidly developing at home and abroad. Since the 1980s, this theory has been widely applied to insurance, banking business, machine diagnostic system, even to buildings fields, traffic controlling, space technology and so on [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. General distribution algorithm can't build suitable multivariate JDF (joint distribution function) which can expresses the relationship between it and marginal single distribution functions. But, copula function can construct copula joint distribution function ( CJDF ) to join the multiple dimensional distribution function to the low dimensional marginal function. At the same time, comparing with BOA, PPCA and other algorithms, CJDF can be applied to non parametric estimation of dependence relation between stochastic variables with less operation and can be preferable to give expression to distribution situation of dominant group [33] In this paper, CJDF will be introduced to research distribution characteristics of the inter dependent linkages between kinematic accuracy and multiple variables distributions. With this ideas, the kinematic precision of bearings will be approximately estimated by establishing math model by measuring parts dimension deviation before assembly.

The meaning of CJDF
As the dependent framework between variables, Copula function includes almost all reliant messages of stochastic variables. When the relationship between variables can't be determined by traditional method, the copula function can be utilized to analyze relevant relationship between variables. CJDF, also called joint distribution function or distribution reliant function, is able to couple the joint distribution function of two or more variables with one dimensional single distribution functions and to distinguish reliant framework of the marginal property from multiple variables distribution. This notion stemmed from Sklar ideas: CJDF is able to be split into one Copula joint function and multiple marginal distribution functions and can be applied to analyze dependence among variables [34].

Constructing CJDF procedure
Constructing CJDF need two steps. Firstly, there need one joint distribution function consisted of two or more variables. Secondly, selecting one appropriate copula function by which the joint distribution function of multiple variables is connected to marginal univariant distribution function [35]: Step one: Joint distribution function Given random variables xi∈ R ( i=1,2,…, n) with continuous marginal univariant distribution function Yi(xi), which n-dimensional joint distribution function can be set as G(x1,x2,…,xn), (x1, x2,…, xn)∈ Rn. These exists one unique Copula function written as followed [36]: Step two: The Type selection of Copula function From the meaning of CJDF, a lot of functions are capable of being used as Copula function, of which there were two kinds functions most researched: ACF(Archimedean Copula function) and ECF(elliptic Copula function) [36]. The most distinct property of the ECP is that its variable has the same type distribution, moreover, the most common type of the ECP is multiple variables Gaussian Copula (MVG) and multiple variables student's Copula (MVT) [36].
According to the equation (2), multivariate Gaussian Copula function can be written as followed: The density function of (3) is written as followed: In building JDF of stochastic variables, CJDF have much superiority: 1)CJDF can be used to build flexibly JDF without limiting marginal function type [37]; 2)when structuring math model by CJDF, marginal distribution function indicates only individual messages of single variable, while the correlative messages among variables is able to be described by CJDF. Therefore, one dimensional distribution of variable and correlation between variables can be researched individually [37].

Influencing factors on the kinematic accuracy of Cylindrical roller bearing
There are many influencing factors on the kinematic accuracy of cylindrical roller bearing. Besides load and machine precision, the kinematic accuracy of the bearing is mainly decided by parts geometrical deviation (dimension deviation of rollers,dimension and form deviation of inner ring, dimension and form deviation of outer ring), radial internal clearance and the number of rollers [38]. When cylindrical bearing running, the kinematic accuracy can also be affected by lubrication, working loading, the assembled pattern, surroundings temperature and pretightening force. In these factors, the parts geometrical deviation is the most important influencing factor. The ending beat and radial run out of bearing are the vital indices of the kinematic accuracy of roller bearings. There must be a certain relationship between the kinematic accuracy and the dimension deviation in the radial and in the axial . To ensure the kinematic accuracy of running parts, the roller bearing need to be controlled in the radial run out and the ending beat of bearings to a certain range. So, in the paper, the relation between the radial and ending beat of bearings and the dimension deviation of bearing parts will be analyzed to control the kinematic accuracy of bearings.
When the parts has finished, the dimension of them has a certain deviation in precision range because of accidental factors such as temperature, machine accuracy, technique of worker and etc. The deviation value is stochastic and follows the statistic law. Before assembling, the parts of bearings are independent in dimension distribution and self-interrelated. While, after being assembled, they can influence and interact on each other. This interact would cause bearing radial run out and ending beat and then influence kinematic accuracy. The influence on the kinematic accuracy is not single action of one part, but synthetic effect of two or more parts. So, the value of the radial run out and ending beat is also stochastic and follows the statistic law, which must have a relation with the joint distributions of parts. The relationship can be studied by the copula theory.

The stochastic variables
From the above section, the radial run out and ending beat of bearing is the main indices influencing on bearing kinematic accuracy. The influence is a integral effect of all bearing parts working. So, without consideration of the influence of other elements, only three elements ( outer ring, inner ring and the rollers) will be taken into account to analyze the kinematic accuracy of bearings. The dimension deviation of them can be measured by precision instrument, so do the the radial run outs and ending beat of bearings. With the statistic theory, the dimension distribution law of them can be gotten by the statistic method and the relation between them can be known by constructing the mathematics model among them with the copula theory and method. Then , the kinematic accuracy of the bearings can be known before assembling.
Take the cylindrical roller bearings (the type of bearing is NU1004) for example. Let variable y1 represented the radial run-out and variable y2 represented the ending beat of bearing respectively. For inner ring , let stochastic variable u represented the upper deviation of outer ring and for outer ring, let stochastic variable v represented the lower deviation of inner ring respectively. Because every bearing consists of a number of rollers, the dimension uniformity of rollers is very important factor influencing on the kinematic accuracy of bearings. So, let stochastic variable w represented the absolute value of the difference which is the biggest deviation minus smallest deviation of rollers.
Took 150 groups parts of the bearing to be marked by the number and measured them respectively and gain 150 groups data about u, v, and w. then, these 150 groups parts were assembled together by the No.and 150 bearings were gained. At the same time, the radial run out and ending beats of bearings were measured in the exclusive instrument. [38] and 150 groups data about y1, y2 were gotten. Only the 10 groups data measured of y1, y2, and u, v, w list in table 1 because of the limited space. Even if the measuring instrument is very precise, it has error itself. The 150 groups data had been processed by the error separation method [38]. The next analysis come from processed data.

parameters estimation
To search for the statistical law of the upper deviation u of the outer ring diameter, the lower deviation v of inner ring diameter and the absolute w of dimension difference of rollers, the statistic software R can be utilized to plot histograms of u , v and w. the deviation distribution of u , v and w were displayed in the histograms forms shown in Figure 1. The dotted curves in the Figure 1 indicated the kernel density of u , v and w. In the same way, the deviation distribution of y1 and y2 were displayed in the histograms forms shown in Figure 2. The dotted curves in the Figure 2 indicated the kernel density of y1 and y2. From the figure 1 histogram and figure 2 histogram, the five stochastic variables u, v, w and y1, y2 all approximately complied with normal distribution, but that is only a hypothesis. With the software R, their distributions law had been tested and verified. The results indicated their normal distribution can't be refused as shown in table.2.
In table 2,  is mean value,  is standard deviation,  ,  are estimation of mean value and standard deviation respectively, P is probability. In normal condition, according to statistics theory, if P>0.05, the above hypothesis is true, or else, P<=0.05, false.
From the table 2, the above hypothesis about normal distribution of u, v, w, and y1, y2 can't be refused.

Analyzing correlation between stochastic variables
Considering different level influence of different variable on the kinematic accuracy, the relativity among the above variables need to be analyzed that can be expressed by the covariance among the variables. With the software R, the covariance matrix between variables u, v, w and y1, y2 is shown in table 3. From table 3, y1 have a very strong relativity with u, v, w respectively and in the same way, y2 have very strong relativity with u, v, w respectively. Because of the deviation in dimension of bearings parts, the bearings working must cause the radial and axial beat of bearings. The bigger the u, v, w , the worse the beat. In addition, from the table 3, the upper deviation of outer ring has a litter stronger relativity with the radial run out than the lower deviation of inner ring, while the upper deviation of outer ring has obviously stronger relativity with the ending beat than the lower deviation of inner ring. The relativity between the rollers uniformity and the radial run out of bearings is the strongest relation comparing to the relativity between the outer ring and the inner ring. The relativity between the rollers uniformity and the ending beat of bearings is also the strongest one.
To understand intuitively the relativity between u, v, w and y1, y2. R-Plot can be used to draw the relativity dots plot between them as shown in Figure 3. From Figure 3, y1, y2 have stronger relations to w than the relation to u, v that illustrates the radial run out and the ending beat of bearings bear very stronger relativity to the uniformity of rollers. In the standard deviation range, the smaller w, the more the uniform and the smaller the beat.

variables distribution function
The radial run out of bearings results from the integral action of bearings parts. The distribution law of the radial run out must have connection with the joint distribution of u 、 v 、 w , then, the distribution function of the radial run out of bearings can be expressed by the joint multivariate distribution function of u 、 v 、 w . So does the ending beat of bearings. A multivariate distribution function is very difficultly to be described. According to the above theory, the CJDF can be utilized to build model to connect the joint distribution function of multiple variables with one-dimensional distribution function. So the distribution function of the radial run out of bearings can be described by the single dimensional marginal distribution function of the upper deviation of outer ring u , the lower deviation of inner ring v and the rollers uniformity w .
Set the single dimensional marginal distribution functions of the variables u 、 v 、 w as F(u) , F(v) , F(w) respectively. From the section 3.2.1, stochastic variables u 、 v 、and w approximately follow the normal distribution. According to mathematical statistics, the distribution function of u 、 v 、and w can be expressed as: values and the standard deviation of u 、 v 、and w respectively, moreover, the density function of u 、 v 、 w can be expressed respectively as: Set the distribution function of the radial run out of bearings as ) y ( F 1 and let w) v, F(u, indicate the joint distribution function of variables u、v and w. From the above stated, there must exist the relation between the single dimensional distribution function ) y ( F 1 and multiple dimensional distribution function w) v, F (u, . Similarly, set the distribution function of the ending beat of bearings as ) y ( F 2 . there also exists the relation between the ending beat distribution function ) y ( . The relations can be expressed in formula as: The high dimensional distribution function w) v, F(u, is very complicated to be solved. In general, many mathematics software such as MATLAB and R software only provides two-dimensional copula functions. To be convenient for analysis, the three dimensional joint distribution function w) v, F(u, can be simplified into three items, each of item is one pair-wise joint distribution function containing two-dimension stochastic variables and multiplied by a parameter i [39]. Three items of w) v, F(u, can be written as three two-dimensional distribution function v) Constructing the math modeling is followed as： 3) are revising parameters, and moreover, the density function modeling of (8) are expressed as followed: Where uv  is the correlation coefficient variables u and v , vw  is the correlation coefficient varia bles v and w , wu  is the correlation coefficient variables w and u .

Constructing copula function
From the last section, the CJDF of three dimension can be solved by dimension reduction. The pair-copula function can use a very flexible and intuitive way to construct high dimensional copula function [40]. the pair-wise joint distribution function v) F(u, can be easily connected to its one-dimensional marginal distribution functions . There exists only one copula function F(v)) C(F(u), set as: From the section 4.2.1, u 、 v 、 w and y1 、y2 satisfied the normal distribution and follow the normal distribution law. From the section 2.2. Gauss copula function can be used to construct mathematical model between the joint distribution function and its marginal distribution function. Now, for Eq. (10) and Eq. (11), each of equation has a Gaussian copula function ,respectively written as: Where    is two-dimensional normal distribution function, F -1 (u) is the inverse function of F(u), F -1 (v) is the inverse function of F(v), F -1 (w) is the inverse function of F(w). The Eq. (13) is spread as followed: The density functions of Gaussian copula functions respectively. So, according to statistics theory the Where:

Parameters estimating
For the specified distribution of stochastic variable, the maximum likelihood method is a widely used parameter estimation method [19]. Therefore, The method of Maximum Likelihood can be introduced to estimate parameters of (15), (16). For easily resolving it, take the logarithms on both sides of the equation and log-likelihood function of the joint distribution function can be gotten as:  (9), the density function can be spread as: From the section 3.2.1, The variables y1, y2 follow the normal distribution and their distribution function and density function can be expressed as: the expression Eq. (20 ). can be spread as: Similarly, Eq. (6) can be spread as:     (8) and the math modeling of the distribution of the radial run-out of the bearings and the ending beat can be obtained as follows: From the coefficient of the above mathematics modeling Eq. (24)and Eq. (25), the interaction with each other between parts of bearings can produce different effect on radial run out and ending beat of bearings. For the first equation of Eq. (24) or Eq. (25), the coefficient of the first item was the biggest one that was 0.2836 (max (0.2836, 0.0028, 0.0097)), that meant the interaction with each other between the inner ring and outer ring can produce the greatest effect on radial run out, then, the effect between the rollers and inner ring, the last, the effect between the rollers and the outer ring. For the second equation of Eq. (24) or Eq. (25), the coefficient of the second item was the biggest one that was 0.0735 (max (0.0086, 0.0735, 0.0418)), that meant the interaction with each other between the inner ring and the rollers can produce the biggest effect on the the ending beat, then , the effect between the rollers and the outer ring, the last, the effect between the outer ring and the inner ring. This shows that when bearing parts are being processed, the upper deviation of outer ring and the lower deviation of inner ring need to be controlled strictly besides controlling dimension in standard tolerance, simultaneously, the rollers must also be strictly selected to guarantee their uniformity.
In general, it is after the parts are assembled into the bearings that can the radial run-out and the ending beat of bearings come into being. According to the above method, as long as the dimensions precision of inner ring, outer ring and rollers are known , the radial and end beat range of the bearing can be calculated before assembling. This meant once the parts precision is determined, the kinematic accuracy of bearing will be sure.
Of course, the above mathematics modeling were built according to the data from NU1004 bearings. For other type bearings, the distribution law of the dimension of parts may be different, the mathematics model may be also different. It need construct corresponding modeling from the relevant parts. But the copula method can be used for any bearings.

Experimental Verifying
To verify the above conclusion, we have done experiment with the same type bearings (NU1004 ). the other 150 groups parts of the bearings were taken randomly, numbered, measured respectively and the other 150 groups data of u、v、w were obtained by the number. Similarly, because of the limited space, only the 10 groups data were listed in table.4. Substituting this 150 groups data into the first equation of (25) and solving it and the mean values 1  and the standard deviation 1  can be gotten. In the same way, Substituting this 150 groups data into the second equation of (25) and solving it and the mean values 2  and the standard deviation 2  can be yielded. According to the parameters, it can be known that the radial run out and the ending beat of bearings will be in a specific scope after assembling. To verify this result, the 150 groups parts were assembled according to their No. After the radial run out and the ending beat of bearings were measured, the mean values  and the standard deviation  can be evaluated by statistic software R. the computing value by (25) and the evaluated value were list in table.5. From table.5, the computing value by modeling (25) and the evaluating value by measured data is approximately unanimous. There was only difference from 0.0001m to 0.0003m between two results.

Conclusions
(1) the coordination among the parts of bearings is very important. The kinematic accuracy of bearings mainly come from the interaction with multiple parts. In general, it is often considered that the high grade precision bearings consist of high precision parts. But that is not the case, the harmonization among its parts is more important.
(2) From the mathematics modeling in the paper, for the radial run out, the harmonization between the inner ring and the outer ring is very important, while for the ending beat, the harmonization between the rollers uniformity and the inner ring should be paid more attention. Therefore, when the radial beat accuracy of bearings must be requested, the harmonization between the inner ring and the outer ring have to be controlled strictly. When the face beat accuracy of bearings must be requested, the harmonization between the rollers uniformity and the inner ring need to be controlled strictly.
(3) Bearing is a system consisted of multiple components. The radial run out and the ending beat are main indices of the kinematic accuracy. The interact with each other among parts exerts joint influence on the radial run-out and the ending beat of bearings and then on kinematic accuracy level. The method of distribution estimation algorithm with CJDF was be introduced to constructing the math modeling between the precision of parts and the radial and ending beat of bearing. With this modeling, the kinematic accuracy of bearing can be forecast by the actual dimension deviation of parts.