Microsoft Word - 271-1905-1-LE-rev ACTA IMEKO  ISSN: 2221‐870X  December 2015, Volume 4, Number 4, 20‐25    ACTA IMEKO | www.imeko.org  December 2015 | Volume 4 | Number 4 | 20  Convergence of a finite difference approach for detailed  deviation zone estimation in coordinate metrology  Ahmad Barari, Saeed Jamiolahmadi  Faculty of engineering and applied science, University of Ontario Institute of Technology, Oshawa, Ontario, L1H 4K7, Canada      Section: RESEARCH PAPER   Keywords: Distribution of Geometric Deviations; coordinate metrology; finite difference method  Citation: Ahmad Barari, Saeed Jamiolahmadi, Convergence of a finite difference approach for detailed deviation zone estimation in coordinate metrology,  Acta IMEKO, vol. 4, no. 4, article 6, December 2015, identifier: IMEKO‐ACTA‐04 (2015)‐04‐06  Section Editor: Franco Pavese, Torino, Italy  Received April 13, 2015; In final form October 31, 2015; Published December 2015  Copyright: © 2015 IMEKO. This is an open‐access article distributed under the terms of the Creative Commons Attribution 3.0 License, which permits  unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited  Funding: This work was supported by Natural Sciences and Engineering Research Council (NSERC) of Canada  Corresponding author: Ahmad Barari, e‐mail: Ahmad.Barari@uoit.ca    1. INTRODUCTION  Although surface coordinate metrology is a key inspection tool in many different engineering and scientific applications, unfortunately it is still influenced by several concerning sources of uncertainties which reduces the reliability of this commonly assumed precise metrology process. The traditional coordinate metrology is completed by sequential performing of three computational tasks of Point Measurement Planning (PMP), Substitute Geometry Estimation (SGE), and Deviation Zone Evaluation (DZE). There are several sources of uncertainties experienced and studied by researchers in conducting each one of the above tasks. In such systems, the inspection accuracy is subject to sorts of “plug-in uncertainty”, i.e. the uncertainty due to estimating some aspects of a probability distribution on the basis of sampling. Sprauel et al showed that uncertainties in CMM measurements follow a normal distribution. They governed the parameters of the surfaces under inspection to characterize the intrinsic error parameters of CMMs [1]. Such sources of uncertainties can adversely affect the accuracy and reliability of the inspection. Recently, an integrated inspection system has been studied as an efficient solution to combat the inherent plug-in uncertainty in coordinate metrology. Choi and Kurfess developed a numerical algorithm based on tolerance zone to interpret the data extracted from CMMs. They showed that a zone fitting algorithm provides more conformance to the tolerance zone if tolerances are characterized by such span [2], [3]. The number of measured points in any coordinate metrology process is limited in increasing the number of measurements, increases the inspection time and it is costly. The typical hard-probing methods usually are very slow in capturing the points. The number of measurement points using these devices in practice is usually even less than one hundred. It is much less expensive to measure more points using optical ABSTRACT  In order to comprehend an entire surface’s deviation zone, infinite measured points are required. Using the common measurement  techniques through coordinate metrology, a limited number of surface actual points can be acquired. However, the obtained points  would  not  provide  sufficient  information  to  examine  the  geometry  thoroughly.  Reliability  and  convergence  of a  novel  solution  to  predict surface behaviour via Distribution of Geometric Deviations (DGD) is studied in this paper. The methodology governs the mean  value property of the harmonic functions to solve the Laplace equation around each measured point. This DGD model can be used to  reconstruct surface deviation values at any unmeasured point of the inspected surface based on a limited number of measured points.  The convergence of the introduced approach is studied in this paper. A complete approach to implement the developed methodology  is described, and the validation process is studied using actual case studies and mathematical functions. This methodology is practical  in closed‐loop inspection and manufacturing processes to form a scheme for compensating the surface errors during manufacturing  process based on the DGD model. ACTA IMEKO | www.imeko.org  December 2015 | Volume 4 | Number 4 | 21  sensors, but unfortunately measurement accuracy in these sensors is still dramatically low. Having this limitation, it can be beneficial to design coordinate metrology processes to utilize estimation techniques that helps to model details of the deviation zone for the locations that are not physically measured. The focus of this paper is on evaluation of a detailed deviation zone when a detailed level of information for the distribution of geometric deviations is required. This level of information is important particularly for an upstream process that needs to be performed on the measured surface. A good example of an upstream is a compensating manufacturing phase that removes or compensates the existing geometric deviations of the measured surface to finish a new surface with the desired and acceptable range of characteristics. A successful inspection process evaluates the minimum geometric deviations of the points and their corresponding points on the best substitute geometry. However, this information is only available for the measured points. Estimation of geometric deviations at the regions that are not measured by coordinate measurements is necessary for a compensation process in the closed-loop of inspection-machining [3]-[7]. The presented method estimates the distribution of geometric deviations based on the pattern of the available data. This method is developed based on the assumption that the detailed deviation zone has a mean-value property similar to a harmonic function. Assigning the surface points to a grid and using the Finite Difference Method derived from the Laplace equation, deviation of the rest of the surface points is approximated. The results show an improvement of approximation by adding more points till a convergence threshold is satisfied. The structure of the paper is as follows: after a brief review of the background information, the developed methodology is presented, and then the implementation of the methodology is explained, followed by a validation case study of an industrial part. 2. BACKGROUND  Three major concepts utilized in developing the presented methodology are: definition of the Laplace equation, Mean- value property, and maximum and minimum principals. A brief review of these concepts is presented. 2.1. Laplace Equation  A second order partial differential equation, which is often in two dimensions, is written as: 2 2 2 2 0, f f x y       (1) is called Laplace equation. If the Laplace equation is solved over a boundary condition, function F which satisfies both the Laplace equation and its boundary conditions, is called a harmonic function. Harmonic functions possess unique properties which are of interest to the developed methodology [8]. 2.2. Mean‐Value Property  If D is a domain of finite measure in the Euclidian space Rn (where n ≥ 2); and, if there exists a point P0 in D for which every H harmonic function in D and integrable over D, the volume mean of H over D is equal to H(P0). Then D is an open ball (or disc for n = 2) with the centre of P0. In other words, if H: U → R is a harmonic function on an open set, it has the mean value property if it satisfies the following relationship [9]: ( ) ( ) ( ) , ( ) . D z h z u w dw D z U   (2) 2.3. Maximum‐Minimum Principal  If H:U → R is a harmonic function on a connected open set U and if there is a point P0 U with the property that 0( ) sup ( )Q Uu P H Q , then u is constant on U. On a corollary of this property the Minimum Principle is derived in the same way for 0( ) inf ( )Q Uu P H Q [9]. In detail, “sup” and “inf” refer to the superior and inferior of the harmonic function. 3. METHODOLOGY  A perfect flat surface can be expressed in two dimensions defined by mutually perpendicular u and v axes. Any out of flatness deviation of this system in surface can be represented in the e direction which is perpendicular to the u-v 2-D plane. Therefore, the flatness coordinate metrology data can be represented in this u-v-e coordinate system, Deviation Space, expressing the detailed deviations of measured points. As illustrated in Figure 1, the deviation of measured point Pi in the u-v-e coordinate system is ei and the point Pi* with coordinate of ui and vi the corresponding point for Pi in the u-v plane. All points Pi belong to the actual measured surface while points Pi* belong to the best substitute geometry that can represent the actual surface. Typically the up-stream design of manufacturing processes uses the best substitute geometry for the required analysis, and considering the overall deviation zones by a deviation boundary over information of the detailed deviation zone can be highly beneficial for processes precision manufacturing and closed-loops of inspection-manufacture. 4. DEVIATION SPACE AND GRID CONSTRUCTION  The substitute geometry is estimated using fitting criteria. The common fitting criteria utilized by industries are the least square fitting and min-max fitting [3]. Other fitting criteria can also be utilized upon the mission of inspection and the up- stream processes. A good example of other fitting criteria is maximum conformance to design tolerances for closed-loop of inspection and machining [6], [7]. Figure 1. u‐v‐e Deviation space, Detailed Deviation Zone representation.  ACTA IMEKO | www.imeko.org  December 2015 | Volume 4 | Number 4 | 22  Upon estimation of substitute geometry for a group of measured points, it is possible to transfer inspection information to the u-v-e coordinate system described above. After the point measurement process and construction of the deviation space, a grid is devised based on the desired resolution of estimation. This resolution can be decided based on the type and nature of utilized coordinate metrology. In the next step, the deviation space corresponding points, Pi*’s, are bijectively assigned to the corresponding grid nodes for further processes. These locations in the deviation space are referred as sites in the rest of this work. Values from the e-axis are considered as potential values for local deviation of the surface. These values are the Euclidian distance from the measured points to the substitute geometry. 4.1. Estimation of Deviations  Estimation of the geometric deviation for an unmeasured point of the actual surface is required when the geometric deviation values for the measured points are available as the attributes of the sites in the deviation space. It is shown in [3], [5] that the distribution of geometric deviations caused by quasi-static manufacturing errors is continuous because it is a direct function of existence continuity in the nominal geometry. Therefore, the geometric deviations can be considered as a continuous variable with known values at sites in the deviation space. With using a proper point measurement process a sufficient number of points from the appropriate locations are measured. Therefore, it is assumed that the sites are available to represent the continuity of geometric deviation function. Geometric deviation of an arbitrary unmeasured location can be estimated by studying its proximity to the known geometric deviations of the measured sites. In order to estimate the detailed deviation zone the conformal mapping method using a harmonic function with Drichlet boundary problem is used. In the current approach, instead of assuming a single shape centre with the potential value of zero, we consider the sites as a group of shape centres and their deviation attributes as the potential function to run the Finite Different Iteration. This way the deviation attributes for the empty nodes between the sites are computed. In order to implement this concept the Laplace equation over the grid is solved. The numerical solution for the Laplace equation can be obtained using Taylor’s series neglecting higher order terms [10]: 2 2 1 12 2 2 1 1 1 1 [ ( , ) ( , ) ( , ) ( , ) 4 ( , )], i i i i i i i i i i f f x y x y x y h x y x y x y                     (3) 2 1 1 2 2 ( , ) 2 ( , ) ( , ) ,i i i i i i x y x y x yf y k        (4) where h and k are the grid step sizes along u and v coordinates respectively, and Ф is the detailed deviation zone function. Considering equal step sizes in u and v directions, the Laplace equation can be rewritten as: 2 2 1 12 2 2 1 1 1 1 [ ( , ) ( , ) ( , ) ( , ) 4 ( , )]. i i i i i i i i i i f f x y x y x y h x y x y x y                     (5) From the subsequent equation a Finite Different Method formula is extracted, if f is considered a harmonic function: 1 1 1 1 1 ( , ) [ ( , ) ( , ) 4 ( , ) ( , )]. i i i i i i i i i i x y x y x y x y x y               (6) To initiate the iteration another equation is defined where j is the iteration index: 1 1 1 1 1 1 ( , ) [ ( , ) ( , ) 4 ( , ) ( , )]. j i i j i i j i i j i i j i i x y x y x y x y x y                (7) In the developed methodology two termination criteria are utilized for the iterations. Considering the typical coordinate metrology accuracy between 10-3 mm and 10-6 mm, a threshold is selected. The ultimate goal of the developed method is to find a smooth distribution of attributes all over the grid. Therefore, the derivative of attribute for every single grid node is calculated and stored in the attribute derivate matrix. The difference in the attribute derivative matrix is calculated in each iteration with its previous stage. The standard deviation of differences in derivative of attributes for all grid nodes is calculated. If the calculated standard deviation reaches to a level below the decided threshold then the iteration process is terminated. The iteration is costly and takes time; therefore, if the iteration has been running more than a specified amount before converging, the iteration process will be terminated. 5. IMPLEMENTATION  The finite difference method can be utilized in the developed methodology only if function Ф is a harmonic function; therefore, it is crucial to prove that there exists a function in the domain of the Laplace equation answer. Since the mean value property is a unique feature of the harmonic functions, if the finite difference method is valid in a domain, then it can be concluded that the governing function in the domain is harmonic and the use of finite difference method is legitimate. This legitimacy is validated for this work by applying the methodology in an actual case study. The Implementation process is as follows: A. Point Measurement Process: an industrial flat surface is measured with either an optical sensor or through hard probing to produce a point cloud; B. Best Plane Fit: the best substitute geometry is constructed using the desired fitting criteria, and their Euclidian distance from the best fitted plane is calculated and assigned to points; C. Mapping to Deviation Space: measured points are mapped to the deviation space. Their corresponding points and their Euclidian distances to the best substitute geometry are considered as their coordinates in the constructed deviation space; D. Gridding: a grid is created based on the desired resolution of the estimation process. Sites are specified in the grid. As a constraint for grid resolution, it is considered that each grid node should not cover more than one site; E. Initialization: the deviation attributes of nodes that contain a site are assigned to them and deviation attributes of all other nodes are equal to zero; F. Percentile Calculation: a series of percentages of points are selected randomly and assigned to the grid node with their potential value. The finite difference method is then applied to the grid; ACTA IMEKO | www.imeko.org  December 2015 | Volume 4 | Number 4 | 23  G. Estimation: the finite difference method is applied to find the deviation attribute of all the empty grid nodes through the iterations. It is worthy to mention that the deviation attribute of the sites remain fixed during the iterations; H. Termination: to compare the results between various percentages used in the model, the method was iterated for 1000 runs in each case. 6. VALIDATOIN AND CASE STUDIES  To validate the process two case studies were studied. The cases consist of one actual manufactured surface and one mathematical surface. The manufactured piece was examined through hard probing to acquire sample points from the surface. In the mathematical case, a sinusoidal function was used to define a 3-D surface; therefore, the deviation error is a function of the position of the points. Using the steps defined in the previous section, the methodology is applied on the surfaces to validate the accuracy and efficiency of the DGD model. 6.1. Case I: Industrial flat surface  As the first case study, an industrial flat surface with over-all dimensions of 30 mm × 30 mm was inspected. The measurement process was conducted using an optical probe on a Faro Coordinate metrology arm. The point measurement planning was based on stratified grid inspection with a step size of 0.5 mm. Through this process 7287 points were captured to generate the point cloud shown in Figure 2. As the fitting process a total least square fitting criteria is utilized. The substitute geometry and corresponding deviations of the measured points are mapped to the deviation space which is illustrated in Figure 3. A grid with size of 320×240 was created and the 7287 sites were attached to it. The iteration process with the maximum number of iterations equal to 1000 was run. Figure 4 to Figure 8 illustrate the developed deviation detailed zone after 50 runs, 100 runs, 500 runs and 736 runs, respectively. As it can be seen in Figure 4 to Figure 8, by increasing the iteration runs the grid is evolving to estimate the deviation attribute of the empty nodes, and it converges to the desired threshold. This convergence is a validation to applicability of the developed methodology for the current problem. Since the finite difference method uses the circular neighbourhood to estimate the centre of the circle, it provides a smooth estimation for the empty nodes through progress of the iteration process. Based on the number of the runs of the finite difference method the grid converges to find the rest of the potential values. As it can be seen in Figure 5 to Figure 8, by increasing the iteration runs which in fact means using more and more the finite difference method, the grid is evolving to find the rest of the nodes based on the centre points. It worthy of mention to note that the potential of the centre points remains fixed during the iteration process. Since the Finite difference method uses the circular neighbour points to estimate the centre of the circle, it provides a smooth estimation of the nodes based on the former points. Figure 9 also demonstrates the convergence process of the estimation. Based on the termination criterion defined, the estimation process is terminated after 735 iterations.   Figure 2. Measured points based on point measurement planning.  Figure 3. Mapping of the measured points to the deviation space.  Figure 4. Deviation zone distribution before iteration runs.  Figure 5. Deviation zone distribution after 50 runs.  ACTA IMEKO | www.imeko.org  December 2015 | Volume 4 | Number 4 | 24  6.2. Case II: Sinusoidal Surface  To test the capability of the proposed methodology in theoretical situations, this case study was considered. The following two directional sinusoidal function was used to generate a hypothetical 3-D surface based upon the assumption that the error deviation is a function of the positions of points as e = f(u,v) [3]. 0.1 sin(2 ) 0.2 sin(2 )e u v       . (8) The surface is constructed on u and v values between -0.5 and 0.5 and then transferred to a grid with the size of 250×250 to build the error model as shown in Figure 10. The error values i.e. sites, vary between -0.2 and 0.2. To initialize the finite difference method 10 % of the original is kept in the grid and the rest were erased. Assuming a metric unit for the coordinate model as millimeter, the iteration process is initiated based on the stages described in Implementation (Section 5) to find the 90 % of the sites that were erased. Similar to the industrial part, the grid starts to converge to the actual surface. This process is illustrated from Figure 11 to Figure 15. The convergence process is terminated at iteration number 202, which is demonstrated in Figure 16. The convergence is ended at an earlier stage regarding to the industrial surface. This basically shows a good efficiency of the algorithm over mathematical functions. 7. CONCLUSIONS  Estimation of the detailed deviation zone is a key task in integrated computational platforms for coordinate metrology systems. The methodology evaluated in this paper is an Figure 6. Deviation zone distribution after 100 runs.  Figure 7. Deviation zone distribution after 500 Runs.  Figure 8. Deviation zone distribution after 735 runs. Figure 9. Convergence of the estimation process for the industrial part.  Figure 10. Two directions sinusoidal surface constructed based on (8).  Figure 11. Deviation zone distribution before iteration runs.  ACTA IMEKO | www.imeko.org  December 2015 | Volume 4 | Number 4 | 25  approach to assess the distribution of geometric deviations on the manufactured surfaces. The convergence of the developed methodology is validated by studying various cases in this paper. The implementation of the detailed deviation zone estimation can be used for upstream manufacturing processes or manufacturing error compensation in closed-loop of inspection-manufacturing. The estimated deviation zone can lead to generate valuable information about the measured surface. This information can be utilized in various design and manufacturing applications. ACKNOWLEDGEMENT  The authors would like to thank the Natural Sciences and Engineering Research Council (NSERC) of Canada for partial funding of this research work. REFERENCES  [1] Sprauel, J. M., et al. "Uncertainties in CMM measurements, control of ISO specifications." CIRP Annals-Manufacturing Technology 52.1 (2003): 423-426. [2] Choi, Woncheol, and Thomas R. Kurfess. "Dimensional measurement data analysis, part 1: a zone fitting algorithm." Journal of manufacturing science and engineering 121.2 (1999): 238-245. [3] A. Barari, H. A. ElMaraghy, and G. K. Knopf, “Search-guided sampling to reduce uncertainty of minimum deviation zone estimation,” J. Comput. Inf. Sci. Eng., vol. 7, no. 4, pp. 360–371, 2007. [4] A. Barari, H. A. ElMaraghy, and W. H. ElMaraghy, “Design for manufacturing of sculptured surfaces: a computational platform,” J. Comput. Inf. Sci. Eng., vol. 9, no. 2, p. 21006, 2009. [5] A. Barari, H. A. ElMaraghy, and P. Orban, “NURBS representation of actual machined surfaces,” Int. J. Comput. Integr. Manuf., vol. 22, no. 5, pp. 395–410, 2009. [6] A. Barari and R. Pop-Iliev, “Reducing rigidity by implementing closed-loop engineering in adaptable design and manufacturing systems,” J. Manuf. Syst., vol. 28, no. 2, pp. 47–54, 2009. [7] H. A. ElMaraghy, A. Barari, and G. K. Knopf, “Integrated inspection and machining for maximum conformance to design tolerances,” in 54th CIRP General Assembly, 2004. [8] B. Epstein and M. M. Schiffer, “On the mean-value property of harmonic functions,” J. d’Analyse Mathématique, vol. 14, no. 1, pp. 109–111, 1965. [9] S. Axler, P. Bourdon, and R. Wade, Harmonic function theory, vol. 137. Springer Science & Business Media, 2001. [10] H. K. Voruganti, B. Dasgupta, and G. Hommel, “A novel potential field based domain mapping method,” in Proceedings of 10th WSEAS Conference on Computers (ICCOMPP06), Athens, Greece, July, 2006, pp. 13–15. Figure 16. Convergence of the estimation process for the sinusoidal surface. Figure 12. Deviation zone distribution after 50 runs.  Figure 13. Deviation zone distribution after 100 runs.  Figure 14. Deviation zone distribution after 150 runs. Figure 15. Deviation zone distribution after 202 run.