Variation propagation modelling in multistage machining processes using dual quaternions

Variation propagation models play an important role in part quality prediction, variation source identification, and variation compensation in multistage manufacturing processes. These models often use homogenous transformation matrix, differential motion vector, and/or Jacobian matrix to represent and transform the part, tool and fixture coordinate systems and associated variations. However, the models end up with large matrices as the number features and functional element pairs increase. This work proposes a novel strategy for modelling of variation propagation in multistage machining processes using dual quaternions. The strategy includes representation of the fixture, part, and toolpath by dual quaternions, followed by projection locator points onto the features, which leads to a simplified model of a part-fixture assembly and machining. The proposed approach was validated against stream of variation models and experimental results reported in the literature. This paper aims to provide a new direction of research on variation propagation modelling of multistage manufacturing processes.


Introduction
In the past two decades, variation propagation modelling in multistage manufacturing processes has been actively researched. A significant contribution has been made in the directions of part quality prediction [1][2][3][4], variation source identification [5,6], variation compensation [7,8], and process-oriented tolerancing [9] in multistage manufacturing processes. These contributions are often studied under the concept of stream of variation (SoV) and model of a manufactured part [10]. The concept of SoV may apply techniques such as differential motion vector [11], equivalent fixture error [6], and kinematic analysis [12]. The concept of the model of the manufactured part applies small displacement torsor (SDT) [10,13].
Despite the wide range of applications and extensions of these models, the fundamental mathematical modelling tools remained the same. All of the aforementioned models use homogenous transformation matrix (HTM), differential motion vector, and/or Jacobian matrix to represent and transform the part, tool and fixture coordinate systems and associated variations [10]. The variation propagation model's complexity, matrix dimensions, and number of coordinate system significantly increase as the number of features and machining stations increase, especially when considering generic fixtures, locating surfaces, and machining error [1,4,14,15]. In the case of SoVrelated methods, for instance, the expression for virtual inspection requires square matrices whose dimensions are the number of features multiplied by 6. In the case of SDT-related methods, the Jacobian matrix that relates geometrical relationship between the functional requirement and functional element pairs is of the size 6 by 6 repeated by the number functional element pairs. These matrices can quickly become unwieldy with increase of the number of features and functional element pairs.
As an alternative, relatively less known mathematical tools, quaternions and their extension dual quaternions, are proposed in this paper. Quaternions are mathematical operators that use four parameters to rotate an object [16]. To translate an object, an extra quaternion is added with the help of dual algebra [17]. The two quaternions, when expressed in dual algebra, are collectively referred as a dual quaternion (DQ). Dual quaternions have 8 parameters, which is significantly lower than the 16 parameters of HTM. The HTMs are generated from a sequence of Euler angle-based operations. Using DQs, however, only a set of an axis and angle is needed to rotate an object. Furthermore, the same mathematical tools can also be used to represent the features of a part. Such advantages lead to a better computational efficiency of DQs relative HTM-based operations [18,19]. Furthermore, since the same mathematical tools are used to represent and transform a part, fixture, and toolpath, the mathematical cumbersomeness associated in the operations is reduced. Given such advantages, dual quaternions have been used in wide range of applications such as in robotic controllers [20], face recognition [21], and object animation [22].
In this paper, we derive a mathematical model required to predict a part quality produced in multistage machining process using DQs. The first step is representing the part with DQs, followed by assembly of the part to a fixture, and machining by a toolpath plane. To assemble the part to a fixture, the distance between locators and the projected point onto corresponding datum feature reduced while considering the 3-2-1 fixture layout constraints.
The use of DQs in transformation, representation, and projection operations, associated with parts and fixtures, leads to a completely different strategy to solve the same modelling problem addressed by variation propagation models such as SoV and SDT. This strategy brings forward a novel application of DQs in variation propagation modelling by utilizing techniques using computational geometry and computer graphics, which branches out from the techniques used in SoV and SDT.
The paper is organized as follows. In Section 2, the background to quaternions and variation propagation in manufacturing are provided. In Section 3, the steps in variation propagation are discussed. Sections 4 and 5 present the case study and discussion. Section 6 presents the conclusion of the results.

Quaternions
Quaternions are mathematical tools that use four parameters to rotate an object in a 3D space [16]. The quaternions have two parts: the real part q 0 ∈ ℝ and the vector part q, which can be denoted by where q 0 , q 1 , q 2 , and q 3 are quaternion quotients; i, j, and k are the unit vectors with properties i 2 = j 2 = k 2 = ijk = − 1.

Dual quaternions
Dual quaternions are extensions of quaternions based on dual algebra [17]. Dual quaternions embed the rotational and translational aspects of rigid transformation, simultaneously. Dual quaternions have two parts: the real component e q r and dual components e q d .
where ϵ is a dual operator such that ϵ ≠ 0, ϵ 2 = 0 The DQ representation of pure rotation is the real component. The rotation operator b R is denoted by, The DQ representation of pure translation is the dual component. The translation operator b T is denoted by where t ∈ ℝ 3 is the vector part of the quaternion representing the distance between corresponding points on two objects. The combined transformation operator is obtained from the multiplication of pure rotation and translation.
where t ∈ ℝ is the norm of t The new pos of an object b A after transformation by b e is denoted by where b e * ¼ e R r −ϵ t 2 Equation (6) is equivalent to the 4 by 4 HTM and Eq. (7) is equivalent to a transformed object to new position when operated by HTM. Table 1 summarizes the quaternion and dual quaternion operations, exemplified by two quaternions e p and e q, and two DQs b p and b q. A more detailed description of quaternions and dual quaternions can be found in [23][24][25][26][27][28].

Screw displacement
The above equations assume the axis of rotation passes through the origin of the reference feature. Screw representation combines the translation and rotation in one formulation. This is achieved by transforming the total transformation operator Eq. (6) is multiplied by translation operator Eq. (5), which is the operator associated with the distance of the axis from the reference origin [28]. Thus, for given point p on the axis of rotation, the combined rotation and translation in the direction of the c RT is given by where n is direction of displacement, and t is subsequent translation following the rotation by the angle between vectors n 1 and n 2 It is worth noting that the total transformation of an object using Eq. (8) and Eq. (6) gives the same result. The choice of the operators is, therefore, a matter of convenience with respect to a specific application. However, when the desired transformation has an axis of rotation and direction of translation that are not collinear, Eq. (6) is the recommended option.

Representation of geometric primitives
Many engineering products can be constructed from a combination of planes, lines, and points. These geometric primitives can be represented by DQs, as illustrated in Fig. 1. A plane is represented by its normal and shortest distance from its coordinate system [29].
For a plane normal π r ∈ ℝ 3 (the vector part of the real component of a dual quaternion), and the shortest distance to the plane from the origin d ∈ ℝ, a plane is represented as A line can be represented by two methods, mainly by applying the Plücker coordinate formulation [30]. In the first method, a line b Λ that has points a ∈ ℝ 3 and b ∈ ℝ 3 is represented by its direction U ¼ ða−b ) and moment around the origin (a × b). A line can also be determined from the direction and a point, as shown in Fig. 1 In the second method, a line is obtained from the intersection of two planes [25]. For two given planes b where π 1;r and π 2;r are the vector components of the real part of b π 1 and b π 2 , respectively; π 1, 0, d and π 2, 0, d are the scalar components of the dual parts of b π 1 and b π 2 , respectively. Moreover, the intersection point of a plane and line can be determined from the direction of a line and any point on the line [25]. The intersection point I given the direction of a line U and a point on the line a can be obtained by where π r ¼ n ; is the dual part of the line b Λ; the real part of the line is Λ r ; π 0, d is the scalar component of the dual part of the plane b π.

Variation propagation modelling
The main purpose of variation propagation modelling is to describe the part quality after machining in multiple stages mathematically, thereby applying the results in reducing the variation in key product characteristics (KPCs). The variation modelling approach proposed in this paper is akin to the model of the manufactured part, which applies the techniques of assembling a part to a fixture. In the same line of thinking, a DQ-represented part is assembled to a fixture, followed by a machining by a DQ-represented plane, and inspection at a virtual inspection station.

Part representation
As aforementioned, many engineering parts can be represented by a set of planes, lines, and points. Specifically, prismatic parts can be constructed by a series of planes. These planes are represented by DQs using a plane's normal and shortest distance to the origin. These variables can conveniently be stored in the form of matrices when implementing. Mathematically, for a part b π 1…M with M, number of features can be represented as where n 1…M and ϵd 1…M are normals and distances to the planes. A rectangular prism with six surfaces, for example, can then be represented as In multistage machining processes, the datum and machining features are likely changed after each stage. A specific feature can be referred by its index. In this paper, b S π , b S σ , and b S τ , with indices (π, σ, τ) ∈ (1…M), are used to refer to the primary, secondary, and tertiary datum features, respectively; similarly, b S μ (μ ∈ (1…M)) is used to refer to the machining feature.

Deriving the direction of motion
The direction of motion of a part can be determined from the intersection of planes that pass through the punctual locators. The intersection of the primary and secondary plane determines the direction of motion of the part towards the tertiary locator. Similarly, the intersection of the primary and tertiary plane determines the direction towards the secondary locators. Figure 2 shows the direction of motion of a part, derived from the intersection of datum fixture planes.
Mathematically, the normal vector of a plane passing through the three points of primary locator L 1…3 ∈ ℝ 3 is obtained by For the secondary datum, the normal of the plane that passes through two secondary locators L 4, 5 ∈ ℝ 3 and normal to the primary datum can be computed by For tertiary fixture, the plane that passes through the last locator, L 6 , is perpendicular to the intersection line between the primary and secondary fixture planes. The tertiary plane is formed by the plane normal of both primary and secondary planes F 3; r . Thus, the dual quaternion representation of primary, secondary, and tertiary locators' planes becomes

Assembling part to fixture
In a 3-2-1 fixture layout, a part is assembled to a fixture in a sequence of primary, secondary, and tertiary datum features.
Hereafter, the rotation, translation, and total transformation operators associated with primary, secondary, and tertiary fix- , and b e h 1…3 , respectively.

Assembling to primary locators
In assembling the part to the primary locators, rotation and translation steps are required. The steps are illustrated graphically in Fig. 3. Mathematically, in assembling the primary datum, the rotation required between the primary fixture b  can be obtained by The translational operator b T h 1 required to close the error between the parallel planes of the primary fixture and primary feature, using Eq. (5), can be denoted as Alternatively, b T h 1 can be obtained as the difference between the primary locators and the projected point on the primary datum feature using Eq. (12). Figure 4a illustrates the projection of primary locators onto the primary datum feature, from which the distance from the locators' points can be computed. Thus, The combined rotational and translational transformation is obtained from the multiplication of rotational operator b R 1 and translational operator b T 1 . Thus, the total transformation operator b e 1 required in assembling the primary datum feature to the primary locators becomes The operator b e h 1 is then used to transform all features of the

Assembling to secondary locators
In assembling the secondary datum feature, the rotation is constrained by the normal vector of the primary locators' plane, which may not necessarily be aligned in the z-direction  The operator b R h 2 helps make the two vectors parallel but not coincident. In this intermediate step, after the rotation, the part can be represented by Since translation of the part towards the secondary locators is not necessarily in orthogonal direction, the direction of the line of intersection of the primary and tertiary locators' planes F h 2;r (Eq. (17)) is used (Fig. 6). Thus, in transformation of the part towards the secondary direction using a single locator, say

Assembling to tertiary locators
In sliding the part towards the tertiary plane, the direction of the line where the planes of the primary and secondary loca-

Including machining variation
The variation propagation models that overlook the machining errors can yield inaccurate prediction. In this regard, Abellan-Nebot et al. [2] extended the classical stream of variation models to include machining errors. The machining error is the cumulative sum of the rotational and translational effects of four factors: (1) geometric, kinematic, thermal variations of the machine tool axes, (2) spindle-thermal variations, (3) the cutting force-induced variation, and (4) cutting tool wear-induced variation [2]. These effects are added to traditional stream of variation as composition of coordinate systems through the application of HTM. The approach can be reformulated using DQs for representation and transformation operations.
To reformulate the approach, it is assumed that the effects of these individual errors are accumulated in a toolpath that can be represented by a plane. Reformulating Eq. (4) presented in [2], the total effect of the individual factors can be computed as the product of the individual screw operators. Thus, where the subscripts a, st, f, and w indicate machine-tool axes, spindle-thermal, cutting-force, and cutting tool wear-induced variations, respectively. The variables θ, δ, n, and p are the rotational deviation, translational deviation, an axis of rotation/translation, and a point on the axis, respectively. Thus, for a plane made by a nominal toolpath b μ 0;h ¼ n 0;h μ þ ϵd 0;h μ , the pos of the resulting toolpath becomes Since the machining feature assumes the toolpath plane's pos after machining, while the part is assembled to the fixture, the initial machining feature can be replaced by the toolpath plane. Thus, the toolpath plane b μ h replaces the machining feature to obtain the new machined feature b S h;′′′ μ , as shown in Fig. 8.
Equation (35) updates the DQ representation of the part b S h;′′′ 1…M while still assembled to the fixture at station h. In the subsequent stations, the part is rotated by an angle, say 180°a round the X-axis, which can be donated as b e hþ1 ¼ cos When the part is transformed for use in subsequent stations, the primary, secondary, and tertiary datum features change. Thus, the indices of the equations above should correspond the new datum features of the stations.

Virtual measurement at inspection station
The machined part is moved to inspection station by transforming the reference datum of the part to a flat granite. A proper constraint of the part can be applied following the  assembly steps discussed in Section 3.3. Focusing on the pri- Once the part is setup on a test station, measurement is conducted using a coordinate-measuring machine (CMM). Given the CMM measurement direction U h 1…M and origin of inspection points Υ h 1…M ∈ℝ 3 at a nominal part b S 0;h 1…M , the actual points can be obtained by projecting inspection points onto their respective planes, in the direction of probing, as shown in Fig. 9. Thus,

Vertices and visualization
As presented above, a part can be represented by a set of planes. Since the planes are infinite in dimension, the dimensions can be truncated at the intersection of the planes using Eqs. (10) and (11). The lines intersect to a third plane at a vertex, which can be obtained using Eqs. (12) and (13). Thus, the intersection line of two planes, say b S 1 (for π =1) and b S 5 , can be expressed as The intersection of the line b Λ In terms of implementation, the intersection points are obtained using a loop that checks all combination of plane pairs that make lines and the intersection of the lines to the rest of the planes.

Case study
A case study was selected to compare and validate the DQbased modelling approach, by treating the part model (Fig. 11) presented in [1][2][3] as a base model. The case study demonstrates variation propagation in multistage machining process, with and without considering machining errors. The case study used the inputs from the work reported in [2] to predict the dimensional error of the part after machining. The main goal of the case study was to show the prediction improvement when including machining errors in multistage   Table 2. Furthermore, in Table 2, the equivalent DQ representation of the features is presented.
The features S 2 , S 1 , and S 3,8 were machined in stations 1, 2, and 3, respectively, as shown in Fig. 12. Three experimental conditions were set such that the locator deviations and machining variation-due to thermal and tool wear-induced deviations-are considered, as shown Table 3. Specifically, the change in spindle temperature of 15°C corresponded to a deviation of − 0.0744 mm on the part [2]. Furthermore, a 0.3mm cutting tool wear on the primary and the secondary cutting edges corresponded to deviations of 0.03744 mm and 0.0404 mm on the part, respectively [2]. Feature S 8 was machined by the secondary cutting edge of the tool.
To assemble the part to fixture, the locator points were projected on the corresponding part features at each station using Eq. (12). For the first station, the primary datum feature k + ϵ(0 1 × 4 ) was used to assemble to the plane made by locators L 1…3 using Eqs. (22)- (25). In assembling the part to the secondary locators, the projection of secondary locator points to corresponding datum, obtained from Eqs. (26)- (29), was used to determine the required total transformation of the part to the second new position. Again, using this new position, the part was transformed by the distance between the tertiary locator and its projection using Eqs. (30) and (31).
Moreover, to compute the machining variation error, the operator b e 4 from Eq. (33) was used. The screw operator for machine tool axes and cutting force-induced variations were set to [1 0 0 0 0 0 0 0]. The translation operators of 1 − ϵ0.0744k for spindle-thermal variations, and 1 + ϵ0.03744k and 1 + ϵ0.0404k for cutting tool wear-induced variation were used. These operators were then multiplied to their corresponding nominal machining plane to yield the results shown in Table 3.
Once the assembly to the fixture and replacement of the machined surface were performed, the part was transformed to an inspection station, to which feature S 1 was assembled to granite, with a DQ representation of k, using Eqs. (37) and (38). Using a virtual CMM points in Table 2 as test points, each of the points was projected onto their respective machined features using Eq. (39). The projected points became the predicted measurement values of S 2 , S 3 , S 8 , which were 44.999, 40.032, 85.008; 45. 003, 40.223, 84.986; and 44.966, 40.193, 85.027 for experiments 1 to 3, respectively. Using these values, the dimensions of KPC 1…3 were extracted by subtracting the points of S 2 and S 3 from that of S 0 , and S 8 from that of S 6 .
These results, when compared with the classical SoV and extended SoV, were equal in experiments 1 and 2, and 20% better in experiment 3 in predicting the actual experimental results reported in [2]. Figure 13 summarizes the prediction error achieved by the SoVs-and DQ-based approaches with respect to the actual experimental results.

Discussion
Dual quaternions are excellent way to combine both rotation and translation motions as well representation of the geometrical objects. In the proposed approach, part features, fixtures, and lines have been represented and transformed by DQs. In the demonstration case study, the two experiments based on SoV, and its extension that improved prediction significantly by including machining error [3], were compared with that of DQ-based approach. The part quality prediction based on DQ was equal to that of the two SoVs' prediction with respect to experimental results when machining errors were not considered. When machining error was considered, the prediction based on DQ was 20% lower prediction error compared to extended SoV with respect to experimental results. However, the rounding errors could be behind the better result achieved by the DQ-  The advantage of DQ-based approach is that locators' heights are assumed independent, while the same assumption cannot apply to stream of variation approach. Furthermore, the implementation of the case study using SoV approach [3] gave the largest matrix size of 54 by 54, while the DQ-based approach kept the initial 9 DQs (9 by 8 matrix) that represented the part. This is a significant reduction in the resulting matrix size.
Beside the mathematical elegance, DQs are computationally efficient. For comparison, the MATLAB implementations of the 3 experiments using SoV- [3] and DQ-based approaches, both of which were implemented without consideration for execution speed when both implementations were executed in an Intel® Core™ i7-8665U 1.90-GHz machine. The DQ-based implementation yielded a result in average speed 0.05 s, whereas the SoV implementation in 0.09 s, when machining error is not considered. When machining error is considered, the DQ-based approach yielded result in 0.07 s and the extended SoV in 0.35 s. This means the average execution time of the DQ-based approach is 32% of that of (extended) SoV-based predictions. The computational advantage can be marked when repeated experiments are conducted, such as in statistical tolerance analysis.
Moreover, the proposed approach applies projection of points onto surfaces, which is in line with the ray tracing methods. Such approaches are independent of shape complexity of an object. Thus, the proposed approach has potential to be used for generic shapes. For generic cases, locator points can be projected onto a tangent plane of a more complex shape, such as cylindrical shapes.
Despite the stated benefits, DQs have weaknesses too. The major inconvenience of working with dual quaternions is that the b q = -b q. The computed axis of rotation between two planes yields an arbitrary direction to leading arbitrary angles. In such cases, the angle of rotation should be rotated in opposite direction by flipping the direction of axis of rotation. Specifically, when comparing two planes, and analysis there on, the direction of the dominate vector should be kept in same sign. Furthermore, care must be taken when using multiple libraries, as they may have different coordinate systems.

Conclusion
This work derived a model of variation propagation in multistage machining processes using DQs. The use of the DQs led to a new strategy in modelling the variation propagation, which branched off from the techniques used in classical stream of variation. The technique involves reducing the distance between locator points and their projections of onto their corresponding DQ-represented datum features, while satisfying the constraints of fixture layouts. Once the part is assembled to a fixture, a DQ-represented toolpath plane replaces the machining feature, followed by transformation to virtual inspection station for measurement at specific points on the features. The DQ-based approach gave a 20% lower average prediction error compared to extended SoV with respect to experimental results. Moreover, for the same inputs, the execution time was reduced by 68% and the largest matrix size was reduced from 54 by 54 to 9 by 8. Since the method is based on projection points on surfaces, it has a potential to apply to more complex shapes, while maintaining the mathematical simplicity. In this direction, it is worth investigating the application of such techniques in quadric shape and NURBS-represented manufactured parts, especially the parts with freeform surfaces.
Funding Open access funding provided by Royal Institute of Technology. This work is financially supported by VINNOVA with grant number 2019-03570. The authors thank VINNOVA for supporting the project.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.