A Novel Geometric Error Compensation Method for Gantry-Moving CNC Machine Regarding Dominant Errors

processes Article A Novel Geometric Error Compensation Method for Gantry-Moving CNC Machine Regarding Dominant Errors Hong Lu 1, Qian Cheng 1, Xinbao Zhang 2, *, Qi Liu 1, Yu Qiao 1 Yongquan Zhang 1 1 School of Mechanical Electronic Engineering, Wuhan University of Technology, Wuhan 430070, China; lzh@whut.edu.cn (H.L.); nhqian@whut.edu.cn (Q.C.); liunuli@whut.edu.cn (Q.L.); xiaoqiao456@whut.edu.cn (Y.Q.); zhangyongquan@whut.edu.cn (Y.Z.) 2 School of Mechanical Science Engineering, Huazhong University of Science Technology, Wuhan 430074, China * Correspondence: zhangxinbao1@hust.edu.cn Received: 7 July 2020; Accepted: 27 July 2020; Published: 31 July 2020 Abstract: Gantry-type computer numerical control (CNC) machines are widely used in manufacturing industry. A novel structure with moveable gantry is proposed to improve traditional gantry-type machine structure s disadvantage of taking up too much space. Geometric s have direct impacts on actual position of tool, which significantly reduces accuracy of machines. Errors of different components are always coupled have uncertain effects on total geometric. Thus, it is essential to find an effective way to identify dominant s do targeted. First, a novel identification method using value leaded global sensitivity analysis (VLGSA) is proposed to find dominant s. In VLGSA, weighting factors which show influence of range are used to modify multi-body system (MBS) model. Results show that dominant s in three directions respectively contribute 80%, 86% 85% of total in ir directions. n, s identified by VLGSA are modeled by least-square linear fitting B-spline interpolation methods, respectively, according to feature of data. Finally, models are applied in a new real-time system developed on Beckhoff TwinCAT servo system. Experimental results from gantry-type CNC engraving milling machine show proposed method can help figure out most dominant s reduce around 90% of total. Keywords: gantry-type CNC machine; geometric ; multi-body system; VLGSA; method 1. Introduction Gantry-type manufacturing equipment has already played an indispensable role in modern manufacturing industry. In traditional gantry-type machines, X-direction feed movement is realized by motion of workbench, while gantry is fixed on machine bed. A vast room of more than twice feed stroke should be reserved to ensure processing capability, which causes waste of space. Thus, a novel structure with a moveable gantry is proposed. Different from traditional structure, workbench is fixed two ball-screw feed systems are symmetrically arranged to drive gantry. Only half space is needed comparing with traditional one. However, machines with moveable gantries require higher accuracy. Many factors are affecting accuracy of gantry-type machine tools, such as geometric, rmal, servo system so forth. rmal s widely appear in engineering machinery devices many researchers have done great jobs on optimizing heat transformation [1 4]. Geometric is anor common component usually Processes 2020, 8, 906; doi:10.3390/pr8080906 www.mdpi.com/journal/processes

Processes 2020, 8, 906 2 of 21 donates most in final inaccuracy [5]. geometric s directly influence final position of tool nose s are generally related to tiny defect of machine itself, like wrong structure of feed system. re are two ways to resolve problem. One is adjusting physical structure, which is complicated expensive may be limited by skill of technical staff. or one is compensating for existing s by numerical control (NC) code or software, which is much more flexible can achieve high quality [6]. However, needs a clear understing of existing s. Error modeling is firstly required to reach goal. technique of geometric modeling has been persistently developed in recent decades. In application of robot industry, Paul, R.P. [7] Veitschegger, W.K. [8] gave methods using homogeneous transformation matrix differential transformation. precise location of robot s terminal point with existence of geometric size s kinematic s was found in ir research. se works established basis of analysis. Chen, J.X. et al. [5] took furr research based on fundamental principle of Paul Veitschegger, by analyzing differential movement differential transformation. y proposed a comprehensive method of both discovering how geometric propagates through every single axis. Zhao, Q.Q. et al. [9] established topology structure of machine tool using multi-body system ory proposed two conceptions, which were relative motion matrices of axis connectors connection matrix of machine tool connectors. Based on homogeneous coordinate transformation ory, Zhou, X.T. [10] built model of classic three-axis machine tool, decomposed final into s existing in components in whole motion chain. se works showed development application of different modeling methods, gave inspirations to furr research. Still, many of m were limited in modeling. Models based on se methods are usually complex mamatical expressions of a vast amount of items, making it a sophisticated work to compensate for all items. However, towards several specific items, which can be defined as dominant s, may provide enough accuracy improvement. Thus, more research could be done to simplify. Sensitivity analysis (SA) is a kind of method to quantify uncertainties of an unknown system using great data, which may help in finding most valuable items to compensate. For a specific model, a vast amount of different input output data can be obtained through all kinds of scientific sampling methods powerful calculation tools. Thus sensitivity analysis can help to identify weight each input has in affecting system outputs. SA consists of two main branches, one is local sensitivity analysis (LSA) or is called global sensitivity analysis (GSA). [11,12] GSA is more commonly used in manufacturing industry with its ability to find out result differences under changes of parameters. It can emphasize interaction power of parameters by using additional methods. Usually, SA is used in economic analysis. However, it is also found suitable in engineering field recently. Fan, J.W. et al. [13] established a prediction model of component tolerance proposed a sensitivity analysis based on component motion. proposed sensitivity analysis was experimentally verified several main parameters that donate to s most were found out. Yan, Y. et al. [14] established model of multi-process assembling built sensitivity analysis model of parameters towards assembling using method of matrix differentiation. Li, J. Xue, F.G. [15] proposed general local sensitivity indices (GLSIs), general global sensitivity indices (GGSIs) general global sensitivity fluctuation indices (GGSFIs) based on traditional sensitivity analysis. A five-axis machine was analyzed using above indices. One common limit of above works is that sensitivity indices calculation is carried out on original model without considering impact of value range towards sensitivity analysis. Besides, sensitivity analysis of machine tools was only done for performance assessment, targeted was missing validity was not verified. Sensitivity analysis tells what to concentrate on; measurement fitting tell how to describe deal with m. With all kinds of advanced measurement equipment, precision data can be obtained. According to data, model can be established for. However, appropriate fitting methods are needed to ensure accuracy. To improve

Processes 2020, 8, 906 3 of 21 manufacturing performance, many researchers have done great works. Fang, K.G. Yang, J.G. [16] considered rmal factors, used orthogonal polynomial method to build total model, compensated 90% of spindle geometric. Wang, W. et al. [17] used Newton interpolation method to construct rmal geometric model of computer numerical control (CNC) milling machine, established a real-time system. experiment based on Fanuc CNC machine was carried out to show system works well. Tuo, Z.Y. et al. [18] using cubic spline interpolation established rmal geometric model improved position accuracy of CNC machine. B-spline curve is commonly used in computer aided design (CAD) areas for its ability to describe continual material. It is also used in trajectory planning for its smooth character. Li, Y.H. et al. [19] used quantic B-spline curve to optimize motion path of pick--place robots. His work showed that B-spline curve is appropriate for smooth machine motion, which is also suitable for relative movement of machine tools. In this paper, multi-body system ory was used to establish geometric model of gantry-moving engraving milling machine. n translation s in all three directions are extracted from matrix to analyze significant items which impact precision most. Value leaded global sensitivity analysis (VLGSA) is proposed used in above analysis process, in which measurement value of each item is considered to modify total geometric model. B-spline curve was developed into B-spline interpolation method used on result of machine measurement. A real-time system based on Beckhoff servo system, TwinCAT control software C# language was built to eliminate s. experiment verified validity of above crucial finding method. 2. Geometric Modeling Based on Multibody System multi-body system (MBS) illustrates a kind of mechanical system in which a series of rigid or flexible bodies interconnect [20]; every single body may undergo translation rotation in different scales. It is a general systematic method for geometric modeling is widely used in analysis of robots, machine tools coordinate measurement. core of MBS is using low-order-body arrays to describe topology structure of system [20], using transformation matrices to calculate motion relationship. Usually, inertial coordinate system is chosen as body B 0, basic component is selected as body B 1, or elements are divided into several branches numbered along direction away from B 1, branch by branch. 2.1. Coordinate System Definition Topology Structure Description of Gantry-Moving CNC Machine In this paper, a gantry-moving engraving milling machine with three-dimension movements is taken to describe ideal of presented geometric modeling method. structure of chosen machine Processes 2020, is shown 8, x FOR PEER in Figure REVIEW 1. 4 of 24 Figure 1. Figure Structure 1. Structure of of gantry-moving engraving milling milling machine. machine. B0-bed, B1-moveable B gantry, 0 -bed, B 1 -moveable B2-Y B3-Z direction slide, B4-spindle with tool, B5-workbench, B6-workpiece. T01 gantry, B 2 -Y direction slide, B 3 -Z direction slide, B 4 -spindle with tool, B 5 -workbench, B 6 -workpiece. T56 represent coordinate transformation between components. T 01 T 56 represent coordinate transformation between components. gantry-moving CNC machine has two servo motors in X direction. Two motors are symmetrically arranged to drive gantry through ball-screw feed systems. X1 X2 represent movement of two feed systems. synchronization between two movements is an important cause of X-direction yaw. Since machine locates on a stable ground coordinate system, machine bed itself can be directly chosen as body B0. As orange lines show in Figure 1, engraving mill machine can be seen as a two-branch topology structure [21]. first branch (tool branch) consists of

Processes 2020, 8, 906 4 of 21 gantry-moving CNC machine has two servo motors in X direction. Two motors are symmetrically arranged to drive gantry through ball-screw feed systems. X1 X2 represent movement of two feed systems. synchronization between two movements is an important cause of X-direction yaw. Since machine locates on a stable ground coordinate system, machine bed itself can be directly chosen as body B 0. As orange lines show in Figure 1, engraving mill machine can be seen as a two-branch topology structure [21]. first branch (tool branch) consists of moveable gantry (B 1 ), Y direction slide (B 2 ), Z direction slide (B 3 ) spindle with tool on it (B 4 ). second branch (workpiece branch) consists of workbench (B 5 ) workpiece (B 6 ). Table 1 shows main sizes of se components. Table 1. Size of main components. Main Parameters of Components Value (mm) Size of workbench (width length) 3000 1500 Stroke of gantry (X-direction stroke) 2400 Stroke of Y slide (Y-direction stroke) 1200 Stroke of Z slide (Z-direction stroke) 240 final topology structure coordinate definition with position motion relationships are shown in Figure 2, where vectors between bodies show ir relations. Take B 0 B 1 as an example, T 01p means position transformation between B 0 B 1 ; T 01s means motion transformation between B 0 B 1 ; T 01pe T 01se represent of position motion transformations respectively. T 01p, T 01pe, T 01s T 01se constitute total coordinate transformation betweenprocesses B 0 2020, B8, x 1 (TFOR 01 ). PEER REVIEW 5 of 24 Figure 2. Figure 2. topology topology structure coordinate definition of of machine. Orange Orange lines show lines show coordinate coordinate transformation transformation between between different bodies, which are are marked with with T01p, T01p, T01pe, T01pe, etc. etc. Table 2 shows low-order-body array of machine, in which L is low-order-body Table 2 shows low-order-body array of machine, in which L is low-order-body operator. operator. For body Bj, its n-th order low-order-body number can be defined as Equation (1): For body B j, its n-th order low-order-body number can be defined as Equation (1): n L(j) = i, (1) where n means that Bi is n-th order low-order-body L n (j) = i of body Bj (1) Similarly, or equation can be inferred as Equations (2) (3): where n means that B i is n-th order low-order-body n n 1 of body B L(j) = L(L (j)) j (2) 0 L(j) = j. (3) Table 2. Low-order-body Array of Chosen Machine. Body 1 2 3 4 5 6 L1(j) 0 1 2 3 0 5

Processes 2020, 8, 906 5 of 21 Similarly, or equation can be inferred as Equations (2) (3): L n (j) = L(L n 1 (j)) (2) L 0 (j) = j. (3) Table 2. Low-order-body Array of Chosen Machine. Body 1 2 3 4 5 6 L1(j) 0 1 2 3 0 5 L2(j) 0 1 2 0 L3(j) 0 1 L4(j) 0 2.2. Coordinate Transformation Model Establishment actual position motion between bodies can be expressed by matrix, more specifically by position matrix, motion matrix, position matrix motion matrix. After binging every single body to ir fixed coordinate, position motion transformation can be demonstrated by multiplication of a series of 4 4 homogeneous eigenmatrix [22]. For a three-dimensional Cartesian coordinate system, translation can be expressed as Equation (4), in which x, y z are translation along three axes, respectively. Rotation around three axes can be expressed as Equations (5) (7). 1 0 0 x 0 1 0 y T = (4) 0 0 1 z Rot(x, θ x ) = Rot(y, θ y ) = Rot(z, θ z ) = 1 0 0 0 0 cosθ x sin θ x 0 0 sin θ x cos θ x 0 cos θ y 0 sin θ y 0 0 1 0 0 sin θ y 0 cos θ y 0 cos θ z sin θ z 0 0 sin θ z cos θ z 0 0 0 0 1 0 However, both structure motion exist in manufacturing process. Take X-direction movement as an example (shown in Figure 3). Errors between two coordinates can be divided into six types, including three linear three angle. Linear s are position straightness along with two directions. Angle s are roll, pitch yaw. Similarly, geometric s in or directions can be defined as Table 3. (5) (6) (7)

sinθ cosθ 0 0 z z Rot(z, θ ) = z 0 0 1 0 However, both structure motion exist in manufacturing process. Take X- direction movement as an example (shown in Figure 3). Errors between two coordinates can be divided into six types, including three linear three angle. Linear s are position Processes 2020, 8, 906 6 of 21 straightness along with two directions. Angle s are roll, pitch yaw. Similarly, geometric s in or directions can be defined as Table 3. (7) (a) (b) (c) (d) Figure 3. Error measurement. (a) Principle of position measurement; (b) Principle of straightness measurement; (c) Position measurement experiment; (d) Straightness measurement experiment. Figure 3. Error measurement. (a) Principle of position measurement; (b) Principle of straightness measurement; (c) Position measurement experiment; (d) Straightness measurement experiment. Table 3. Geometric s in chosen caving milling machine. Error Type Error Item Error Symbol Position x x Y-direction straightness y x X-direction feed s Z-direction straightness z x Roll α x Pitch β x Yaw γ x Position y y X-direction straightness x y Y-direction feed s Z-direction straightness z y Roll β y Pitch α y Yaw γ y Position z z X-direction straightness x z Z-direction feed s Y-direction straightness y z Roll γ z Pitch β z Yaw α z X-Y perpendicularity γ xy Structure s X-Z perpendicularity β xz Y-Z perpendicularity α yz According to principle of homogeneous matrix transformation, general matrix between B 0 B 1 which containing all six s can be calculated by multiplying six homogeneous eigenmatrices, as Equation (8) shows (where c for cos s for sin ). T 01e = T x x T y x T z x T θx x T θy x T θz x c θy x c θz x c θy x s θz x s θy x x x c θx x s θz x + s θx x s θy = x c θz x c θx x c θz x s θx x s θy x s θz x c θy x s θx x y x s θx x s θz x c θx x c θz x s θy x c θz x s θx x + c θx x s θy x s θz x c θx x c θy x z x (8)

Processes 2020, 8, 906 7 of 21 In translation movement, three angle s are usually tiny, so approximation that sin θ θ cos θ 1 can be introduced. Putting approximation into Equation (8) ignoring second-order high-order elements, equation can be simplified as Equation (9). T 01e = 1 θz x θy x x θz x 1 θx x y θy x θx x 1 z According to method above, transformation matrices of chosen engraving mill machine are listed in Table 4. Table 4. transformation matrix of chosen caring mill machine. Neighbor Bodies Ideal Position Motion Matrix Position Motion Error Matrix 1 γ xy β xz 0 γ T 01p = I 4 4 T 0 1 01pe = xy 1 α yz 0 β xz α yz 1 0 1 0 0 x 1 γ x β x x x 0 1 0 0 γ T 01s = T 0 0 1 0 01se = x 1 α x y x β x α x 1 z x (9) 1 2 2 3 3 4 0 5 5 6 T 12s = T 23s = T 12p = I 4 4 T 12pe = I 4 4 1 0 0 0 0 1 0 y T 0 0 1 0 12se = 1 γ y β y x y γ y 1 α y y y β y α y 1 z y T 23p = I 4 4 T 23p = I 4 4 1 0 0 0 0 1 0 y T 0 0 1 0 23se = 1 γ z β z x z γ z 1 α z y z β z α z 1 z z T 34p = I 4 4 T 34pe = I 4 4 T 34s = I 4 4 T 34se = I 4 4 T 05p = I 4 4 T 05pe = I 4 4 T 05s = I 4 4 T 05se = I 4 4 T 56p = I 4 4 T 56pe = I 4 4 T 56s = I 4 4 T 56se = I 4 4 coordinate system eigenmatrix of two neighboring bodies can be written as Equation (10). Eigenmatrix with s can be written as Equation (11). T 01 = T 01p T 01pe T 01s T 01se (10) T 01e = T 01 T 01i. (11) Coordinate system eigenmatrix of any two bodies can be written as Equation (12). T SG = K=1 K=m,L m (S)=Q T K L (S)LK 1 (S) 1 J=1 J=n,L n (G)=Q T J L (G)LJ 1 (G) (12) As for topology structure of chosen engraving milling machine mentioned above, coordinate transformation eigenmatrix in both ideal situation actual situation (with s) should be calculated as Equation (13) T 64 = ( k=1 1 j=1 T L k (6)L (6)) k 1 T L j (4)L j 1 (4) k=2 j=4 (13) = (T 05 T 56 ) 1 T 01 T 12 T 23 T 34

Processes 2020, 8, 906 8 of 21 geometric of machine can n be calculated by subtracting ideal transformation eigenmatrix actual transformation eigenmatrix between tool workpiece, like Equation (14). result T 64e is final matrix containing angle s linear s between tool coordinate system workpiece coordinate system. In planar milling, angle s represent small posture fault linear s along three axes have most significant impact on misalignment between tool nose target point on workpiece. 3. Dominant Errors Identification with VLGSA T 64e = T 64 T 64i (14) From perspective of black box, system models can always be demonstrated as a function Y = f(x), in which x is vector of several uncertain inputs of system y is one of system outputs which is chosen to be analyzed. In this paper, model targets total geometric of engraving milling machine. Y represents final X equals vector of each item like X-direction position. By doing VLGSA to model function, effects of inputs on a particular output can be quantified into indices. In or words, impact of each item on total geometric can be described with numbers dominant s, which cause most of total, can be found. 3.1. Decomposition of Variance Indices Calculating In system, model function Y = f(x) is under assumption that inputs distribute independently within a unit uniform space. refore, actual input scales should be transformed into unit space [23]. function can be decomposed as Equation (15). Y = f 0 + d f i (X i ) + i=1 1 i<j d f ij (X i, X j ) +... + f 12...d (X 1, X 2,..., X d ), (15) where f 0 is a component of constant; f i is function of X i, which represents effect X i donates alone ( geometric caused by one specific component); f ij is function with two inputs X i X j, representing effect to output given by X i X j toger. Similarly, or higher parts can be defined as above. This composition must be in condition that sub-components are independent of each or. In or words, all components in equation are orthogonal. It can be shown in following Equation (16). 1 f i1 i 2...i s (X i1, X i2,..., X is )dx ik = 0, (k = 1, 2,..., s). (16) 0 It makes result that components of decomposition function can be defined in terms of expected values, as Equations (17) (19) show. f 0 = E(Y) (17) f i (X i ) = E(Y Xi ) f 0 (18) f ij (X i, X j ) = E(Y Xi, X j ) f 0 f i f j (19) Or components containing more inputs can be defined as equations above. f i which shows effect of change input X i can be known as primary interaction or first-order interaction. f ij shows effect of changing both X i X j at same time. It can be known as second-order interaction. higher-order interactions can be defined similarly.

Processes 2020, 8, 906 9 of 21 Square integrate decomposition Equation (15), new equation is got can be organized as Equations (20) (21). 1 0 f2 (X)dX = f 0 2 + d 1 0 i=1 1 0 f i 2 dx i + 1 i<j d f 2 (X)dX f 0 2 = 1 1 0 0 f ij 2 dx i dx j +... + 1 0... 1 0 f 12...d 2 dx 1 dx 2... dx d (20) d j=1 1 i<j d 1 0 1... f 2 i...j dx i... dx j (21) 0 It is found that left part of equation is variance of Y(X) right part is variance terms. This leads to decomposition with variance expression, shown as Equations (22) (24). Var(Y) = d V i + i=1 1 i<j d V ij +... + V 12...d (22) V i = Var Xi (E X i (Y Xi )) (23) V ij = Var Xij (E X ij (Y X i, X j )) V i V j, (24) where X i represents whole series of variables except X i. n first-order sensitivity index, which is direct sensitivity factor based on variance, can be defined as Equation (25). S i = V i Var(Y). (25) After estimating, component Var Xi (E X i (Y Xi )) can be written as Equation (26). Var Xi (E X i (Y X i )) 1 N N f(b) j (f(ab i ) j f(a) j ) (26) j=1 second-order higher-order sensitivity index can be defined in same way, which represents effect of output variance causing by change of several inputs, simultaneously. sensitivity indices satisfy following principle. d S i + i=1 1 i<j d S ij +... + S 12...d = 1 (27) For some functions which are simple in structure, first-order sensitivity indices can be obtained by calculating integrals of decomposition function. However, more common way for most cases is estimating, which provides an easy solution to all kinds of functions it is often finished by Monte Carlo method. Monte Carlo method needs to generate a set of rom points in unit hypercube space, n use points to calculate use statistic method to estimate character of model. In actual work, rom points are usually pseudorom. most used points sequence is Latin hypercube sequence Sobol sequence. Several steps should be taken as follow to calculate sensitivity index using estimation method. 1. Generate a matrix of N 2d, in which N is number of sample points d is dimension of model (number of inputs). Points in se matrices should be romly distributed. 2. Generate matrix A using first d columns of above matrix matrix B using rest d columns. se two matrices distribute first two samples, which have n points in d-dimension hypercube space.

Processes 2020, 8, 906 10 of 21 3. Generate d more matrices AB i by replacing i-th column of matrix A with same column in matrix B. 4. Use above d + 2 matrices as inputs of model get d + 2 result matrices. ses matrices contain f(a) j, f(b) j f(ab i ) j mentioned in Equation (26). 5. Use matrices above to calculate sensitivity index with Equation (27). 3.2. Machine Geometric Error Measurement for Value Leaded Factors Construction sample points of global sensitivity analysis distribute in unit hypercube space. However, items are different in scale usual sample method may cause inaccuracy in estimating actual impact of each item. Thus, a value-leaded factor based on measurement is proposed to modify geometric model, which normalizes s effects to same level. Keysight 5530 dual-frequency laser interferometer measurement system is used in measurement, cooperating with several kinds of mirror groups, 21 geometric s of engraving milling machine can be detected. Measurements are done with stard of ISO230-1 processes are shown in Figure 3. This work includes many time-consuming steps. For each component, measurement is done in different feed speeds, with range from 500 mm/min to 6000 mm/min interval of 500 mm/min. measurement stroke is x = 2000 mm, y = 1000 mm, z = 200 mm, which covers main processing area. start points are set in x = 200 mm, y = 100 mm z = 20 mm in machine coordinate system to leave enough space for installing measurement instruments. All measurements are done for three times. scale of each ir ratios are obtained from measurements. y are used for normalization of sample of global sensitivity analysis. measured s are shown in Table 5. Table 5. Max value in whole stroke. No. Item Max Value No. Item Max Value 1 x x 0.1212 mm 12 γ y 0.0034 mm/1000 mm 2 y x 0.0456 mm 13 z z 0.0218 mm 3 z x 0.0805 mm 14 x z 0.0213 mm 4 α x 0.0037 mm/1000 mm 15 y z 0.0337 mm 5 β x 0.0032 mm/1000 mm 16 γ z 0.0033 mm/1000 mm 6 γ x 0.0034 mm/1000 mm 17 β z 0.0035 mm/1000 mm 7 y y 0.0728 mm 18 α z 0.0036 mm/1000 mm 8 x y 0.1068 mm 19 γ xy 0.0125 mm/1000 mm 9 z y 0.1538 mm 20 β xz 0.0082 mm/1000 mm 10 β y 0.0036 mm/1000 mm 21 α yz 0.0063 mm/1000 mm 11 α y 0.0037 mm/1000 mm 0.0037 mm/1000 mm means deviation caused by α x in stroke of 1000 mm is 0.0037 mm. Or angle s are described in same way. 3.3. Indices Calculation Global Sensitivity Analysis analysis program is made in Matlab according to steps of proposed VLGAS method. geometric model of engraving milling machine built with multi-body system ory is adjusted for normalization based on measurement. Putting model into analysis program setting analysis parameters, first-order sensitivity indices in three directions are calculated listed in Table 6. angle s are at a very low level, high-order elements linking to angle s can be eliminated from MBS model. Thus, some elements sensitivity indices are zero in a particular direction. For X-direction movement, as shown in Figure 4, with geometric s changing in a certain range, sensitivity indices of x x, x y x z are largest contribute 80% of total value toger. For Y direction, most dominant element is y y, which has highest sensitivity index value of 0.424, or two dominant elements are y x y z, se three elements contribute about

Processes 2020, 8, 906 11 of 21 86% of total value. Similarly, three linear s z x, z y z z have largest impact on Z direction, with proportion of around 85%. Number Table 6. First-order sensitivity index of each item. Item First-Order Sensitivity X-Direction Y-Direction Z-Direction 1 x x 0.333 0.002 0.010 2 y x 0.004 0.277 0.011 3 z x 0.006 0.003 0.092 4 α x 0.000 0.018 0.014 5 β x 0.029 0.003 0.007 6 γ x 0.039 0.007 0.000 7 y y 0.003 0.424 0.007 8 x y 0.306 0.027 0.008 9 z y 0.010 0.000 0.460 10 β y 0.028 0.000 0.007 11 α y 0.000 0.018 0.008 12 γ y 0.005 0.005 0.000 13 z z 0.004 0.004 0.297 14 x z 0.160 0.001 0.003 15 y z 0.006 0.160 0.007 16 γ z 0.000 0.000 0.000 17 β z 0.000 0.000 0.000 18 α z 0.000 0.000 0.000 19 γ xy 0.002 0.058 0.000 20 β xz 0.027 0.000 0.058 21 α yz 0.000 0.018 0.014 data is drawn shown in Figure 4 to demonstrate clearly. In conclusion, for chosen engraving milling machine, dominant s are mainly Processes 2020, 8, x FOR PEER REVIEW 14 of 24 linear items, specifically speaking position two straightness s from or two feed directions. feed directions. In orin words, or words, aiming at at se se three three s s can provide can provide remarkable remarkable improvement improvement to chosen to chosen machine machine with with minimum cost. (a) (b) (c) Figure 4. Figure Image4. Image of first-order of sensitivity index. index. (a) First-order (a) First-order sensitivity of sensitivity X-direction ; of (b) X-direction Firstorder sensitivity of Y-direction ; (c) First-order sensitivity of Z-direction. ; (b) First-order sensitivity of Y-direction ; (c) First-order sensitivity of Z-direction. 4. Error Compensation In this section, data in whole stroke are listed modeled based on ir character. Different models are applied to a novel real-time system compared. linear in X direction is chosen as an example. Errors in or directions can be dealt same way. As shown in last section, main components, which have largest impact on X-direction linear, are X-direction position, X-direction straightness along with Y- direction movement X-direction straightness along with Z-direction movement. y all have a close relation to feed movement, which means y can be compensated by

Processes 2020, 8, 906 12 of 21 4. Error Compensation In this section, data in whole stroke are listed modeled based on ir character. Different models are applied to a novel real-time system compared. linear in X direction is chosen as an example. Errors in or directions can be dealt same way. As shown in last section, main components, which have largest impact on X-direction linear, are X-direction position, X-direction straightness along with Y-direction movement X-direction straightness along with Z-direction movement. y all have a close relation to feed movement, which means y can be compensated by adjusting feed trajectory. 4.1. Position Error Model Based on Least-Square Method Position s are measured three times on each point average values are calculated to reduce measurement. values are listed in Tables 7 9. Position (mm) Table 7. Position in X direction. Error (µm) First Second Third Average 0 0 0.012 0.030 0.014 200 4.797 6.652 6.051 5.833 400 21.881 19.252 18.362 19.831 600 32.453 31.848 29.204 31.168 800 45.158 40.832 41.461 42.483 1000 59.489 59.900 59.560 59.649 1200 77.401 68.519 63.663 69.861 1400 85.961 80.098 78.703 81.587 1600 95.536 90.457 87.197 91.063 1800 111.129 106.270 98.881 105.427 2000 126.542 118.889 118.035 121.155 Origin of position located on x = 200 mm in machine coordinate system. Position (mm) Table 8. Position in Y direction. Error (µm) First Second Third Average 0 0.015 0.035 0.017 0.022 100 14.051 14.713 17.019 15.261 200 12.107 11.935 15.522 13.188 300 17.304 17.404 22.221 18.976 400 25.420 25.753 30.710 27.294 500 21.872 22.118 27.867 23.952 600 37.037 37.911 43.779 39.575 700 43.048 43.633 49.924 45.535 800 52.870 54.024 61.151 56.015 900 66.971 68.642 72.338 69.317 1000 71.151 72.285 74.859 72.765 Origin of position located on y = 100 mm in machine coordinate system. According to tables, position s show an obverse trend of linear distribution in whole stroke. value accumulates through feed movement maximum value appears at end of stroke. maximum values in three directions are 121.155, 72.765 19.215, respectively. Thus, least square method is chosen to establish mamatical model of position s to describe linear trend.

Processes 2020, 8, 906 13 of 21 Position (mm) Table 9. Position in Z direction. Error (µm) First Second Third Average 0 0.001 0.053 0.026 0.0267 20 2.384 2.422 2.032 2.279 40 4.349 4.725 4.561 4.545 60 6.437 6.542 6.443 6.474 80 8.771 8.713 8.739 8.741 100 10.944 11.021 10.867 10.944 120 13.196 13.022 12.879 13.032 140 15.562 15.469 15.264 15.432 160 17.372 17.551 17.411 17.445 180 19.198 19.131 19.316 19.215 200 21.542 22.037 21.693 21.757 Processes 2020, 8, x FOR PEER REVIEW 16 of 24 Origin of position located on z = 20 mm in machine coordinate system. end of stroke. maximum values in three directions are 121.155, 72.765 19.215, respectively. Thus, least square method is chosen to establish mamatical model of least-square method is a commonly used fitting method, in which fitting curve not position s to describe linear trend. necessarily cross least-square sample points method but is a makes commonly se used points fitting as method, closein aswhich possible fitting to curve curve. not For sample points necessarily (x i, y i ), cross a series sample of equations points but can makes bese built points as used close toas create possible to function curve. For of points, sample points ( x residual sum of squares n i, y i ), a series δ 2 of equations can be built used to create function of points, is chosen to decide best coefficients of unknown function. i residual sum of squares i=1 n 2 δ is chosen to decide best coefficients of unknown function. i i= 1 fitting functions are shown as Equations (28) (30). fitting curves of position s are fitting functions are shown as Equations (28) (30). fitting curves of position s are shown in Figure 5. shown in Figure 5. y = 0.0612x + 4.09 (28) y = 0.0612x + 4.09 (28) y = 0.071x 0.117 (29) y = 0.071x 0.117 (29) y = 0.108x y = 0.0958 (30) (30) (a) (b) (c) Figure 5. Image Figure of 5. fitting Image of curves. fitting curves. (a) Fitting (a) Fitting curve of X-direction position position ; ; (b) Fitting (b) curve Fitting of Y-curve of direction position ; (c) Fitting curve of Z-direction position. Y-direction position ; (c) Fitting curve of Z-direction position. 4.2. Straightness Error Model Based on B-Spline Interpolation. straightness s are also measured three times to calculate average value. curve of s is drawn in Figure 6 to show intuitively.

Processes 2020, 8, 906 14 of 21 4.2. Straightness Error Model Based on B-Spline Interpolation. straightness s are also measured three times to calculate average value. curve of Processes 2020, 8, x FOR PEER REVIEW 17 of 24 s is drawn in Figure 6 to show intuitively. (a) (b) Figure 6. Figure Image 6. of Image straightness of straightness. (a). Straightness (a) Straightness in X-direction in X-direction movement; movement; (b) Straightness (b) Straightness in Y-direction movement; (c) Straightness in Z-direction movement. in Y-direction movement; (c) Straightness in Z-direction movement. According to Figure 6, max value of straightness s is lower than position s. According to Figure 6, max value of straightness s is lower than position s. curves fluctuate with total trend of increasing first n decreasing, except Y-direction curves fluctuate along with with X-direction totalmovement. trend of increasing maximum first values respectively n decreasing, appear in except point Y-direction of along x1 = 1600 withmm, x2 X-direction = 1000 mm, movement. y1 = 400 mm, y2 = 500 maximum mm, z1 = values 100 mm respectively z2 = 100 mm appear for insix point of x 1 = straightness 1600 mm, xs. 2 = 1000 mm, distortion y 1 = 400 of mm, guides y 2 = during 500 mm, manufacturing z 1 = 100 mm is leading z 2 = 100 cause mmof for six straightness s. distortion of guides during manufacturing is leading cause of B-spline curve is a kind of unique curve, which is developed from Bezier curve is straightness s. constructed by linear combination of B-spline base curves. B-spline curves have many excellent properties, B-spline such curve as geometric is a kind invariance, of unique convex curve, hull, which convexity is developed retention, from variation reduction Bezier curve is constructed local support by linear so forth. combination B-spline surface of B-spline reconstruction base curves. is one of B-spline research curves hotspots have many critical excellent properties, technologies such asof geometric reverse engineering. invariance, refore, convexin hull, this paper, convexity B-spline retention, interpolation variationis reduction used to local support establish so mamatical forth. B-spline surface model of reconstruction straightness is, onewhich of can research describe hotspots continual critical technologies material of property reverse of engineering. guide. B-spline refore, base in function this paper, curve construction B-spline interpolation function are listed is used to in Appendix A. establish mamatical model of straightness, which can describe continual Exping measurement point number doing simulation with material interpolation property of result, guide. it is shown B-spline that values base function of residual curve keep construction at a low level function after are listed in Appendix interpolation A.. Interpolation results are shown in Figure 7. Exping measurement point number doing simulation with interpolation result, it is shown that values of residual keep at a low level after interpolation. Interpolation results are shown in Figure 7. (c)

Processes 2020, 8, 906 15 of 21 Processes 2020, 8, x FOR PEER REVIEW 18 of 24 (a) (b) (c) (d) (e) (f) Figure Figure 7. Compensation 7. Compensation simulation simulation with with B-spline B-spline interpolation model. model. (a) (a) Y-direction straightness in X-direction in X-direction movement; movement; (b) Z-direction (b) Z-direction straightness straightness in X-direction in X-direction movement; movement; (c) X-direction (c) direction straightness in Y-direction movement; (d) Z-direction straightness in Y-direction straightness in Y-direction movement; (d) Z-direction straightness in Y-direction movement; movement; (e) X-direction straightness in Z-direction movement; (f) Y-direction straightness (e) X-direction straightness in Z-direction movement; (f) Y-direction straightness in Z-direction movement. in Z-direction movement. 4.3. Compensation 4.3. Compensation System System Experiment Experiment Aiming to, many researchers in previous works try to change G code Aiming to, many researchers in previous works try to change G code according to modeling result, which performs well in reducing value [24]. according However, to this kind of modeling method need result, pre-work which for performs each program well in cannot reducing be real-time. With value [24]. However, advantages this kind of Beckhoff s of method open need CNC pre-work system, for eachmodel program can help to cannot improve bemachine real-time. accuracy With advantages in a more of Beckhoff s effective open real-time CNCway. system, can help to improve machine accuracy in a more effective real-time way. strategy combines virtual axis electronic cam technology based on TwinCAT control system of Beckhoff. core of strategy is uploading model points into cam table to help change actual tool path offsetting trajectory deviation. Firstly,

Processes 2020, 8, x FOR PEER REVIEW 19 of 24 strategy combines virtual axis electronic cam technology based on Processes 2020, 8, 906 16 of 21 TwinCAT control system of Beckhoff. core of strategy is uploading model points into cam table to help change actual tool path offsetting trajectory deviation. Firstly, a virtual axis is established to accept motion comm. virtual axis is a kind of electronic symbol in control software. It It does not connect to any physical component only has digital motion simulation. refore, it itrepresents represents ideal ideal motion motion state. state. Secondly, Secondly, set set virtual virtual axis axis as as master master create create anor anor axis axis symbol symbol connecting connecting to to physical physical axis axis motor motor as as slave. slave. Thus, Thus, slave slave axis axis can can move move with with master master axis, axis, which whichis is virtual virtualaxis axis has hasno physical output, according to to specific specific rules. rules. Thirdly, Thirdly, set set rule rule between between two axes two as axes electronic as electronic cam, n cam, n physical axis physical can do axis superimposed can do superimposed movements movements of both ideal of both movements ideal movements cam-curve cam-curve movements movements based on based on models. models. cam cam table table can can be created be created in in two two ways: ways: cam cam design design tool tool assembled in TwinCAT software software user user interface interface linking linking to softto programmable soft programmable logic controller logic controller (PLC) program. (PLC) Soft program. PLC method Soft PLC has method better real-time has better property real-time property table value can table bevalue changed can in be changed feed process. in feed process. system can be described system in can Figure be described 8. hardware in Figure of 8. hardware system of is shown in Figure system 9. is shown in Figure 9. Figure 8. structure of system. Figure 9. hardware of system. Suppose that x is ideal position when feed system moves along X direction xδx is is position generated in in feed feed stroke, stroke, n n cam cam table table value value in this in this position position should should be set be as set x. as Δx Give. Give motion motion comm comm to to virtual virtual axis, axis, electronic electronic cam cam offset offset will will compensate compensate in position x. With help of Beckhoff ErCAT communication system, massive data can be transmitted in real-time. refore, coordinate information in all directions can be read cam

Processes 2020, 8, x FOR PEER REVIEW 20 of 24 Processes 2020, 8, x FOR PEER REVIEW 20 of 24 Processes 2020, 906 17 data of 21 in 8,position x. With help of Beckhoff ErCAT communication system, massive can transmitted refore, coordinate information in all directions be read be in positioninx.real-time. With help of Beckhoff ErCAT communication system,can massive data cam table can be modified with calculated value based on models, can can be transmitted real-time. refore, coordinate information in all can be read table be modifiedinwith calculated value based on directions models, without delay. without process expressed in Figure 10.value based on models, camdelay. table can be modified with is calculated process is expressed in Figure 10. without delay. process is expressed in Figure 10. Figure 10. process of proposed method. Figure processof of proposed proposed method. Figure 10. 10. process method. system is composed of both hardware software. hardware platform composedof ofboth bothhardware hardware software. software. hardware platform system system is is composed hardware platform of system is provided by Beckhoff industrial personal computer (PC) ofof by Beckhoff Beckhoffindustrial industrialpersonal personal computer (PC) system system is provided provided by computer (PC) software is written in C# language. With Beckhoff ADS communication protocol, C# language. language.with With Beckhoff Beckhoff ADS communication protocol, software software is is written written in C# ADS communication protocol, data in PLC Servo controller areare easy toto bebe dealt with algorithm inin software data PLC Servo controller easy dealt with algorithm software data inin PLC Servo controller are easy to be dealt with algorithm in software n transmitted back. n ntransmitted transmittedback. back. Straightness is measured verify effect. Straightness after after is effect. Straightness after ismeasured measuredtoto toverify verify effect. models built by three or usual fitting methods (Polynomial Fitting, Newton Interpolation Cubic models built by three or usual fitting methods (Polynomial Fitting, Newton Interpolation models built by three or usual fitting methods (Polynomial Fitting, Newton Interpolation Spline are usedare to compare with with chosen method ofmethod B-spline interpolation. As for CubicInterpolation) SplineInterpolation) Interpolation) are used to chosen of of B-spline interpolation. Cubic Spline used to compare compare with chosenmethod B-spline interpolation. total effect, a space movement of x = 2000 mm, y = 1000 mm z = 200 mm carried for total total effect, effect, aa space mm, y =y 1000 mm z = z200 mmmm AsAsfor spacemovement movementofofx x= =2000 2000 mm, = 1000 mm =is200 out X-direction total is measured. values are measured by absolute grating is carried out X-direction total is measured. values are measured by absolute is carried out X-direction total is measured. values are measured by absolute grating rulers. total isatcurbed at a low level, both modeling method rulers. Ifrulers. total is curbed a low level, both method sensitivity grating IfIf total is curbed at a low level, bothmodeling modeling method sensitivity analysis method can be verified to be effective, for only dominant components analysis method can be verified to be effective, for only dominant components analyzed by sensitivity analysis method can be verified to be effective, for only dominant components analyzed by method are compensated. Compensation resultsin are shown11. in Figure 11. method are compensated. Compensation results are shown Figure analyzed by method are compensated. Compensation results are shown in Figure 11. (a) (a) (b) Figure 11. Cont. (b)

Processes Processes 2020, 2020, 8, 8, 906 x FOR PEER REVIEW 21 of 1824 of 21 (c) Figure 11. 11. Compensation experiment result. (a) Residual Y-direction straightness in in X-direction movement before after ; (b) Residual Z-direction straightness in in X-direction movement before after ; (c) (c) Residual Residual X-direction X-direction total total before before after after. In In Figure 11a,b, straightness s along X-direction movement are are shown. five five data data curves curves are are respectively residual residual of four of chosen four chosen models, models, including including proposed B-spline proposed model B-spline three model commonly three used commonly fitting used models fitting models without anywithout. any Three. comparisonthree models comparison are Newton models are interpolation Newton model, interpolation cubicmodel, spline interpolation cubic spline interpolation polynomial model. polynomial From image, model. From it is obvious image, that it is obvious that ability of ability B-spline interpolation of B-spline model interpolation is most significant. model is Newton most significant. interpolation Newton model interpolation has good performance model has good in performance in middle part of stroke but it also shows Runge s phenomenon, which leads middle part of stroke but it also shows Runge s phenomenon, which leads to huge residual to huge residual in beginning end of movement. polynomial model has a in beginning end of movement. polynomial model has a good ability to good ability to follow trend of curve had has a smooth residual curve. With follow trend of curve had has a smooth residual curve. With polynomial model, polynomial model, process avoids steep fluctuation of velocity acceleration. process avoids steep fluctuation of velocity acceleration. cubic spline cubic spline interpolation model has same capability of eliminating with B-spline interpolation model has same capability of eliminating with B-spline interpolation model interpolation model in most positions but still has defects in beginning end of in prediction. most positions B-spline but still interpolation has defects in method beginning is a kind of sectional end interpolation, of prediction. which means B-spline interpolation shape of method interpolation is a kindcurve of sectional only depends interpolation, on nearby which points. means Thus, shape curve of will interpolation not be curve influenced only depends by Runge s on nearby phenomenon. points. Thus, massive curve jump willof not both be velocity influenced byacceleration Runge s phenomenon. does not exist massive in jump interpolation of both velocity curve because acceleration of local doesadjustable not exist inproperty, interpolation B-spline curve interpolation because ofcan local fit irregular adjustable property, data without B-spline influence interpolation fitting canof fitor irregular flat areas. data In conclusion, without influence fitting of orbased flat areas. on In conclusion, B-spline interpolation model based improves on B-spline machine s interpolation straightness model improves accuracy to a machine s remarkable straightness level. Y-direction accuracy to a remarkable Z-direction level. straightness Y-direction s in X-direction Z-direction straightness movement s reduce into 2.92 X-direction μm 5.27 movement μm from reduce 46.55 to μm 2.92 µm80.62 μm. 5.27 µm total fromstraightness 46.55 µm 80.62 accuracy µm. has total been straightness improved by accuracy around 95%. has been improved by around 95%. Figure Figure 11c 11c shows total total of of X-direction position in space movement with a stroke of x = of 2000 x = 2000 mm, mm, y = 1000 y = 1000 mm mm z = 200 z = 200 mm. mm. Two Two curves curves respectively respectively show show before before after. after. max value max value before before appears appears at at end end of of stroke, stroke, with with around around 110 µm. 110 After μm., After, total total shows shows a trend a of trend fluctuating of fluctuating enlargement, enlargement, with with max of 9.88 μm. In conclusion, after aiming at dominant max of 9.88 µm. In conclusion, after aiming at dominant components components calculated by VLGSA method, total of X direction in space movement calculated by VLGSA method, total of X direction in space movement reduces around reduces around 90%. remarkable accuracy improvement of machine shows validity of 90%. remarkable accuracy improvement of machine shows validity of using VLGSA using VLGSA in finding key s. Meanwhile, considering VLGSA method is applied to in finding key s. Meanwhile, considering VLGSA method is applied to MBS model, MBS model, model accuracy is also verified. model accuracy is also verified. 5. Conclusions

Processes 2020, 8, 906 19 of 21 5. Conclusions Structural motional geometric s widely exist in CNC machines. y are most common factors reducing precision machine accuracy, affecting quality of products. refore, characterizing trends, finding most dominant items finally decreasing total geometric play an essential role in development of modern industry. However, considering complexity of precision CNC machine tools, process of targeted can be sophisticated. following works have been done in this paper to achieve goal of improving machine accuracy. 1. A geometric model was established with MBS method. By building topological structure relationship, sophisticated transmitting was decoupled into calculation of several matrices, which can be measured separately. 2. Based on MBS model, items dominant most to total geometric were identified by value leaded global sensitivity analysis (VLGSA), in which actual range of s was considered in process of sensitivity analysis. interaction of different parameters was analyzed in proposed analysis method. 3. Appropriate fitting methods were applied to measured data. data was measured by a high-precision dual-frequency laser interferometer. Considering its excellent ability to describe continual material significant usage in CAD, B-spline interpolation method developed on B-spline curve was used to establish straightness model. 4. A real-time online system was built based on Beckhoff TwinCAT servo system. A experiment was carried out using designed system. straightness in X-direction movement was compensated with four models to verify ability of B-spline interpolation model. X-direction in a space movement of x = 2000 mm, y = 1000 mm z = 200 mm was compensated as a sample. remarkable reduction of geometric in X-direction showed above research was effective. However, this research is based on CNC machine operating at a single speed in cold start temperature. working condition is more sophisticated in actual manufacturing process. Thus, furr works need to concentrate more on relationship of geometric s working condition, such as acceleration changing workforce. Author Contributions: Conceptualization, H.L. Q.C.; methodology, X.Z. Q.C.; software, Q.C. Y.Q.; validation data curation, Q.C. Q.L.; writing-original draft preparation, Q.C.; writing-review editing, Q.C. Y.Z. All authors have read agreed to published version of manuscript. Funding: This research was funded by National Natural Science Foundation of China (Grant no. 51675393), Hubei Province Intellectual Property Guidance Development Project (Grant no. 2019-1-76) Excellent Dissertation Cultivation Funds of Wuhan University of Technology (Grant no. 2018-YS-035). Conflicts of Interest: authors declare no conflict of interest. Appendix A For n+1 control points P i (i = 0, 1,..., n) a knot vector U = {u 0, u 1,..., u m }, a feature polygon can be formed by connecting control points. A k + 1 order B-spline curve can be expressed as follow (2 k n + 1): n C k (u) = N i,k (u)p i (A1) i=0 where N i,k (u) is k order B-spline base function, whose recursive function is shown in Equation (A1). N i,0 (u) = N i,k (u) = (u u i)n i,k 1 (u) u i+k u i (A2) { 1, ui u u i+1 0, else + (u i+k+1 u)n i+1,k 1 (u) u i+k+1 u, u i+1 k < u < u k+1

Processes 2020, 8, 906 20 of 21 References 1. Sarafraz, M.M.; Peyghambarzadeh, S.M. Experimental study on subcooled flow boiling heat transfer to water diethylene glycol mixtures as a coolant inside a vertical annulus. Exp. rm. Fluid Science 2013, 50, 154 162. [CrossRef] 2. Goodarzi, M.; D Orazio, A.; Keshavarzi, A. Develop nano scale method of lattice Boltzmann to predict fluid flow heat transfer of air in inclined lid driven cavity with a large heat source inside, Two case studies: Pure natural convection & mixed convection. Phys. Stat. Mech. Appl. 2018, 509, 210 233. 3. Sarafraz, M.M. Experimental investigation on pool boiling heat transfer to formic acid, propanol 2-butanol pure liquids under atmospheric pressure. J. Appl. Fluid Mech. 2013, 6, 73 79. 4. Goodarzi, M.; Amiri, A.; Goodarzi, M.S. Investigation of heat transfer pressure drop of a counter flow corrugated plate heat exchanger using MWCNT based nanofluids. Int. Commun. Heat Mass Transf. 2015, 66, 172 179. [CrossRef] 5. Chen, J.X.; Lin, S.W.; Zhou, X.L. A comprehensive analysis method for geometric of multi-axis machine tool. Int. J. Mach. Tools Manuf. 2016, 106, 56 66. [CrossRef] 6. Bai, H.; Shen, J.X. Synsis modeling of CNC machine tool based on multi-body system ory. Aeronaut. Manuf. Technol. 2017, 1, 117 121. 7. Paul, R.P.; Zhang, H. Computationally efficient kinematics for manipulators with spherical wrists based on homogeneous transformation representation. Int. J. Robot. Res. 1986, 5, 32 44. [CrossRef] 8. Veitschegger, W.K.; Wu, C.H. Robot accuracy analysis based on kinematics. IEEE J. Robot. Autom. 1986, 2, 171 178. [CrossRef] 9. Zhao, Q.Q.; Jun, H.; Liu, Z.G. Modeling method on motive axes transfer chain for machine tool of arbitrary topological structure. J. Mech. Eng. 2016, 52, 130 137. [CrossRef] 10. Zhou, X.T. Research on Machining Accuracy Prediction Based on Machine Tool Error Modeling. Master s sis, Shong University, Shong, China, 20 March 2016. 11. Chhatre, S.; Francis, R.; Zhou, Y.H.; Titchener, H.N.; King, J.; Keshavarz, M.E. Global sensitivity analysis for determination of parameter importance in bio-manufacturing processes. Biotechnol. Appl. Biochem. 2008, 51, 79 90. [CrossRef] 12. Cossarini, G.; Solidoro, C. Global sensitivity analysis of a trophodynamic model of Gulf of Trieste. Ecol. Model. 2008, 212, 16 27. [CrossRef] 13. Fan, J.W.; Tao, H.H.; Wu, C.J.; Tang, Y.H. Method of sensitivity analysis of machine tools. J. Beijing Univ. Technol. 2019, 45, 314 320. 14. Yan, Y.; Wang, G.; Zhang, F.P. Precision assembly geometric sensitivity analysis based on transformation model for precision assembly. Trans. Beijing Inst. Technol. 2017, 37, 682 686. 15. Li, J.; Xue, F.G.; Liu, X.J. Geometric modeling sensitivity analysis of a five-axis machine tool. Int. J. Adv. Manuf. Technol. 2015, 82, 2037 2051. [CrossRef] 16. Fan, K.G.; Yang, J.G.; Yang, L.Y. Orthogonal polynomials-based rmally induced spindle geometric modeling. Int. J. Adv. Manuf. Technol. 2013, 65, 1791 1800. [CrossRef] 17. Wang, W.; Zhang, Y.; Yang, J.G. Geometric rmal for CNC milling machines based on Newton interpolation method. Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. 2012, 227, 771 778. [CrossRef] 18. Tuo, Z.Y.; Huang, Y.Q.; Shen, M.W. Modeling of geometric rmal complex position of CNC machine tools based on cubic spline interpolation. J. Shanghai Jiao Tong Univ. 2016, 50, 668 672. 19. Li, Y.H.; Huang, T.; Derek, G.C. An approach for smooth trajectory planning of high-speed pick--place parallel robots using quantic B-splines. Mech. Mach. ory 2018, 126, 479 490. [CrossRef] 20. Guo, C.; Yang, L.; Li, Q.Y. Modeling of volumetric s for NC machine tools on multi-body system ories. Mach. Des. Manuf. 2005, 3, 123 125. 21. Liu, Z.F.; Li, D.D.; Liu, Z.Z. Gantry machining tool assembly method predictive optimization based on multi-body system. Comput. Integr. Manuf. Syst. 2014, 20, 394 400. 22. Fan, J.W.; Tang, Y.H.; Wang, H.L. Error modeling identification technology for a CNC camshaft grinder. J. Vib. Shock 2017, 36, 257 264.

Processes 2020, 8, 906 21 of 21 23. Yu, D.J.; Li, R. Application of Sobol method to sensitivity analysis of a non-linear passive vibration isolators. J. Vib. Eng. 2004, 17, 210 213. 24. Gao, X.; Tong, H.; Zhou, L. Geometric of CNC machine tools based on G-code modification. Manuf. Technol. Mach. Tool 2015, 1, 57 62. 2020 by authors. Licensee MDPI, Basel, Switzerl. This article is an open access article distributed under terms conditions of Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).