Statistical Analysis of Industrial Grinding Brush Force Characteristics Based on Finite Element Approach

Abstract

This paper presents an in-depth study of the helical grinding brush force characteristics aiming at developing a mobile robot system to perform rust removal and other surface processing tasks. Based on an off-line Finite Element model that can calculate brush filament deformation and force behaviors, a mathematical regression model has been developed to summarize brush force changes subjected to varying conditions into a series of mathematical equations. The predictions of the mathematical model are well converged with the Finite Element modeled results and the R-squared value is up to 0.95. The paper presents the model form and calibrated coefficients, which may provide an advantageous tool to predict the grinding brush force changes in real time and contribute well to an automatic grinding control application.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

With rapid development of industrial automation, surface treatment processes such as grinding, polishing, deburring, rust, and paint removal have become emerging robot utilization areas and received considerable research attentions [1–4]. Steel grinding brushes are currently widely used on hand hold grinders to perform these tasks without comprehensive control techniques. Intelligent grinding robots with brush tools are still rare in the market [5, 6] due to difficulties involved in the end effector control, and more often the surface processing machines are large equipment with other effectors such as those based on magnetic assisted finishing methods [7–9].

This study addresses analysis on the grinding brush force characteristics aiming at design of small mobile robotic grinding system. The goal is to eventually develop a mathematical end effector model that can real time predict grinding brush force under varying conditions so that the process can be controlled based on the model outputs.

Studies on brush tools so far include empirical researches on tooth brushes [10, 11], experimental and mathematical analysis of abrasive brushes [12–14], Finite Element (FE) analysis of the behaviors of road sweeping brush with straight rectangular tines [15–17], and recently exploratory analysis on helical brush filament characteristics using the FE method [18]. The studies on the road sweeping brushes [15–17] have originated and developed the methodology of modeling brush deformation using Finite Element method to analyze brush contact geometry and contact force. Brush tines were meshed with ANAYS BEAM-189 elements and the deformation subjected to position constraint, rotational effect, and frictional effect was calculated using iterative approach [15]. The results were further summarized into mathematical models [16] and utilized in experiments and dynamics analysis [17]. To expand the technique for helical grinding brush a new model for helical filament analysis has been established [18].

This study is focused on the cup type grinding brush with helical filaments as shown in Figure 1. In the previous helical brush analysis [18], FE platform ANSYS and APDL programming were utilized to create models for varying filament geometries and solved for different operating variables. However, ANSYS models can only be run off-line which is not competent for real-time control purpose. In the earlier sweeping brush study [16], a regression modeling approach has been developed which can summarize the FE modeling results into a mathematical model with limited number of model coefficients. This methodology can be employed in the current study to create statistical model for helical grinding brush, so that further control applications can directly utilize the mathematical model outputs in a real-time environment.

[figure omitted; refer to PDF]

2. Finite Element Modeling and Results Validation

The study firstly assumes that the brush filaments are equal in geometry and can be investigated using one example. The important geometric parameters for modeling include filament length L, filament radius r, radius of helical lap R, mounting inclination angle φ, distance between neighboring laps δ, and the number of laps N, as shown in Figure 2(a). The model firstly generates key points with the specified parameters and then generates the helical line using APDL SPLINE commend through the key points. The helical line is then meshed with 100 BEAM-189 elements, each has 3 nodes, and each node has 6 degrees of freedom.

[figures omitted; refer to PDF]

As shown in Figure 2(b), when the brush is rotating and compressed onto a target surface, there are two operational parameters, namely, the rotational speed ω and the vertical penetration Δ. The loading conditions applied to a helical filament include full constraints at the top, i.e., UX=0, UY=0, UZ=0, ROTX=0, ROTY=0, ROTZ=0, at the top node, the penetration Δ applied in an upward manner, i.e., UY=Δ at the bottom node, frictional force f on the tip, and centrifugal force $F_c$ on the filament body, which are applied in a manner as explained below. The end support force Fv is the value to be found from the modeling results.

Since the filament exhibits large deformation, the location of a specific point on the filament body will not be the same as that before deformation, which affects the directions and parameters of the centrifugal force and the frictional force. A ramped iterative loading approach has been used to solve the deformation equilibrium in multiple steps, and the geometry results from the last step are retrieved using ANSYS APDL codes to update geometry and parameters for solving the next equilibrium step.

The centrifugal force on a specific node in the i^th step in the ANSYS global coordinate system as shown in Figure 2(b) can be written as$$\begin{array}{}\text{(1)}& F{c}_i=m·r_i·{\omega }^{\mathrm{2}}=m·{\omega }^{\mathrm{2}}·\sqrt{{x_{i-\mathrm{1}}}^{\mathrm{2}}+{z_{i-\mathrm{1}}}^{\mathrm{2}}}\end{array}$$where m represents the mass of an element containing this node, ω is the brush rotational speed, and $x_{i-1}$ and $z_{i-1}$ are the node coordinate parameters retrieved from the previous step. The frictional effects on filament tips are assumed to be analyzable using frictional coefficient u, although more complex abrasive effects should be investigated in the future. The frictional force on the filament tip in the i^th step can be given by$$\begin{array}{}\text{(2)}& f_i=\mu ·F{v}_{i-\mathrm{1}}\end{array}$$where $F{v}_{i-1}$ is the vertical contact force from the target surface retrieved from the previous step. The friction is along a tangential direction and it can be applied using two components as below$$\begin{array}{c}\text{(3)}& f_{ix}=\mu ·F{v}_{i-\mathrm{1}}·\sin \left[\arctan \left(\frac{z_{i-\mathrm{1}}}{x_{i-\mathrm{1}}}\right)\right]f_{iz}=\mu ·F{v}_{i-\mathrm{1}}·\cos \left[\arctan \left(\frac{z_{i-\mathrm{1}}}{x_{i-\mathrm{1}}}\right)\right]\end{array}$$The model solving procedure is a ramped and iterative process as developed in the previous study for sweeping brush [15]. The model gradually applies rotational speed in a number of substeps until the geometry of a free-rotating filament is obtained. Then the position constraint and frictional force are applied and solved by adding 1mm penetration in each substep. At each substep the APDL codes automatically recalculate the deformation and update geometric parameters until the maximum difference between neighboring iterations are no larger than 0.01%. The convergence of the final model is therefore ensured automatically by the APDL codes. An example showing the modeled filament deformation and stress is shown in Figure 3.

[figure omitted; refer to PDF]

By changing parameters in the model using ANSYS APDL codes, results have been obtained and statistically compared to assess the effects of different filament geometric variables and operating parameters on the characteristics of the brush tool, as shown in Figure 4. The results were obtained for a 22mm long helical filament with 3 laps and the filament radius is 0.15mm.

[figure omitted; refer to PDF]

A test rig as shown in Figure 5 has been constructed in the lab to verify the FE modeling results.

[figure omitted; refer to PDF]

The test rig is equipped with four coin-sized FSR load cells below the test plate, and each FSR load cell consists of a flexible substrate with printed semiconductor and a flexible substrate with printed interdigitating electrodes, resulting in linearly reduced electric resistance when they are pressed. The grinding brush can be rotated by a 1KW permanent magnet DC motor and driven downward by a lead screw onto the test plate, and the vertical displacement can be recorded by reading the angle of rotation of the turning wheel. Using the test rig, the brushing force can be measured from the load cells and then plotted against different operating parameters such as rotational speed and vertical displacement, as presented in Figure 6. The grinding brush used in the test consists of 1200 steel helical filaments whose radius is 0.15mm and the average length is 22mm.

[figure omitted; refer to PDF]

The experimental results showed well consistent relationships comparing with the FE modeled results. The difference may be caused by vibration and interactions between neighboring filaments, which can be investigated in later studies. Above all, FE modeling results can be utilized to explore optimized brush design and grinding brush control issues.

3. Statistical Modeling and Results Validation

The FE model has to be calculated off-line and take remarkable CPU resources which make it unsuitable to use directly on a real time control system. In a foreseeable robot grinding system, the controller would need to control the contact force by varying displacement and rotational speed in real time. This results in the requirements to develop a mathematical regression model over the FE modeling results, so that the end effector force as shown in Figure 4 can be quickly predicted from mathematical equations.

In order to obtain enough datasets to construct a mathematical regression model, random changes in the geometric parameters and operating variables have been introduced, and the FE model was updated with these changes and rerun with APDL Marcos to eventually produce 5000 different datasets. The variables and the ranges of changes are listed in Table 1.

Table 1

List of FE model explanatory factors and the ranges of changes.

Variable	Range of changes	Description
L	20mm ~ 40mm	Filament length
r	0.1mm ~ 0.2mm	Filament radius
R	0.5mm ~ 2mm	Filament lap radius
N	2 ~ 4	Number of laps
φ	20deg ~ 60deg	Filament mounting angle
Δ	0mm ~ 15mm	Vertical penetration
ω	100rpm ~ 10000rpm	Brush rotational speed

Although the FE model produces extensive results such as stress, strain, deformed shapes, etc., the most interested value is the vertical contact force produced on the filament tip, since this dominates the grinding effects. The force values from the 5000 datasets were regarded as observations for further statistical modeling. In order to discover the distribution pattern, the observation values have been imported into SPASS for Quantile-Quantile (Q-Q) plot verifications, as shown in Figure 7.

[figure omitted; refer to PDF]

As demonstrated in Figure 7 and similar to the previous research results on metal sweeping brush [16], the contact force values of helical grinding brush filament best obey Gamma distribution, which has below probability density function$$\begin{array}{}\text{(4)}& f\left(x,k,\theta \right)=\frac{x^{k-\mathrm{1}}{e}^{-x/\theta }}{{\theta }^k\Gamma \left(k\right)}\end{array}$$where k is a shape parameter and θ is a scale parameter. The likelihood function for observations (x₁, x₂,…x_N) is$$\begin{array}{}\text{(5)}& L\left(k,\theta \right)={{\displaystyle \prod }}_{i=\mathrm{1}}^{N}f\left(x_i,k,\theta \right)\end{array}$$Based on the maximum likelihood estimation (MLE) approach, the log-likelihood function that should be maximized can be derived as follows:$$\begin{array}{}\text{(6)}& l\left(k,\theta \right)=\left(k-\mathrm{1}\right){{\displaystyle \sum }}_{i=\mathrm{1}}^N\ln \left(x_i\right)-{{\displaystyle \sum }}_{i=\mathrm{1}}^N\frac{x_i}{\theta }-Nk\ln \left(\theta \right)-N\ln \left(\Gamma \left(k\right)\right)\end{array}$$For model calibration purpose, the log-likelihood value was calculated for each individual force observation, so that the variations of explanatory factors in each dataset can be taken into account. The individual log-likelihood function is given by $$\begin{array}{}\text{(7)}& l_i\left(k,\theta \right)=\left(k-\mathrm{1}\right)\ln \left(x_i\right)-\frac{x_i}{\theta }-k\ln \left(\theta \right)-\ln \left(\Gamma \left(k\right)\right)\end{array}$$where x_i is the predicted filament force, k is the observed filament force, and θ is the scale factor that can be determined from calibration.

The regression model forms often mentioned in literatures [19, 20] include polynomial regression, logistic regression, multiplicative model, hazard function model, etc. The model suitability varies in different applications whilst for brush tool analysis purpose, the previous study [16] has addressed that the multiplicative model form is superior, for its less coefficients and easily configurable model functions that can quickly reveal the contribution of each explanatory variable. The general multiplicative model form can be given in (8), where f_i (x_i) represents an individual subfunction for one of the model explanatory factors, such as filament length, rotational speed, etc.$$\begin{array}{}\text{(8)}& f=f_{\mathrm{1}}\left(x_{\mathrm{1}}\right).f_{\mathrm{2}}\left(x_{\mathrm{2}}\right).f_{\mathrm{3}}\left(x_{\mathrm{3}}\right)\text{...}\end{array}$$By investigating the effect of each variable as demonstrated in Figure 4 and following comprehensive calibration testing, the final model form has been determined to be$$\begin{array}{}\text{(9)}& f=C·\left(\mathrm{1}+k_L·{X_L}^{n_L}\right)·\left(\mathrm{1}+k_r·{X_r}^{n_r}\right)·\left(\mathrm{1}+k_R·X_R\right)·\left(\mathrm{1}+k_{\phi }·X_{\phi }\right)·\left(\mathrm{1}+k_N·X_N\right)·\left(\mathrm{1}+k_{\omega }·{X_{\omega }}^{n_{\omega }}\right)·\left(\mathrm{1}+k_{\Delta \mathrm{1}}·\left(X_{\Delta /L}-k_{\Delta \mathrm{2}}·X_{\omega }\right)\right)\end{array}$$where C is a leading constant and k and n denote the gradient and exponent, respectively, with respect to the seven explanatory factors listed in Table 1, whilst X represents the explanatory attribute values processed or normalized in a form for calibration convenience.

Herein $X_{\Delta /L}$ indicates the processed attribute values are penetration (mm) divided by the filament lengths (mm), which are therefore dimensionless values. The rotational speed value $X_{\omega }$ is included twice in the equation and associated with coefficient $k_{\Delta \mathit{\text{2}}}$, which is because a revolving brush would deflect due to centrifugal effects and therefore reduce penetration on the target surface. The units of measure of the processed attribute values are shown in Table 2.

Table 2

Units of measure of the attribute values.

Attribute Value	Unit of measure
$$X_L$$	100 / mm
$$X_r$$	10$\ast$mm
$$X_R$$	mm
$$X_{\phi }$$	degree / 10
$$X_N$$	Integer 1, 2, 3, 4, 5
$$X_{\omega }$$	rpm / 10000
$$X_{\Delta /L}$$	mm / mm

The model calibration was carried out using the MLE approach. The final model equations in (9) were applied to each attribute dataset to calculate the predicted values, and the likelihood estimate for each dataset was then calculated using (7). The summed total of the likelihood values were then maximized using Newton-Raphson iteration method to find the most appropriate model coefficients. Model calibration was carried out with half of each data set, referred to as the calibration data, chosen by allocating a random number to each observation and then selecting half the records based on the value of the random number. The rest data were used for validation.

The model fitting performance was firstly evaluated against each explanatory factor, as the one-dimensional graphs presented in Figure 8. Since the dataset includes large number of random combinations of explanatory factors, the one-dimensional plot does not necessarily explain how each factor affects the model behavior. Instead, the trends can be noticed from the model coefficients, whilst the one-dimensional plots represent how the fitted model results look like from one angle of consideration. From Figure 8 it can be seen that the model well captures the effects of filament length, filament radius, and penetration. For the rotational speed, mounting angle, and filament radius, the trends are not clear in the graphs due to combined other effects; however the model fitting performance still demonstrates well converged results.

[figure omitted; refer to PDF]

In order to evaluate model fitting performance for different designs made by different manufacturers, since there is not a universal standard, five randomly selected filament samples have been selected and validation of predicted values has been carried out as presented in Figure 9. It can be seen that in most cases the FE modeled results and the predictions from the statistical model are well converged.

[figure omitted; refer to PDF]

The attribute values used for each case validation are presented in Table 4. In case 2 where the filament length is longer than usual, there is noticeable deviation in Figure 9, which indicates an imperfection in the model fit for long filament. This deviation may be associated with the less fitted section of the length graph in Figure 8, likely due to a constrained integer exponent used in the model function. However, the benefit of using simpler functions and fewer coefficients is to avoid an overfitted model which may be less useful for prediction purposes.

The model validation has also been carried out using X-Y plot for each observed, i.e., FE modeled, and predicted filament force values, as demonstrated in Figure 10. If the predictions of the statistical model are ideally accurate, each point should be distributed along a diagonal line. In Figure 10 the results are well distributed and the trend line is nearly diagonal. The R-squared value of the comparison results is 0.9589.

[figure omitted; refer to PDF]

The comparison of regression modeled brush force with FE modeled values and experimental results are shown in Figure 11.

[figure omitted; refer to PDF]

It can be seen that the consistency is good although the difference increases when the rotational speed is large, which is due to model fitting changes with rotational speed as can be seen from Figure 8. The comparison indicates a largest error 2.3% between the regression results and the FE modeled results. Nevertheless, in practice it is fairly reliable to use the regression model form (9) along with the calibrated coefficients in Table 3 to predict brush force in real time for a specific brush design and operation scenario.

Table 3

Model coefficient values.

Coefficient	Value	Coefficient	Value
θ	0.9572	$$k_{\phi }$$	-0.0809
C	3.44E-5	$$k_N$$	-1.1126
k L	0.0239	$$k_{\omega }$$	-0.3223
n L	3	$$n_{\omega }$$	2
k r	20.987	$$k_{\Delta 1}$$	133.657
n r	4	$$k_{\Delta 2}$$	-0.09687
k R	-0.3027

Table 4

Attribute values of individually validated filament samples.

Case	Parameters
1	φ=45°, r =0.15mm, L=24mm, N =3, R =1mm, ω=2000rpm
2	φ=45°, r =0.15mm, L=40mm, N =3, R =1mm, ω=2000rpm
3	φ=45°, r =0.18mm, L=28mm, N =4, R =1.9mm, ω=2700rpm
4	φ=43°, r =0.18mm, L=29mm, N =2, R =1.9mm, ω=7700rpm
5	φ=50°, r =0.14mm, L=20mm, N =3, R =0.7mm, ω=4900rpm

4. Outcomes from the Grinding Brush Force Analysis

During the modeling process, several patterns of the grinding brush force characteristics have been noticed and are addressed below.

(i)

The effect of penetration to the brush force is slightly nonlinear and affected by other factors. A multivariable regression model has captured most of the effects and for the penetration effect alone, a simple linear function can give well converged results.

(ii)

The rotational speed affects brush force by adding centrifugal effects to the filaments. If the penetration remains constant then a higher rotational speed results in lower contact force, since the centrifugal effects contribute more to the filament deflection.

(iii)

The radius of the filament has significant effect on the brush force and the regression model suggests a biquadratic relationship.

(iv)

The filament length affects brush force in a declined manner, i.e., longer filaments tend to soften the brush stiffness, since less filament deflection can produce the same amount of penetration.

(v)

The filament mounting angle, number of laps and radius of laps also affect brush force in declined manners, and these have been modeled using linear functions.

5. Conclusion

This study carries out analysis on the characteristics of helical abrasive brush filaments and presents a statistical regression model that can predict brush force in real time. The brush filament deformation was first calculated in ANSYS and validated through experiments as the observed data, and the statistical model was then calibrated using multiplicative model function and the maximum likelihood estimation (MLE) approach. The statistical model can predict well fitted results and the R-squared value is over 0.95. Comparison also indicates the largest error between the regression modeled results and the FE modeled results is within 2.3% for a commonly used grinding brush.

The model also decoupled the interactive effects between brush operating variables such as the penetration and the rotational speed. The relationships have been well explained by model coefficients and are valuable in automatic grinding control applications in future studies. Effects due to brush design factors such as filament radius, length, inclination angle, lap radius, and number of laps have all be included in the regression model, which can be used in optimal grinding brush selection or redesign. Altogether the statistical grinding brush model can make considerable contributions in the industrial grinding applications.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study is supported by the National Natural Science Foundation of China (Grant no. 61703192), Shandong Province Higher Institution Science and Technology Plan Project (Grant no. J17KA028), Major Scientific and Technological Innovation Project of Shandong Province (Grant no. 2017CXGC0907), Shandong Province Natural Science Foundation (Grant no. ZR2016FL13/ZR2014EL023), Scientific Research Fund of Liaocheng University (Grant no. 3180500/318011519), and the Special Fund Plan for Local Science and Technology Development Lead by Central Authority, China.

References

[1] M. Zhang, J. Wang, W. Zhang, X. Wang, "Robotic belt grinding method for the surface of whole propeller blade," Jiqiren/Robot, vol. 37 no. 3, pp. 318-368, 2015.

[2] X. Xie, L. Sun, "Force control based robotic grinding system and application," Proceedings of the 2016 12th World Congress on Intelligent Control and Automation (WCICA), pp. 2552-2555, DOI: 10.1109/WCICA.2016.7578828, .

[3] A. Sato, K. Shen, M. Minami, T. Matsuno, "Improvement of force-sensorless grinding accuracy with resistance compensation," Artificial Life and Robotics, vol. 22 no. 4, pp. 509-514, DOI: 10.1007/s10015-017-0385-y, 2017.

[4] K. Shin, "A study on the process optimization of brush deburring grinding system," Journal of Manufacturing Engineering & Technology, vol. 21 no. 3, pp. 394-400, DOI: 10.7735/ksmte.2012.21.3.394, 2012.

[5] S. W. Jeon, W. Jeong, D. Park, S.-B. Kwon, "Design of an intelligent duct cleaning robot with force compliant brush," Proceedings of the 2012 12th International Conference on Control, Automation and Systems, ICCAS 2012, pp. 2033-2037, .

[6] V. F. Makarov, A. V. Vinogradov, A. V. Nurtdinov, "Automated polishing of sharp edges on gas-turbine components by abrasive polymer brushes," Russian Engineering Research, vol. 32 no. 1, pp. 102-107, DOI: 10.3103/S1068798X12010182, 2012.

[7] V. Ganguly, T. Schmitz, A. A. Graziano, H. Yamaguchi, "An analysis of polishing forces in magnetic field assisted finishing," Proceedings of the ASME 2012 International Manufacturing Science and Engineering Conference,DOI: 10.1115/MSEC2012-7256, .

[8] R. S. Mulik, P. M. Pandey, "Experimental investigations into the finishing force and torque in magnetic abrasive finishing process," Proceedings of the ASME 2013 International Mechanical Engineering Congress and Exposition,DOI: 10.1115/IMECE2013-62365, .

[9] V. V. Shanbhag, K. Naveen, N. Balashanmugam, P. Vinod, "Modelling for evaluation of surface roughness in magnetic abrasive finishing of flat surfaces," International Journal of Precision Technology, vol. 6 no. 2,DOI: 10.1504/IJPTECH.2016.078190, 2016.

[10] S. R. Jefferies, "Abrasive finishing and polishing in restorative dentistry: a state-of-the-art review," Dental Clinics of North America, vol. 51 no. 2, pp. 379-397, DOI: 10.1016/j.cden.2006.12.002, 2007.

[11] A. Wiegand, J. P. M. Burkhard, F. Eggmann, T. Attin, "Brushing force of manual and sonic toothbrushes affects dental hard tissue abrasion," Clinical Oral Investigations, vol. 17 no. 3, pp. 815-822, DOI: 10.1007/s00784-012-0788-z, 2013.

[12] T. P. Ton, H. Y. Park, S. L. Ko, "Experimental analysis of deburring process on inclined exit surface by new deburring tool," CIRP Annals - Manufacturing Technology, vol. 60 no. 1, pp. 129-132, DOI: 10.1016/j.cirp.2011.03.124, 2011.

[13] N. Raymond, M. Soshi, "A study on the effect of abrasive filament tool on performance of sliding guideways for machine tools," Proceedings of the 3rd CIRP Conference on Surface Integrity (CIRP CSI), CIRP45, pp. 223-226, DOI: 10.1016/j.procir.2016.02.169, .

[14] G. Mathai, S. Melkote, "Effect of process parameters on the rate of abrasive assisted brush deburring of microgrooves," The International Journal of Machine Tools and Manufacture, vol. 57, pp. 46-54, DOI: 10.1016/j.ijmachtools.2012.02.007, 2012.

[15] M. A. Wahab, G. Parker, C. Wang, "Modelling rotary sweeping brushes and analyzing brush characteristic using finite element method," Finite Elements in Analysis and Design, vol. 43 no. 6-7, pp. 521-532, DOI: 10.1016/j.finel.2006.12.003, 2007.

[16] C. Wang, Q. Sun, M. A. Wahab, X. Zhang, L. Xu, "Regression modeling and prediction of road sweeping brush load characteristics from finite element analysis and experimental results," Waste Management, vol. 43, pp. 19-27, DOI: 10.1016/j.wasman.2015.06.027, 2015.

[17] M. M. Abdel-Wahab, C. Wang, L. V. Vanegas-Useche, G. A. Parker, "Experimental determination of optimum gutter brush parameters and road sweeping criteria for different types of waste," Waste Management, vol. 31 no. 6, pp. 1109-1120, DOI: 10.1016/j.wasman.2010.12.014, 2011.

[18] C. Wang, Q. Sun, L. Xu, "Analysis of helical abrasive brush filament force characteristics for surface treatment processes," Advances in Mechanical Engineering, vol. 9 no. 7,DOI: 10.1177/1687814017715980, 2017.

[19] D. G. Kleinbaum, L. L. Kupper, A. Nizam, E. S. Rosenberg, Applied Regression Analysis and Other Multivariable Methods, 2014.

[20] A. C. Rencher, W. F. Christensen, Methods of Multivariate Analysis,DOI: 10.1002/9781118391686, 2012.