Keywords: Geometrical Tolerancing, Skin Model, Tolerance Analysis, Assembly Simulation, Coordinate Metrology, GeoSpelling

Acknowledgement

I have spent three and a half years in I2M, Université de Bordeaux, to do some research on assembly simulation and deviation evaluation. The work was funded by China Scholarship Council. I am extremely grateful to my advisor Alex Ballu. Three years ago, he provided me the opportunity to do my PhD study in France. During these years, he taught me systematically how to carry out research, and the way of critical thinking from a global perspective. Besides, he gives me sufficient flexibility on the topics that I am interested in and provides inspirational suggestions when needed. I would like to express my sincere gratitude for his care about my life, my career, and my whole stay in Bordeaux.

Moreover, I would like to thank Prof. Michel Mesnard, for his in depth sharing of the French culture and language learning which enables me to enjoy fully in a foreign environment. I also want to thank Denis Teissandier, Vincent Delos and Yann Ledoux for their endless help during these years. And to thank the members of thesis defense jury, Nicolas Perry, Nabil Anwer, Antoine Tahan, and Jean François Fontaine, for their great support on the improvement of my thesis.

And my thanks goes to Santiago and Ricardo, not just my colleagues but friends. We entered into the laboratory at the same time. Their companion of the tough days, and the time we spent together for travelling, surfing, skating, and barbecue have enriched my memories. The friendships grows stronger with the days and it would be cherished for a lifetime. Thanks to all the friends in the laboratory, Doriane, Ulises, David, Rachel, Soukaina, Jorge, Sasha, who have brought enormous joy to me during these days.

I also want to thank Yingying, who contacted me three years ago. Thanks to Siyu, Yan and Qiqi for their help when I first come to France. Thanks to Cheng, Yana, Zhengyang, Qiu, Si, Jie, Zuowei, Ting, who have helped me a lot during these years. I want to thank my parents and my wife for their unconditional love and support.

Figures

Résumé Etendu

La géométrie d'une pièce fabriquée réelle diffère de la pièce virtuelle modélisée à l'aide d'un système de Conception Assistée par Ordinateur (CAO). Cette différence est due à l'accumulation d'écarts lors du processus de fabrication, qui sont inévitables. Une mauvaise gestion de cette variation géométrique par rapport à l'intention du concepteur est généralement la source de problèmes d'assemblage et de dysfonctionnement du produit ainsi que de défauts d'ordre esthétique, ce qui fera perdre au produit sa compétitivité sur le marché. D'autre part, une exigence trop serrée sur les écarts de fabrication entraîne une augmentation des coûts, ce qui rend le produit non rentable. Afin de gérer cette problématique, une activité, nommée tolérancement géométrique, est menée afin de définir et gérer les écarts géométriques et d'atteindre un compromis entre qualité et coût.

La méthode classique, la plus simple pour effectuer une simulation des tolérances, consiste à utiliser des chaînes de cotes unidimensionnelles. Cependant, avec des exigences accrues pour une gestion plus fine des écarts géométriques du produit, les méthodes actuelles, basées sur des surfaces de forme idéale, ne sont plus satisfaisantes en termes de précision. Les défauts de forme doivent être pris en compte car de plus en plus de pièces mécaniques nécessitent une grande précision géométrique. De plus, la gestion des tolérances doit être étendue au cycle de vie complet du produit et doit être envisagée le plus tôt possible.

Dans cette mouvance, les normes évoluent vers un plus haut degré de description de la spécification pour améliorer la communication entre les différents intervenants et organisations.

Les modèles et les termes utilisés pour décrire les exigences en termes de précision des produits sont définis dans la série des normes relatives aux spécifications géométriques des produits (ISO GPS

Introduction

The geometry of a real manufactured part differs from the virtual part created using Computer-Aided Design (CAD) systems. This difference is due to the accumulation of deviations during manufacturing processes, which are unavoidable. Careless management of this geometry divergence from design intent will ultimately cause problems in product assemblability, product functioning and even esthetic requirements, causing the product to lose its competitiveness in the market. On the other hand, too tight control of manufacturing deviations causes a surge in costs, which make the product unprofitable. In this scenario, work called tolerancing is conducted in order to define and manage geometry deviations and reach a compromise between quality and cost.

The simplest method for conducting tolerance stack up is to use the 1D tolerance chain, which has already been used in industry for more than 50 years [START_REF] Charpentier | A scientific point of view of a simple industrial tolerancing process[END_REF]. This method calculates extreme limits for defects, which could generate too strict requirements on product quality. Statistical methods like root-sum-square [START_REF] Bender | Statistical tolerancing as it relates to quality control and the designer[END_REF] and system moment [START_REF] Evans | Statistical tolerancing: The state of the art. I-Background[END_REF] are introduced to enhance the tolerancing method. Statistical-based methods allow a certain number of parts to have defects outside the accepted ranges (e.g. 6𝜎), which is more economical and practical than the worst-case method.

To study the underlying mechanisms in 3D, Bourdet and Clément introduced the concept of the Small Displacement Torsor (SDT) to evaluate the defects on a part surface [START_REF] Bourdet | Controlling a complex surface with a 3 axis measuring machine[END_REF] and conduct tolerance analysis. Later, more ideal surface-based tolerance analysis methods were developed [START_REF] Requicha | Toward a Theory of Geometric Tolerancing[END_REF][START_REF] Gao | Generalized 3-D tolerance analysis of mechanical assemblies with small kinematic adjustments[END_REF][START_REF] Clément | Theory and Practice of 3-D Tolerancing for Assembly[END_REF][START_REF] Desrochers | Application of a Unified Jacobian Torsor Model for Tolerance Analysis[END_REF][START_REF] Teissandier | Operations on polytopes: application to tolerance analysis[END_REF][START_REF] Arroyave-Tobón | Applying screw theory for summing sets of constraints in geometric tolerancing[END_REF][START_REF] Davidson | A New Mathematical Model for Geometric Tolerances as Applied to Round Faces[END_REF][START_REF] Zou | A gap-based approach to capture fitting conditions for mechanical assembly[END_REF]. They are based on the assumption that the influence of form defects on assembly deviation is negligible.

However, with the increasing demands for a more comprehensive management of product geometrical deviation, current ideal surface-based methods cannot fulfill the requirements. Form defects need to be considered since more and more mechanical parts require high geometrical precision. Moreover, tolerance management should be extended to whole product life cycle management, and it should be considered as early as possible [START_REF] Ballu | A new Design Method based on Functions and Tolerance Specifications for Product Modelling[END_REF].

In this movement, standards are evolving towards a higher degree of description of the specification to improve communication between various stages and organizations. In the series of Geometrical Product Specifications (GPS) standards [START_REF]Geometric Product Specification-General concepts-Part 1: Model for geometrical specification and verification n[END_REF][START_REF]Geometric Product Specification-General concepts-Part 2: Basic tenets, specifications, operators, uncertainties and ambiguities[END_REF][START_REF]Geometrical product specifications (GPS) --Matrix model n[END_REF], models and terms to describe product precision requirement are standardized. These symbolic specifications meet their limits in avoiding ambiguities and generalizing concepts for new functional and manufacturing needs [START_REF] Charpentier | A scientific point of view of a simple industrial tolerancing process[END_REF][START_REF] Dantan | Geometrical product specifications -model for product life cycle[END_REF][START_REF] Mathieu | A Model for a Coherent and Complete Tolerancing Process[END_REF][START_REF] Dantan | Geometrical variations management for additive manufactured product[END_REF]. To generalize and to communicate geometrical specifications without ambiguities, Ballu and Mathieu proposed a language (GeoSpelling) and a conceptual model (Skin Model). These were selected as the basis of standard IS0 17450.

GeoSpelling [START_REF]Geometric Product Specification-General concepts-Part 1: Model for geometrical specification and verification n[END_REF][START_REF] Mathieu | GEOSPELLING: a common language for Specification and Verification to express Method Uncertainty[END_REF] is a language based on limited operations. The design of GeoSpelling considered semantic coherence, the demand for digitalization and automation, readability for both humans and machines, and the learning cost for engineers. In the work of Ballu et al. [START_REF] Ballu | Formal Language for GeoSpelling[END_REF], a formal language for GeoSpelling is proposed. Basic instructions like logical expressions and loops are given with examples. Implementation of GeoSpelling on gear parts was researched by Zhang et al. [START_REF] Zhang | Geometric Product Specification of Gears: The GeoSpelling Perspective[END_REF]. In their work, the geometrical characteristics for gears are deduced and expressed using the GeoSpelling approach.

The aim of the Skin Model [START_REF]Geometric Product Specification-General concepts-Part 1: Model for geometrical specification and verification n[END_REF][START_REF] Ballu | Analysis of dimensional and geometrical specifications: standards and models[END_REF] is to represent the geometry of actual parts. In CAD systems, the product is designed and viewed in the nominal condition, which satisfies the design constraints and contains no form defects. In contrast, the skin model does contain manufacturing defects. It is an infinite model defined as a "model of the physical interface of the work piece with its environment" [START_REF]Geometric Product Specification-General concepts-Part 1: Model for geometrical specification and verification n[END_REF]. Therefore it is a link between the real physical part and the virtual nominal part. In the work of Anwer et al. [START_REF] Anwer | The skin model, a comprehensive geometric model for engineering design[END_REF], the skin model concept was analyzed at different levels. The skin model shape, which is the discrete instance of the skin model, is also proposed for computation and visualization requirements [START_REF] Zhang | A discrete geometry framework for geometrical product specifications[END_REF][START_REF] Schleich | Skin Model Shapes: A new paradigm shift for geometric variations modelling in mechanical engineering[END_REF][START_REF] Yan | Generation of consistent skin model shape based on FEA method[END_REF].

At this stage, the generation of skin model shape and related activities for tolerancing does not meet the needs of industry and much more research must be carried out. The objective of this work is to complement earlier research and further develop new approaches for a better response to demands. In Chapter 3, consistent skin model shapes are generated. Form defects on surfaces are simulated using the methods introduced in Chapter 2. The Finite Element Analysis (FEA) based method and mesh smoothing method are proposed to generate the complete skin model shape.

These skin model shapes containing manufacturing defects will be used for simulation and visualization in later works.

In Chapter 4, assembly simulation is conducted. Differences between assembly using the nominal model, ideal surface-based deviation model, and skin model shapes are pointed out. The importance of boundary conditions for assembly simulation are emphasized. A simulation method based on Linear Complementarity Condition (LCC) and Small Displacement Torsor (SDT) is developed. Examples and analysis are given at the end.

In Chapter 5, the virtual metrology method is developed. Deviation on single skin model shape or skin model shape assembly can be evaluated. Meanwhile, the method is customized and integrated into the development of an online Virtual Laboratory (VL) for e-learning.

Form Defect Simulation

In this chapter, manufacturing defect simulation methods are reviewed and classified.

Examples are given based on a unique model to show their differences and potential applications.

Criteria are proposed to help researchers and engineers choose suitable simulation methods.

These simulation methods will be used to generate skin model shapes.

Introduction

Concurrent Engineering is a solution to introduce downstream knowledge, such as manufacturing and maintenance, into product design [START_REF] Yassine | Complex Concurrent Engineering and the Design Structure Matrix Method[END_REF]. In early design stages, modification of geometry and structure is less expensive as no real product has yet been manufactured. Many Computer-Aided Engineering (CAE) tools are used during the design stage to anticipate possible product malfunctions. Among these are Computer Aided Tolerancing (CAT) systems, which are developed to analyze and control geometrical problems in the design stage [START_REF] Clément | Theory and Practice of 3-D Tolerancing for Assembly[END_REF][START_REF] Teissandier | Operations on polytopes: application to tolerance analysis[END_REF][START_REF] Geis | Use of Vectorial Tolerances for Direct Representation and Analysis in CAD-systems[END_REF][START_REF] Desrochers | A dimensioning and tolerancing assistance model for CAD/CAM systems[END_REF].

Compared with the fact that geometrical quality is influenced by various stages of manufacturing and design, the current methods used for engineering and research are usually based on simple variation of the ideal features [START_REF] Schleich | Comparative Study on Tolerance Analysis Approaches[END_REF], which are not able to represent the real situation of the product or cover the whole product lifecycle. Meanwhile, the ideal feature-based expression of tolerance specifications leads to ambiguity [START_REF] Ballu | Formal Language for GeoSpelling[END_REF] and may cause problems when integrated with the Product Lifecycle Management (PLM) system.

In recent research on CAT, more and more studies have tried to introduce the idea of skin model shape and non-nominal geometry [START_REF] Zhang | Discrete shape modeling for skin model representation[END_REF][START_REF] Franciosa | Simulation of variational compliant assemblies with shape errors based on morphing mesh approach[END_REF][START_REF] Wagersten | Using Morphing Techniques in Early Variation Analysis[END_REF][START_REF] Formosa | Modal expression of form defects[END_REF][START_REF] Huang | Mode-based Decomposition of Part Form Error by Discrete-Cosine-Transform with Implementation to Assembly and Stamping System with Compliant Parts[END_REF][START_REF] Wilma | Manufacturing Signature for Tolerance Analysis[END_REF][START_REF] Armillotta | Tolerance Analysis Considering form Errors in Planar Datum Features[END_REF][START_REF] Bruyère | Statistical tolerance analysis of bevel gear by tooth contact analysis and Monte Carlo simulation[END_REF]. Several new methods have been proposed to simulate manufactured parts or reconstruct measurement data. General problems have been raised, such as how to simulate manufactured surfaces with or without detailed manufacturing information, and how to integrate non-nominal skin model shapes into current CAT systems, etc. The first step in solving these problems is to generate skin model shapes that reflect real parts.

In this study, features of the part are treated independently, emphasizing feature segmentation and combination processes [START_REF] Yan | Toward an Automatic Generation of Part Models with Form Error[END_REF]. Simulation of the skin model shape is divided into three steps, as described below:

 Segmentation. To be able to treat each feature independently and simulate deviations with different precision requirements, segmentation of the feature is conducted in advance.

 Manufacturing deviation generation. Geometric manufacturing defects in the skin model shape are simulated. They are saved as deviation data for each vertex. The segmentation and combination are to take the influences of multi-stage manufacturing process into consideration. On an actual workpiece, particularly a machined workpiece, there are independencies between isolated features or groups of fixtures due to specific machining conditions and clamping positions; different clamping positions introduce orientation and position errors. To be more realistic, the simulation is applied independently to each feature of the part or to each group of features; this enables location errors to be introduced according to the manufacturing process and part setup.

The combination method in Figure 2-1(d) will be introduced in Chapter 3. In this chapter, we focus on the form defect generation, which corresponds to Figure 2-1(c). The main purpose of this chapter is to help researchers to choose suitable methods to generate realistic skin model shapes which could be used for research and engineering. To this end, a classification and comparison of different existing simulation methods were conducted.

First, we give an overview of the generation methods, then we classify the different methods.

Based on this classification, non-mode based methods and mode based methods are introduced respectively. A comparison of the different methods is provided in a later section. Finally, conclusions and further work are described. 2.2 Form Error Simulation Methods

Overview of Different Methods

The one obvious difference between skin model shape and nominal model is that the former contains more detailed shape information than the latter. Thus, to generate skin model shapes, the related form error simulation methods should be considered in detail.

To simulate manufacturing deviations on the part, there are two quite different strategies [START_REF] Yan | Toward an Automatic Generation of Part Models with Form Error[END_REF].

The first is manufacturing process oriented, which simulates the detailed manufacturing process and generates realistic parts. These methods are often considered as CAE methods, like Finite Element Analysis (FEA) and Virtual Machining (VM). The simulation result covers not only the part model with manufacturing deviations, but also the process parameters which could be used to improve the manufacturing processes.

The other kind of strategy is geometry oriented methods. Here, the detailed manufacturing process and factors, like the path of tools and vibrations, are not considered. The aim is only to generate a shape with form defect which is similar to the manufactured parts. In this chapter, we will focus on the geometry oriented methods.

As form error is an important factor that influences tolerancing, different authors have produced studies with geometry oriented methods. In [START_REF] Zhang | A discrete geometry framework for geometrical product specifications[END_REF], both 1D and 3D random noise are used to generate the skin model shapes with form error. Because the variations in vertex coordinates are generated independently, distances between two connected vertices may be dramatic, which causes non-smooth results. The mesh morphing approach, a method that edits the meshed model directly, is used to simulate form error. Both Franciosa et al. [START_REF] Franciosa | Simulation of variational compliant assemblies with shape errors based on morphing mesh approach[END_REF][START_REF] Das | Physics-driven Shape Variation Modelling at Early Design Stage[END_REF] and Wagersten et al. [START_REF] Wagersten | Using Morphing Techniques in Early Variation Analysis[END_REF] used the morphing method to generate form error in compliant assembly analysis, and Zhang et al. [START_REF] Zhang | Discrete shape modeling for skin model representation[END_REF] used it to generate certain form errors (e.g. paraboloid, cone and ellipsoid). Mode based methods [START_REF] Formosa | Modal expression of form defects[END_REF][START_REF] Huang | Mode-based Decomposition of Part Form Error by Discrete-Cosine-Transform with Implementation to Assembly and Stamping System with Compliant Parts[END_REF][START_REF] Wilma | Manufacturing Signature for Tolerance Analysis[END_REF][START_REF] Kingslake | Basic wavefront aberration theory for optical metrology[END_REF][START_REF] Cao | Geometrical Simulation of Multiscale Toleranced Surface With Consideration of the Tolerancing Principle[END_REF][START_REF] Ganti | Generalized fractal analysis and its applications to engineering surfaces[END_REF][START_REF] Hasegawa | Calculation of the fractal dimensions of machined surface profiles[END_REF][START_REF] Nguyen | Modeling cutter tilt and cutter-spindle stiffness for machine condition monitoring in face milling using high-definition surface metrology[END_REF] generate surfaces based on several typical modes, which come from modal decomposition of measurement data or other modal analysis methods, like FEA. The random field method is used to simulate both irregular forms and random parameters in structural mechanics [START_REF] Mollon | 3D generation of realistic granular samples based on random fields theory and Fourier shape descriptors[END_REF], and it has also been used to simulate form error [START_REF] Schleich | Skin Model Shapes: A new paradigm shift for geometric variations modelling in mechanical engineering[END_REF]. In addition, statistical shape analysis based on Principal Component Analysis (PCA) is proposed for modeling 3D surfaces [START_REF] Pacella | Different Formulations of Principal Component Analysis for 3D Profiles and Surfaces Modeling[END_REF] and generating skin model shapes [START_REF] Zhang | Discrete shape modeling for skin model representation[END_REF] from training data sets.

Considering that there are so many methods, and that sometimes they share basic ideas, we classified these different methods into three categories. The classification is based on the principle of the methods and simulation results, they are: random noise, mesh morphing, and mode based methods. The mode based methods were subdivided into trigonometric function methods and spectral methods. Detailed methods and their categories can be seen in Figure 2-2, theories and applications will be discussed in the following sections.

Model Used for Simulation

To evaluate and compare different simulation methods, here we take an example and apply these methods to the same model. The nominal model we use is a rectangular plane. We assume that the generated form deviation is along the normal direction of vertices. The discrete model of the rectangular surface can be seen in Figure 23. The generation of skin model shape is concluded as the combination of the nominal model and manufacturing defects 𝝃𝝐:

𝑆𝐾 = 𝑆𝐾 ̅̅̅̅ + 𝝃𝝐 (1)
where 𝝃 is the 𝑛 × 𝑚 matrix. 𝑛 is the number of variables, which is equal to the number of vertices. 𝑚 is the number of column vectors that contain manufacturing defects. For each column in 𝝃, the defect could be simulated by random noise methods, mesh morphing methods or mode based methods. The coefficients or deviation value in 𝝐 could be generated based on former experience, or on analysis of the measurement data when manufactured parts are available.

Because we simulate only the form error and we do not decompose the measurement data, several typical modes are carefully chosen for mode based methods. Since the principle of different methods could be different, it is hard to set uniform control parameters or coefficients and compare the simulation result. Thus, we considered them as random variables, and the results aimed at showing the potential applications of these methods. The simulation results are presented after the illustration of each method. For methods which have specific applications, we also provide additional simulation examples.

Non-Mode Based Methods

Random Noise Method

Simulation with random noise to represent the form error on skin model shapes was conducted by Zhang et al. [START_REF] Zhang | A discrete geometry framework for geometrical product specifications[END_REF]. In their method, both 1D and 3D Gaussian methods were introduced to simulate the form error.

In the 1D Gaussian random noise method, the Probability Density Function (PDF) has the following form:

𝑓(𝑣) = 1 𝜎√2𝜋 𝑒 - (𝑥-𝜇) 2 2σ 2 (2)
where 𝜎 is the variance, 𝜇 is the mean value, and 𝑥 is the Gaussian variable.

The idea of the 1D Gaussian random noise method is illustrated in Figure 2-4(a). The deviation of vertex on the model is along its local normal direction. The 3D Gaussian method is the extension of the 1D Gaussian method, and contains deviations in three directions, as can be seen in Figure 2-4(b). We use the 1D Gaussian method for simulation, every vertex is given a random deviation along its normal direction, which follows the Gaussian distribution with mean value 𝜇 = 0 and standard deviation 𝜎 = 0.01. The simulation result is shown in Figure 4. To visualize the shape after simulation, the form error is amplified. The color change indicates the different values of deviation. As can be seen from the result, regions of similar color tend to be limited in a very small area and the size of these small areas is related to the size of the mesh. There is no relation between neighboring vertices and the surface is very chaotic. Compared with the size of feature and the number of vertices, this method is not sufficient to represent various kinds of form errors on the whole model. Nevertheless, because of its local property, it is suitable for simulating random noise and complementing other form errors.

Mesh Morphing Methods

Random Shape Morphing Using Shape Function

With the wide application of CAE tools, the mesh-based model is used in many applications, and modifying methods that are applied to the mesh directly have been developed to accelerate the design and simulation speed [START_REF] Raffin | Curvilinear displacement of free-form-based deformation[END_REF][START_REF] Staten | A comparison of mesh morphing methods for 3D shape optimization[END_REF]. Morphing is one of the methods that works on the meshed models directly. By specifying the control points and their influence parameters (e.g. influence range, shape functions), users can modify the model efficiently.

Wagersten et al. [START_REF] Wagersten | Using Morphing Techniques in Early Variation Analysis[END_REF] applied the morphing method to simulate the assembly of non-rigid sheet metals. Their morphing process is conducted by Altair HyperWorks. Meanwhile, Franciosa et al. [START_REF] Franciosa | Simulation of variational compliant assemblies with shape errors based on morphing mesh approach[END_REF] introduced morphing and finite element analysis to simulate the classic welding assembly process, which consists of Place, Clamp, Fasten and Release (PCFR). The following details of the morphing method are based on the work of Franciosa et al. [START_REF] Franciosa | Simulation of variational compliant assemblies with shape errors based on morphing mesh approach[END_REF].

The morphing method is illustrated in Figure 2-6. In the simulation of the form error, deviation is given to the 𝑖-th control point 𝑣 𝑐,𝑖 at first. The subscript index "𝑐" means control point. The influence range of control point 𝑣 𝑐,𝑖 is called the influence hull 𝑙 𝑖 .

Considering a vertex 𝑣 𝑗 inside the influence hull 𝑙 𝑖 , its deviation is controlled by a shape function 𝑓. Different shape functions may be used to define the morphing shape. Some of these functions guarantee smooth shapes and zero slope at the boundary of influence hull. As an example, the shape function in Equation ( 3) is used by Franciosa et al. [START_REF] Franciosa | Simulation of variational compliant assemblies with shape errors based on morphing mesh approach[END_REF]:

𝜉 𝑐𝑖,𝑗 = 𝑓 ( ‖𝑣 𝑗 -𝑣 𝑐,𝑖 ‖ 𝑙 𝑖 ) = 𝑓(𝑑) = { 1 -3𝑑 2 + 2𝑑 3 𝑑 < 1 0 𝑑 ≥ 1 (3)
where ‖𝑣 𝑗 -𝑣 𝑐,𝑖 ‖ is the distance between 𝑣 𝑗 and 𝑣 𝑐,𝑖 . The calculation result 𝜉 𝑐,𝑖,𝑗 indicates the influence of control point 𝑣 𝑐,𝑖 on point 𝑣 𝑗 . Assume the control point has deviation 𝜖 𝑐,𝑖 , then the deviation on point 𝑣 𝑗 is calculated as:

𝜖 𝑗 = 𝜉 𝑐,𝑖,𝑗 • 𝜖 𝑐,𝑖 (4)
When applying the morphing process, there may be several control points, and the deviation values for all the points on the mesh should be calculated. In this case, Equation ( 4) is in matrix form, the same as Equation (1). The influence hull mentioned above is a sphere with radius 𝑙 𝑖 . However, the shape of the influence hull may be different. In [START_REF] Franciosa | Simulation of variational compliant assemblies with shape errors based on morphing mesh approach[END_REF], an ellipsoid influence hull is given as an example (Figure

𝑑 is the variable of the shape function in Equation (3), 𝑣 𝑗,𝑥 , 𝑣 𝑗,𝑦 and 𝑣 𝑗,𝑧 are coordinates of 𝑣 𝑗 along the 𝑥, 𝑦, 𝑧 directions. 𝑣 𝑐,𝑖,𝑥 , 𝑣 𝑐,𝑖,𝑦 and 𝑣 𝑐,𝑖,𝑧 are coordinates of 𝑣 𝑐,𝑖 . 𝑙 𝑖,𝑎 , 𝑙 𝑖,𝑏 , 𝑙 𝑖,𝑐 are the length of semi-principal axes of the ellipsoid. To generate form errors which are more realistic, different kinds of influence hulls and shape function may be introduced.

For example, we set the number of control points at a constant 4, and they are chosen from all the vertices randomly. Their deviations along the normal direction follow normal distribution. To simplify, the influence hull used here is a sphere which has been introduced before, and its radius is a random value ranging from 50𝑚𝑚 to 100𝑚𝑚. The shape function in Equation ( 3) is used to control the deviation shape. The result is shown in Figure 234567.

This is a very flexible method, as the parameters in the simulation have specific meaning and can be controlled. A given form error could be generated by adjusting the parameters, like the position of control points and their influence hull. Meanwhile, the simulation result could be more random if we set all the parameters randomly.

Second order shape morphing

In manufacturing, some errors on the parts tend to show similar characteristics in a batch of products. These errors are usually caused by kinematic errors by machine tools, and they influence every manufactured part in a similar way and vary slowly [START_REF] Kiridena | Kinematic modeling of quasistatic errors of three-axis machining centers[END_REF][START_REF] Uddin | Prediction and compensation of machining geometric errors of five-axis machining centers with kinematic errors[END_REF]. Because these errors share some common characteristics, like shape, pattern or amplitude, compensations are usually introduced based on the study of these.

Second order shapes have been widely used in surface fitting, simulation and reconstruction [START_REF] Bolle | On three-dimensional surface reconstruction methods[END_REF]. To simulate kinematic errors on a plane and cylinder, Zhang et al. [START_REF] Zhang | Discrete shape modeling for skin model representation[END_REF] used typical second order shapes to express them. In their work, the deviation 𝜖 along the vertex normal direction could be expressed by a second order function 𝜖 𝑖 = 𝑓(𝑥 𝑖 , 𝑦 𝑖 , 𝑧 𝑖 ), where 𝑥 𝑖 , 𝑦 𝑖 , 𝑧 𝑖 are the coordinates of vertices in local coordinate associated to the considered surface. Using this method, they constructed paraboloid, cone, sphere, cylinder and ellipsoid for plane features, and taper, concave, convex and banana for cylinders. Figure 2345678shows some examples of simulations on a cylinder.

Complex shapes could be generated by the superposition of several different shapes, and also of random noises. The translation and rotation of the whole surface were also introduced, based on rigid body movement, to move all the vertices of the feature together.

To simulate form errors on our example, two kinds of shape (paraboloid and cylinder) are simulated, as can be seen in Figure 2-9. For Figure 2-9(a), deviation of vertices along the 𝑧 direction is calculated from Equation 6, and the methods are similar for other second order shapes:

𝜖 𝑖 = 𝑣 𝑖,𝑥 2 𝑎 2 + 𝑣 𝑖,𝑦 2 𝑏 2 (6)
where 𝑎, 𝑏 are constants related to the curvature of the paraboloid. Even kinematic error shows some rules or patterns; these patterns are still changing slowly and parts are different from each other. Thus, to simulate a batch of parts with kinematic error, variations should be introduced to the parameters of second order shape functions.

Mode Based Methods

Trigonometric Function Modes

To model manufacturing deviations on a plane surface, trigonometric functions are simple harmonic terms that could be used, and they are easy to modify by varying the coefficients. In the work of Wilma et al. [START_REF] Wilma | Manufacturing Signature for Tolerance Analysis[END_REF], a model based on trigonometric terms is used to simulate the manufacturing deviations on a circle model. The Weiestrass-Mandelbrot fractal function, which contains trigonometric terms, is used to simulate the manufactured surfaces [START_REF] Armillotta | Tolerance Analysis Considering form Errors in Planar Datum Features[END_REF][START_REF] Ganti | Generalized fractal analysis and its applications to engineering surfaces[END_REF][START_REF] Hasegawa | Calculation of the fractal dimensions of machined surface profiles[END_REF]. There are also other models that are based on trigonometric functions [START_REF] Homri | Tolerance analysis --Form defects modeling and simulation by modal decomposition and optimization[END_REF][START_REF] Qiao | Approach to the Deviation Representation of Non-ideal Cylindrical Surfaces Based on the Curvilinear Coordinate System[END_REF][START_REF] Song | Shape deviation modeling for fused deposition modeling processes[END_REF]. In the following, we describe in detail the most used methods: Zernike polynomials and the Discrete Cosine Transform (DCT).

Zernike Polynomials

Zernike polynomials are used to decompose the form errors of optical parts since they form an orthogonal basis that could represent typical defects in these parts [START_REF] Kingslake | Basic wavefront aberration theory for optical metrology[END_REF][START_REF] Zernike | Phase contrast, a new method for the microscopic observation of transparent objects[END_REF]. Zernike polynomials represent the typical defects in a unit disc by polar coordinate (𝜌, 𝜃), and they take the form of even and odd respectively [START_REF] Lê | Influence des défauts de forme sur le comportement des liaisons: étude expérimentale et théorique[END_REF]:

{ 𝑍 𝑛 𝑚 (𝜌, 𝜃) = 𝑅 𝑛 𝑚 (𝜌)cos (𝑚𝜃) 𝑒𝑣𝑒𝑛 𝑍 𝑛 -𝑚 (𝜌, 𝜃) = 𝑅 𝑛 𝑚 (𝜌)sin (𝑚𝜃) 𝑜𝑑𝑑 (7)
where 𝜌 is the radius ranging from 0 to 1, 𝜃 is the angle from 0 to 2𝜋. 𝑛 is the number of highest order polynomials and 𝑚 is the natural numbers which index the frequency of components ( 𝑚 ≤ 𝑛 ). The definition of 𝑅 𝑛 𝑚 (𝜌) is: As can be seen from Equation 7and Figure 2-10, the simulation based on Zernike polynomials fits the disc shape parts better, like a lens, compared with our usage, which usually has more general shapes. Thus it is not applied in our comparison.

𝑅 𝑛 𝑚 (𝜌) = { ∑ (-1) 𝐾 (𝑛-𝑘)! 𝐾!( 𝑛+𝑚 2 -𝑘)!( 𝑛-𝑚 2 -𝑘)! (𝑛-𝑚) 2 ⁄ 𝑘=0 𝜌 (𝑛-2𝑘) 𝑛 -𝑚 𝑖𝑠 𝑒𝑣𝑒𝑛 0 𝑛 -𝑚 𝑖𝑠 𝑜𝑑𝑑 (8)

Discrete Cosine Transform

DCT is a signal processing method. Because of the availability of fast algorithms and its energy compaction property, it has been widely used in image compression. Huang et al. [START_REF] Huang | Mode-based Decomposition of Part Form Error by Discrete-Cosine-Transform with Implementation to Assembly and Stamping System with Compliant Parts[END_REF] used forward and inverse transformation to decompose and reconstruct form error based on measurement data.

In the case of a two-dimensional (𝑥, 𝑦 direction) rectangle plane surface, form error is expressed by deviation along the third direction (𝑧 direction). To conduct decomposition and reconstruction, 𝑘 sample points are taken equally along the 𝑥 and 𝑦 axes, thus there are 𝑘 2 sample points in total.

Let 𝑛 , 𝑚 indicate the numbers of sample points along the 𝑥 , 𝑦 directions, which are independent with the actual coordinates, and 𝑢, 𝑣 be the frequency along the 𝑥, 𝑦 axes, the modes of DCT method could then be expressed as:

𝛿(𝑛, 𝑚, 𝑢, 𝑣) = { 1 𝑘 𝑢 = 𝑣 = 0 2 𝑘 𝑐𝑜𝑠 ( 2𝜋𝑛𝑢 2𝑘-1 ) 𝑐𝑜𝑠 ( 2𝜋𝑚𝑣 2𝑘-1 ) 𝑜𝑡ℎ𝑒𝑟𝑠 (9)
where the value of 𝛿(𝑛, 𝑚, 𝑢, 𝑣) indicates the deviation of sample point (𝑛, 𝑚) under the mode (𝑢, 𝑣). (0,1) The advantage of the DCT method is that the explicit form makes it easy to control the shape, and any shape of form error could be reconstructed. However, it is usually limited to plane features; to apply this method to other general shapes, a projection process may be required before simulation [START_REF] Wilma | Manufacturing Signature for Tolerance Analysis[END_REF][START_REF] Homri | Tolerance analysis --Form defects modeling and simulation by modal decomposition and optimization[END_REF][START_REF] Qiao | Approach to the Deviation Representation of Non-ideal Cylindrical Surfaces Based on the Curvilinear Coordinate System[END_REF][START_REF] Song | Shape deviation modeling for fused deposition modeling processes[END_REF].

(1,0) (1,1) (1,2) (2,1) (2,2)

Spectral Methods

Simulation Based on Spectral Graph Theory

Considering the geometry we handle is discretized by triangle mesh, it is easy to link the topology of the mesh to a graph. A graph is defined as 𝐺 = (𝑉, 𝐸), where 𝑉 is the set of vertices and 𝐸 is the set of edges [START_REF] Bondy | Graph theory with applications[END_REF]. Thus for a mesh representing geometry, we can easily compare the vertices and edges of the mesh with the vertices and edges of the corresponding graph.

Based on the topology of the graph, a matrix could be associated with it. Two widely used matrices associated to a graph are the adjacency matrix 𝐴 and the Laplacian matrix 𝐿 [START_REF] Chung | Spectral graph theory[END_REF]. These matrices could be considered as an operator or a quadratic form. In [START_REF] Ohtake | Mesh regularization and adaptive smoothing[END_REF], the authors smooth the mesh by combining the Laplacian flow with the function of mean curvature, while in [START_REF] Vallet | Spectral geometry processing with manifold harmonics[END_REF] the Laplacian Beltrami operator is used for filtering and simplifying the geometry.

The analysis related to the eigenvalue and eigenvector of the Laplacian matrix is called the spectral graph theory. The matrices and their eigenvectors could show some interesting properties of the original graph, and have been widely applied to computer graphic research. By improving the Laplacian matrix, the spectral based mesh segmentation has gained lots of research attention due to its simplicity and efficiency [START_REF] Wang | Spectral 3D mesh segmentation with a novel single segmentation field[END_REF][START_REF] Liu | Mesh segmentation via spectral embedding and contour analysis[END_REF].

Let 𝑖 and 𝑗 be two vertices on the graph (and also the corresponding geometry) which contains 𝑛 vertices, the unweighted adjacency matrix of a graph is defined as:

𝐴 = { 1 𝑖𝑓(𝑖, 𝑗) ∈ 𝐸 0 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 (10)
A diagonal matrix 𝐷 is defined so that its diagonal entry 𝐷(𝑖, 𝑖) is the degree of the vertex 𝑖.

For an unweighted graph, the degree of a vertex is equal to the number of edges connected to it.

Based on the two matrices, the graph Laplacian is [START_REF] Zhang | Spectral Mesh Processing-COMPUTER GRAPHICS forum[END_REF]:

𝐿 = 𝐷 -𝐴 (11)
The graph Laplacian is a real symmetrical matrix. Based on the Spectral Theorem [START_REF] Lay | Linear algebra and its applications[END_REF], its eigenvectors corresponding to different eigenvalues are orthonormal. The eigenvalue and eigenvector of a graph Laplacian reveals some deeper characteristics of the graph. To understand certain characteristics, different geometry related weights could be introduced to adjacency matrix 𝐴, like distance and curvature.

One important characteristic is that the eigenvector of a Laplacian matrix contains oscillation, and tends to partition the vertices. This can be illustrated by a simple example of line mesh. The Using this property, modes can be generated, as in the DCT method. By modifying the weights in adjacency matrix 𝐴, the geometry property is introduced to the Laplacian matrix. Here we define the weight as:

𝐴(𝑖, 𝑗) = e (-𝑑 𝑖,𝑗 ′ ) , 𝑖𝑓(𝑖, 𝑗) ∈ 𝐸 (12)
where 𝑑 𝑖,𝑗 ′ = 4

𝑑 𝑖,𝑗 𝑑 𝑚𝑎𝑥
is the scaled mesh edge length, 𝑑 𝑖,𝑗 is the edge length and 𝑑 𝑚𝑎𝑥 is the maximum value of 𝑑 𝑖,𝑗 . The diagonal entry in 𝐷 is:

𝐷(𝑖, 𝑖) = ∑ 𝐴(𝑖, 𝑘) 𝑛 𝑘=1 (13)
In our 3D case with the discrete plane surface, typical eigenvectors are shown in

(c) (d) (a) (b) (c) (d) (e) (f)
It also should be noticed that even though they have something in common, their theoretical origins and backgrounds are quite different. We will start from their theoretical origins, then conclude and compare their numerical application from the perspective of spectral graph theory.

Natural Vibration Mode Analysis

To characterize the form errors, Formosa et al. [START_REF] Formosa | Modal expression of form defects[END_REF][START_REF] Ledoux | Integrating form defects of mechanical joints into the tolerance studies[END_REF][START_REF] Goic | Multi scale modal decomposition of primary form, waviness and roughness of surfaces[END_REF] proposed a method based on structural modal analysis using FEA. Orthogonal modes are generated, and manufacturing deviations are simulated by linear combination of these modes. Measurement data can also be decomposed in an inverse process.

Generally, a linear dynamic system with mass matrix 𝑀 and stiffness matrix 𝐾 can be written as [START_REF] Formosa | Modal expression of form defects[END_REF]:

𝑀𝒒̈+ 𝐾𝒒 = 0 (14)
𝒒 is the displacement vector. The modes of the system can be written in a general form as:

𝒒 𝒊 = 𝑸 𝒊 • 𝑐𝑜𝑠(𝜔 𝑖 • 𝑡) (15)
where 𝒒 𝒊 indicates the 𝑖-th mode, 𝑸 𝒊 is the amplitude vector and 𝜔 𝑖 is the pulsation. The characteristic equation of the problem is deduced by introducing Equation [START_REF]Geometric Product Specification-General concepts-Part 2: Basic tenets, specifications, operators, uncertainties and ambiguities[END_REF] to Equation ( 14):

(𝐾 -𝜔 2 𝑀)𝒒 = 0 (16)
Assume each edge in the mesh model is a 3D bar element, the material is isotropic, and the mass matrix [START_REF] Zeng | Fundamentals of Finite Element Analysis[END_REF] is:

𝑀 = 𝑚𝑰 (17)
Mass 𝑚 is decided by material properties and 𝑰 is the identity matrix. Let 𝑚 = 1, 𝑘 𝑖,𝑗 is the stiffness between node 𝑖 and 𝑗, we can compare the stiffness matrix to the Laplacian in Equation ( 11) and get:

{ 𝐴(𝑖, 𝑗) = 𝑘 𝑖,𝑗 , 𝑖𝑓(𝑖, 𝑗) ∈ 𝐸 𝐷(𝑖, 𝑖) = ∑ 𝐴(𝑖, 𝑘) 𝑛,𝑘≠𝑖 𝑘=1 (18)
In application [START_REF] Ledoux | Integrating form defects of mechanical joints into the tolerance studies[END_REF], the CAD model is first meshed to shell elements by finite element software, related material parameters are given to guarantee a reasonable simulation result. The mode bases are then generated by calculating the vibration modes.

We selected several typical modes and specified the coefficients randomly. The modes we (

) (3) (7) (8) (9) 1

Random Field Method

As defined in [START_REF] Adler | Random fields and geometry[END_REF], a random field is a continuous function of random process with multidimensional variables, and the relationships between variables are described by autocorrelations.

Since the random field is a continuous function, to apply it to discrete models and facilitate their manipulation, the random field is discretized [START_REF] Li | Optimal Discretization of Random Fields[END_REF].

The random field method is introduced to CAE simulations to simulate the variation of properties (e.g. loads, material, geometry, manufacturing process, etc.). Stochastic FEA using the random field method is conducted to evaluate structure reliability [START_REF] Sudret | Stochastic finite element methods and reliability: a state-of-theart report[END_REF]. In [START_REF] Schleich | Skin Model Shapes: A new paradigm shift for geometric variations modelling in mechanical engineering[END_REF], the random field method based on series expansion is used to simulate form error on skin model shapes. Let 𝜹 be the vector of vertex deviation along the normal direction, then:

𝜹 = 𝝁 + 𝝃𝝀 𝟏 𝟐 𝝈𝝐 (19)
𝝁 and 𝝈 are the mean value vector and standard deviation matrix for the random field respectively. 𝝀 and 𝝃 are matrices formed by the 𝑘 largest eigenvalues and eigenvectors of autocovariance matrix 𝐶𝑜𝑣. The coefficient vector 𝝐 is chosen randomly, thus a skin model shape with random form error is generated.

The element of auto-covariance matrix 𝐶𝑜𝑣 can be expressed as:

𝐶𝑜𝑣 𝑖,𝑗 = 𝜌 𝑖,𝑗 = 𝑒 (- ‖𝑣 𝑖 -𝑣 𝑗 ‖ 2 𝑙 𝜌 2 ) (20)
where 𝜌 𝑖,𝑗 is the squared exponential autocorrelation function between two vertices. 𝑙 𝜌 is the correlation length which determines the influence range of the vertex. In the prediction stage, 𝑙 𝜌 could be specified from experience, while in the observation stage, it could be calculated from measurement or simulation data.

The random field method could also be regarded as a spectral method, as the eigenvectors and eigenvalues of auto-covariance matrix 𝐶𝑜𝑣 are calculated. Compared to the Laplacian matrix, which depends on the topology of the mesh, the random field method established the relation between all the vertices, and thus is a complete graph. The auto-covariance matrix could be rewritten with matrices 𝐴 and 𝐷 like the Laplacian matrix:

{ 𝐴(𝑖, 𝑗) = 𝜌 𝑖,𝑗 , 𝑖, 𝑗 ∈ 𝑉 𝐷(𝑖, 𝑖) = 1 (21)
Eigenvectors of the auto-covariance matrix show similar properties to other spectral methods.

In From the figure, we can see that the deformation on each surface is changing gradually. With the increasing influence of length 𝑙 𝜌 , the average area of single deformation range on the surface increases at the same time.

(

) (2) (3) (4) 1

Principal Component Analysis

Compared with methods that generate form error that is totally dependent on random variables, in some situations process simulation or measurement data are available and more realistic skin model shapes can be generated.

Statistical Shape Analysis (SSA) is different from the methods described above, and is used to analyze a set of training data by statistical methods. SSA based on PCA has been widely used for image processing in medical research [START_REF] Aguirre | Statistical shape analysis via principal factor analysis[END_REF]. Zhang et al. [START_REF] Zhang | Discrete shape modeling for skin model representation[END_REF] used this method to generate skin model shapes with form defects based on simulation data.

PCA is a linear multivariate analysis method. By calculating the eigenvector and eigenvalue of covariance matrix, PCA aims to find a projection that reduces the data dimension and reveals inner properties. The projection could be realized by changing the original data coordinates to new ones, which correspond to the eigenvector direction [START_REF] Lay | Linear algebra and its applications[END_REF]. The eigenvectors represent the principal direction, and the eigenvalues indicate the variance of data along their corresponding eigenvectors.

In the method developed by Zhang et al. [START_REF] Zhang | Discrete shape modeling for skin model representation[END_REF], a discrete skin model shape containing 𝑁 vertices is described as:

𝑆𝐾 = [𝑣 1,𝑥 , 𝑣 2,𝑥 … 𝑣 𝑁,𝑥 , 𝑣 1,𝑦 , 𝑣 2,𝑦 … 𝑣 𝑁,𝑦, 𝑣 1,𝑧 , 𝑣 2,𝑧 … 𝑣 𝑁,𝑧 ] 𝑇 (22)
where (𝑣 i,𝑥 , 𝑣 i,𝑦 , 𝑣 i,𝑧 ) are the coordinates of the 𝑖-th vertex. Given 𝑚 skin model shapes of one part (𝑆𝐾 1 , 𝑆𝐾 2 … 𝑆𝐾 𝑚 ), which are generated by simulation or measurement as the training set, the mean model 𝑆𝐾 ̅̅̅̅ is:

𝑆𝐾 ̅̅̅̅ = 1 𝑚 ∑ 𝑆𝐾 𝑖 𝑚 𝑖=1 (23)
Then the covariance matrix could be calculated:

𝑙 𝜌 = 5 𝑙 𝜌 = 10 𝑙 𝜌 = 20 𝐶𝑜𝑣 = 1 𝑚 ∑ (𝑆𝐾 𝑖 -𝑆𝐾 ̅̅̅̅ )(𝑆𝐾 𝑖 -𝑆𝐾 ̅̅̅̅ ) 𝑇 𝑚 𝑖=1 (24)
Eigenvalues and corresponding eigenvectors which reflect the principal components of the form error could be evaluated from the covariance matrix. Assuming the largest 𝑘 eigenvalues and eigenvectors are selected, then new skin model shapes could be generated from these eigenvectors:

𝑆𝐾 = 𝑆𝐾 ̅̅̅̅ + 𝝃𝝐 (25)
where 𝝃 is the matrix of first 𝑘 eigenvectors and 𝝐 is the random coefficient vector. The precision of reconstruction depends on how many eigenvectors (the size of 𝑘) are employed. The coefficient vector 𝝐 is floating around its mean value randomly and models with similar form defects are generated.

The covariance matrix 𝐶𝑜𝑣 could be considered as a Laplacian matrix as in Equation ( 26),

where 𝑑 is the relative deviation from the mean value. However, one dramatic difference between PCA and the other spectral method is that the geometry information of the nominal model is not used. The 𝐶𝑜𝑣 measures the relation between vertices from the point of view of deviation value (referred to mean value), not geometrical or topological structures of the nominal mesh model.

{ 𝐴(𝑖, 𝑗) = ∑ -𝑑 𝑖 𝑑 𝑗 𝑚 𝑘=1
, 𝑖 ≠ 𝑗 𝐷(𝑖, 𝑖) = ∑ 𝑑 𝑖 2 𝑚 𝑘=1 [START_REF] Schleich | Skin Model Shapes: A new paradigm shift for geometric variations modelling in mechanical engineering[END_REF] Unlike former methods which could be conducted without measurement data, the SSA could not generate new models without measurement or simulation data. When training data is available, however, the form error could be predicted and more models could be generated. using various methods, and several typical decomposition modes generated by the PCA method are shown in Figure 2-20. We recommend that readers look at [START_REF] Zhang | Discrete shape modeling for skin model representation[END_REF][START_REF] Ballu | Comparison of mode decomposition methods tested on simulated surfaces[END_REF] for more details.

Analysis and Comparison of Methods

In this section, we analyze the characteristics of different methods, then compare them. Since each method has its application condition, it is not possible to distinguish the best one from only one criterion. We will first consider the simulation methods from different criteria independently.

Then a global review is provided for each method, taking its engineering application into consideration. These will then help researchers find the best solution for them.

Criteria

The selected criteria are: multi-scale, surface complexity, measurement data integration, parametric control, calculation efficiency.

Multi-Scale

A machined surface contains irregular deviations of different orders. These deviations are classified as form, flaw, waviness and roughness at different scales. They may affect the mechanical properties such as real area of contact, friction, wear and lubrication. Depending on the problem to analyze, different levels of deviation are concerned. Therefore the ability to express multi-scale manufacturing defects shows the potential use of the method.

For the 1D Gaussian random noise method, because the deviation between connected vertices is independent, it is limited to simulate defects at a certain level. This deviation level is also decided by the point cloud density. Therefore its simulation results are more local and scattered. With second order shape morphing, we encounter the opposite problem. It is easy to simulate global form defects with simple functions. However, for some high frequency local defects, the shape function will be complex, and there is no such research as yet. Thus both 1D Gaussian random noise and second order shape morphing are not flexible enough to simulate multi-scale manufacturing defects.

The rest of the methods are sufficient to simulate multi-scale defects. For random shape morphing, this could be done by changing the control point number and the size of influence hull.

For mode based methods, the deviation frequency is easily chosen.

Surface Complexity

In tolerance analysis of mechanisms, the planar and cylindrical mating pairs are the most common ones, while for tolerance analysis of the assembly of sheet metal parts, more complex shapes are encountered. It will be better if the simulation method could be used for any general cases without having to consider the shape of the part. Methods using shape functions encounter difficulties for this criterion. Second order shape morphing established the function globally. For some complex surfaces, the construction of the shape function for morphing will be a difficult job. Trigonometric function modes are usually constructed on a Cartesian coordinate system, then mapped to other coordinate systems to fit the surface shape. In the work of Song et al. [START_REF] Song | Shape deviation modeling for fused deposition modeling processes[END_REF], mapping between a Cartesian coordinate system and a polar coordinate system was used.

As an exception, random shape morphing could handle complex surfaces without additional operations. In fact, it can be applied to any kind of mesh without considering their genus. Spectral mode based methods could also be applied to general surfaces.

Measurement Data Integration

In the early design stage, in order to conduct tolerance analysis, the manufacturing defects were simulated randomly. For many situations, however, former process data and measurement data may be available. Full utilization of these data will generate a more realistic simulation result and benefit the optimization of design and manufacturing.

The data integration property is linked with multi-scale defect simulation ability. The 1D Gaussian random noise method uses the mean value and standard deviation of the measurement data. This is compact, but its simulation ability is also limited. The second order shape morphing method is used more to simulate global shape deformations, which ignore local defects like roughness. Therefore the measurement data is filtered.

Methods that have the ability to simulate multi-scale defects are also able to fully utilize the measurement data. Mode based methods have the advantage of decomposing measurement data by mode basis. Suitable construction of the mode basis could generate a better approximation of the measurement data. After measurement data decomposition, the weight coefficients for each mode are calculated. The reconstruction of a similar model is based on the combination of modes using weight coefficients.

The random shape morphing method could decompose measurement data, but the choice of control parameters may be difficult. For random shape morphing, too many parameters could be customized, like the shape of basic function, the size of influence hull, and the number and position of control points. This will be a problem if the user does not have any prior knowledge about the measurement data. The mode based method is easier as orthonormal bases with similar patterns are used.

Parametric Control

The shape of deviation is the main concern in manufacturing defect simulation. A flexible and parametric control of the deviation shape will be helpful in further applications.

The random shape morphing method may not be an easy solution for measurement data decomposition, but it is the most flexible method. The shape of simulated defects is controlled by the shape of basic function, which can be modified explicitly and easily. The method could simulate almost all the shape defects needed. This is the same for second order shape morphing.

The shape of defects could be precisely controlled by second order functions.

The 1D Gaussian random noise method uses only mean value and standard deviation, thus it only controls the statistical characteristics rather than the detailed defect shape. In mode based simulation methods, the user can choose the modes and their coefficients, which control the defect shape indirectly.

Calculation Complexity

Calculation complexity is important when users need to simulate a batch of models, or the density of mesh/point cloud is high. Simulation of skin model shape is only the linear combination of deviation data, as explained in Equation [START_REF] Charpentier | A scientific point of view of a simple industrial tolerancing process[END_REF]. Therefore the complexity of deviation data calculation plays a critical role. Spectral mode based methods calculate eigenvectors before simulation. With the increase in mesh density, the time used to calculate the same frequency eigenvector also increases dramatically. Moreover, the available mathematical software usually calculates from the smallest or largest eigenvalue. This means that if we need the eigenvector with medium frequency, we always need to calculate half of all eigenvectors.

Considering the matrix we used for the calculation of eigenvectors, for example the matrix of natural vibration analysis, its entries correspond to the topology of the mesh. The matrix will be sparse when it is large in size. For the PCA method, relations between all vertices are established and the matrix is dense. The matrix from random field may vary between sparse and dense depending on the correlation length 𝑙 𝜌 and autocorrelation function 𝜌 𝑖,𝑗 . Depending on the method we choose and the modes we want, suitable algorithms may be chosen to accelerate the calculation of eigenvectors.

Based on the discussions above, we provide a comparison of these methods, as shown in

Advantages, Drawbacks and Application of Each Method

Among the methods mentioned above, the 1D Gaussian method is quite different from the rest. From morphing to mode based methods, relations are established between neighboring vertices to generate a smooth shape deviation. With the 1D Gaussian method, however, independent and chaotic deviations are simulated.

For a real manufactured surface, both the chaotic deviation (random noise) and smooth deviation could be found, but on a different scale. On an actual surface, random noise corresponds to roughness or waviness based on experience. The smooth deviation could be considered as the result of a filtration operation, which filters random noise from the actual manufactured surface [START_REF]Geometric Product Specification-General concepts-Part 1: Model for geometrical specification and verification n[END_REF][START_REF] Ballu | Global view of geometrical specifications[END_REF]. Generation of realistic surfaces needs both smooth and chaotic deviations. In a first approximation, the 1D Gaussian method is able to introduce roughness or waviness. However, depending on the process, roughness or waviness are not really random noises. For machining, they are irregular deviations based on the geometry of machine tool vibration. With the 1D Gaussian method it is also possible to introduce measurement uncertainty of the surface.

Both mode based methods and the 1D Gaussian method could be used to simulate random noise. For mode based methods, with the increasing frequency, the deviation of the mode tends to be more local and irregular. By taking high frequency modes, it is possible to simulate the random noise, but this needs to take more modes and chose the coefficients carefully to eliminate artifacts (as in Figure 2-17), which may require a large amount of computation. This is different with the 1D Gaussian method. Because the deviation between each vertex is independent, the calculation complexity has a linear relation with the number of vertices, which could be much faster than mode based methods when the mesh of the surface is fine.

On the other hand, to simulate smoother and continuous surfaces, the 1D Gaussian method is no longer suitable. Instead, mesh morphing or mode based methods could be used.

Using a second order shape morphing method, manufacturing deviations caused by the kinematic error during machining could easily be simulated, as in Figure 2-8 and Figure 23456789. These kinematic errors are sometimes called "systematic errors", due to their connections with particular manufacturing processes [START_REF] Zhang | Discrete shape modeling for skin model representation[END_REF]. It should be noted that even this systematic error is close to some kind of second order shape, it may also vary between each part around a mean value. To simulate these deviations, random parameters should be introduced to the functions of second order shape. Meanwhile, most of the simulation cases are dealing with plane and cylinder features [START_REF] Zhang | Discrete shape modeling for skin model representation[END_REF][START_REF] Kurokawa | Evaluation of shot peened surfaces using characterization technique of three-dimensional surface topography[END_REF], and the shape of deviation is simple. This makes it easy to establish the second order functions. For more general non-primitive surfaces with multi-scale complex form error, it will be hard to choose a suitable function.

Random mesh morphing with basic functions is a flexible method that could be used for any kind of mesh surfaces. The explicit form of basic function makes it easy to modify the simulation parameters, and could simulate both smooth and chaotic deviations. For example, to simulate chaotic deviations like the 1D Gaussian method, we could assume that the influence hull of the basic function is smaller than the edge length of the mesh, thus each vertex is independent. To simulate smoother deviations, we can just increase the size of influence hull, and choose suitable basic functions and parameters. In this respect, it is the most general and efficient simulation method. However, the link between the simulation parameters and engineering measurement is not yet established. A unified process is needed to decide the parameters and basic functions.

Further research could be done to link the simulation parameters with measurement results to bring the simulation result closer to engineering practice.

Trigonometric function modes provide a compact expression of deviations. With only a few variables, they express deviations from form to roughness [START_REF] Huang | Mode-based Decomposition of Part Form Error by Discrete-Cosine-Transform with Implementation to Assembly and Stamping System with Compliant Parts[END_REF][START_REF] Armillotta | Tolerance Analysis Considering form Errors in Planar Datum Features[END_REF]. The limitation of this method concerns the complexity of the surface. Geometrical mapping between the complex surface and the deviations of a plane should be established. This is easy with primitive geometries, like a cylinder. To extend this method to general and complex surfaces, the mapping process may be complicated or infeasible.

Spectral theory based methods established the relations between vertices (based on topology of the mesh or measurement data), and generate both smooth and chaotic deviations. The advantages of these kinds of method are that they can use the unified processes for multi-scale deviation simulation, and are not limited by the shape of the nominal model. In [START_REF] Schleich | Skin Model Shapes: A new paradigm shift for geometric variations modelling in mechanical engineering[END_REF][START_REF] Samper | Modeling of 2D and 3D Assemblies Taking Into Account Form Errors of Plane Surfaces[END_REF], methods to evaluate the coefficients of modes are established, and they are connected to statistical parameters like standard deviation. This is an easy and practical way to analyze measurement data without prior knowledge, and could be used to describe certain types of manufacturing defect [START_REF] Dantan | Geometrical variations management for additive manufactured product[END_REF].

One problem of spectral theory based methods is the calculation of modes. To simulate deviations at different scales, or decompose measurement data, the more modes we use, the more realistic the result could be. However, the calculation of eigenvectors could be quite time consuming.

Conclusion

Geometrical quality control starts from the product design stage, continues into manufacturing, inspection, until recycling. To guarantee geometrical quality and the related functionality of the product, CAT systems have been developed to simulate manufacturing deviations and virtual mockup. In most commercial CAT systems, deviations are generated by transforming nominal features, which do not contain more details of the manufactured surfaces.

This approximation could not satisfy the increasing demand for product precision, and is not sufficient to cover the whole product lifecycle. Some researchers have started to introduce virtual parts into simulation, with detailed manufacturing deviations, like form error. These are still laboratory prototypes. Two problems are encountered when introducing form error into the CAT systems: 1) how to simulate or regenerate the manufacturing deviations on the part, and 2) how to assemble parts with form errors under complex processes.

The aim of this chapter was to find suitable simulation solutions for the skin model shape, which is a discrete representation of the real manufactured part. We classified the current simulation methods into three categories, random noise methods, mesh morphing methods and mode based methods. This classification is not based entirely on their origins, but also on their mathematical and computational approaches. This is to give the reader an idea of what kind of result one method could generate. For example, even the natural vibration modes and random fields are quite different methods considering their origins, but in application, both of them combine eigenvectors.

Moreover, existing simulation methods already cover a lot of ranges and could satisfy most applications. It is not possible to find a method which is best for all aspects, therefore we defined criteria to compare them from different points of view. These criteria included multi-scale, surface complexity, measurement data integration, parametric control, and calculation complexity. This review and comparison could help researchers and engineers to have a better understanding of these methods before application, and to avoid reinventing similar methods.

With the simulated manufacturing defects, methods to generate complete and consistence skin model shapes will be explained in next chapter.

Skin Model Shape Generation

In this chapter, using the simulated form defects, skin model shapes are generated.

Inconsistency between non-ideal surfaces is highlighted, and solutions are proposed to solve this problem. Examples are given to compare the differences between results from local and global solutions.

Introduction

Given the limitations of nominal model-based methods, many methods have been developed

to represent detailed geometric deviations. In Chapter 2 we have provided a review of form error simulation methods, and the method classification is shown in Figure 2-2. Concerning these simulation methods, all of them may be applied to simple canonic surfaces (plane, cylinder, sphere, and torus). Some of them may be applied to a set of canonic surfaces or to the whole surface of a complex part: 3D random noise, morphing.

After the independent simulations, the generated surfaces must be combined to create the skin model shape. As the deviations are simulated independently for each feature, and as the corresponding features have connections at their boundaries, combining the independent deviation data is a problem. Geometric inconsistencies appear at the edges when the whole part surface is reconstructed from the individual ones.

The objective of this chapter is to present geometrical issues involved when combining the individual features. From these issues, different methods are proposed to obtain fully consistent skin model shapes for any kind of surface, particularly a method based on Finite Element Analysis (FEA). Finally, the method is applied to a trial part and to a representative part from the mechanical industry. This corresponds to the simulation process in Figure 2-1(d).

Geometric Issues

This section highlights the difficulties encountered in combining deviation data from each feature and adding them to the nominal model. The following sections are illustrated with 2D examples, as in Figure 3-1(b) (assume the vertices lie in the same section plane), but they are representative of the main problems in 3D space.

The vertex normal is estimated by Mean Weighted Equally algorithm [START_REF] Jin | A comparison of algorithms for vertex normal computation[END_REF]. Figure 3-2 shows the difference when we consider each feature independently or not. To simulate manufacturing deviation, which has independence between features, the situation in

Non-Connection

Unlike the features in the nominal model, which are connected at the edges, independent features with form deviations are no longer connected when we put them together. As can be seen from the figures, depending on the configuration type and deviation value, sometimes the two features are not connected while sometimes they intersect. The 3D case is much more complex. Figure 3-4 illustrates on a 3D view the fact that the surfaces may have zones of connection and non-connection on one and the same edge. To guarantee the topological coherence of skin model shape, this issue has to be solved.

Face Connection

It is possible to calculate the intersection between features and delete the intersect parts, or calculate the split and fill in more vertices. In [START_REF] Chen | A sharpness-dependent filter for recovering sharp features in repaired 3D mesh models[END_REF][START_REF] Ju | Fixing geometric errors on polygonal models: a survey[END_REF][START_REF] Yang | Shape-controllable geometry completion for point cloud models[END_REF], mesh repairing methods are developed. This is not easy to carry out and the topology of the mesh needs to be changed. Meanwhile, it leads to more problems, such as the deviation simulation of new added vertices, the choice of new edge or corner, etc. Thus a simpler method is preferred.

The simplest way to connect the faces is to add the simulated deviation to the vertices of the nominal model (discrete but not segmented into separate features) directly. The deviations are added along the vertex normals. This makes a vertex on an edge, which belongs to two faces, has two deviations along the normal direction of the two faces at the vertex, and the final vertex is the result of the addition of the two deviation vectors.

Obtuse and Acute Dihedral Angles

However, when the connected features are at an obtuse or acute dihedral angle, there are new issues. The precision of each feature is influenced during the combination process by vector addition. Figure 34567(a) shows the situation when two connected features are at an obtuse dihedral angle. By the principal of vector addition, the position of the edge vertex after combination can be calculated easily. However, this point does not correspond to the desired point. It is clear that by the addition of vectors, shape error for both single features is increased. For the acute dihedral angle in Figure 34567, the deviation of the edge vertex decreases in a similar way.

The desired point is the intersection of two offset lines parallel to the two nominal lines; the two values of the offsets are equal to the two deviations of the edge vertex according to the two lines. Thus, the intersection point keeps the deviation of the edge vertex along the normal direction the same as before the deviation combination. In the proposed method, the objective will be to obtain this intersection point.

Ratio Mesh Size/Deviation Magnitude

The visualization of the form error shape is also important, as it could help engineers to understand potential defects of parts, make optimal decisions and improve the manufacturing process [START_REF] Ballu | Virtual Metrology Laboratory for e-Learning[END_REF]. Besides using the contour map of deviation, an amplified form error shape can show the deviation directly. Thus, we amplify the form errors on the skin model shape. However, this amplification also causes problems when a combination is conducted.

This can be seen in Figure 3-8. The amplified deviation may become larger than the mesh size, as shown in Figure 3-8(a). When the combination method is used by direct deviation addition to a nominal model, mesh triangles connected to the edge vertex will be over-stretched or intersect with each other.

Deviation Combination Methods

Several methods have been studied to combine the geometric deviations in a single mesh without intersection or stretching, while guaranteeing the precision of the feature. Based on their principles, they are classified as a local method or a global method.

Local Method

In local methods, the deviations are first added to the nominal model, then the mesh regularization process can be conducted to eliminate any intersection or split on the mesh [START_REF] Schleich | Skin Model Shapes: A new paradigm shift for geometric variations modelling in mechanical engineering[END_REF].

The general process used for local methods is explained below:  Computation of deviations for each feature independently.

 Addition of the deviations to the nominal model. For vertices inside the faces, the deviations are added along the normal direction directly. For vertices on the edges or corners, the positions of the vertices are calculated as the intersection of local tangent planes, as explained before.

 Regularization of the mesh to handle mesh problems with the skin model shape. Two methods are used, one based on Laplacian mesh regularization and the other on spring analogy.

Next, the principles of the two regularization methods are explained.

Laplacian Mesh Regularization

The Laplacian mesh regularization method [START_REF] Ohtake | Mesh regularization and adaptive smoothing[END_REF], which is simple and fast, has been used to regularize the mesh in many different applications [START_REF] Riehl | A staggered approach to shape and topology optimization using the traction method and an evolutionary-type advancing front algorithm[END_REF][START_REF] Yamauchi | Feature sensitive mesh segmentation with mean shift[END_REF].

Laplacian regularization is based on the umbrella operator [START_REF] Kobbelt | Interactive multi-resolution modeling on arbitrary meshes[END_REF] to adjust the position of the vertices of the mesh repeatedly. As explained in Figure 3456789, the umbrella operator 𝑼(𝑷) of a vertex 𝑷 is a vector, which is defined by the following expression:

𝑼(𝑷) = 𝟏 ∑ 𝒘 𝒊 𝒊 ∑ 𝒘 𝒊 𝑸 𝒊 𝒊 -𝑷 (27)
where 𝑷 is the vertex to be adjusted and 𝑸 𝒊 are the neighbors of 𝑷. The 𝑤 𝑖 are positive weights which can be adjusted, but one simple way is to set 𝑤 𝑖 = 1. The position of vertex 𝑷 is adjusted by adding a scaled umbrella operator 𝝀 • 𝑼(𝑷):

𝑷 𝒏𝒆𝒘 = 𝑷 𝒐𝒍𝒅 + 𝝀 • 𝑼(𝑷 𝒐𝒍𝒅 ) (28)
𝑷 𝒐𝒍𝒅 and 𝑷 𝒏𝒆𝒘 are the coordinates of vertex 𝑷 before and after adjustment, 𝜆 > 0 is a small positive number. Every vertex coordinate is adjusted by iteration with Equation ( 28), and parameter 𝜆 decides the adjustment speed in each iteration. The operator 𝑼(𝑷) can be decomposed according to the normal vector 𝒏 and tangent vector 𝒕 at vertex 𝑷, then:

𝑷 𝒏𝒆𝒘 = 𝑷 𝒐𝒍𝒅 + 𝝀 𝟏 • 𝑼 𝒏 (𝑷 𝒐𝒍𝒅 ) + 𝝀 𝟐 • 𝑼 𝒕 (𝑷 𝒐𝒍𝒅 ) (29)
where 𝑼 𝒏 (𝑷) and 𝑼 𝒕 (𝑷) are the normal and tangential components of 𝑼(𝑷). By optimizing the calculation method of 𝑼(𝑷) and regularization speed 𝜆, the method is improved by Ohtake et al. [START_REF] Ohtake | Mesh regularization and adaptive smoothing[END_REF]. The Laplacian mesh regularization method is efficient and easy to apply. As shown in Equation [START_REF] Geis | Use of Vectorial Tolerances for Direct Representation and Analysis in CAD-systems[END_REF], it is possible to adjust the position of the vertices along a specific direction. Thus, to conserve the geometric precision of features, we adjust the vertices inside a surface only by the tangential component of the umbrella operator. This is illustrated in Figure 3-10(a), where vectors 𝒕 and 𝒏 correspond to the tangential and normal vectors respectively. After mesh regularization, the manufacturing defects, which are expressed by the deviations along vertex normal directions, are not influenced.

For vertices on the edge, the umbrella operator is projected to a direction which is perpendicular to the two normals of the local tangent planes of vertex 𝑷. As explained in Figure 3-10(b), 𝒏 𝟏 and 𝒏 𝟐 are the normals of local tangent planes when vertex 𝑷 is on the edge of two connected faces 𝑺 𝟏 and 𝑺 𝟐 . Direction 𝒏 𝟑 is perpendicular to both 𝒏 𝟏 and 𝒏 𝟐 , thus adjustment along 𝒏 𝟑 will not influence deviation along these two normal directions.

For vertices on the corner, they already have deviations along three non-parallel directions, thus their position will be fixed during the mesh regularization process.

Q 1 Q 2 Q 3 P Q i U(P)
Figure 3-10 Projection of umbrella operator during mesh regularization process.

Spring Analogy

The spring analogy is commonly used in Computational Fluid Dynamics (CFD) for moving boundary problems and mesh regularization [START_REF] Ait-Ali-Yahia | Anisotropic mesh adaptation: towards user-independent, mesh-independent and solver-independent CFD. Part II. Structured grids[END_REF][START_REF] Blom | Considerations on the spring analogy[END_REF][START_REF] Zeng | A semi-torsional spring analogy model for updating unstructured meshes in 3D moving domains[END_REF]. Due to its simple and basic assumption, modification of the mesh can be done efficiently.

Spring analogy applies Hook's law to adjust the mesh; the principle is explained in Figure 3-11.

Based on the definition of equilibrium length, the method can be classified into vertex springs and segment springs [START_REF] Blom | Considerations on the spring analogy[END_REF], and their equations are generalized. Let point 𝑷 be the vertex to be adjusted, and 𝑸 𝒊 its neighbor vertices, the force applied to vertex 𝑷 is:

𝑭 𝑷 = ∑ 𝒌 𝒊 (𝑸 𝒊 -𝑷 -𝒅 𝒊 ) 𝒊 ( 30
)
𝑘 𝑖 is the stiffness of the spring, and 𝒅 𝒊 is the equilibrium vector. For vertex springs, the equilibrium vector is zero, while for segment springs, 𝒅 𝒊 is defined as:

𝒅 𝒊 = 𝑸 𝒊,𝒐𝒍𝒅 -𝑷 𝒐𝒍𝒅 (31)
where the subscript "old" indicates the initial position of the vertices. The displacement vector of vertex 𝑷 is:

𝜹 𝑷 = 𝑭 𝑷 ∑ 𝒌 𝒊 𝒊 (32)
By solving the equations iteratively, the position of the vertices can be calculated:

𝑷 𝒏𝒆𝒘 = 𝑷 𝒐𝒍𝒅 + 𝜹 𝑷 (33)
For vertex springs, if we set the stiffness 𝑘 𝑖 = 1, it is the same as the simplest umbrellaoperator 𝑼(𝑷), which has been explained above. Both the vertex spring and segment spring methods can be improved by modifying the stiffness 𝑘 𝑖 .

As with the 2D spring analogy in its application to our problem, 3D springs are considered to calculate vertex displacement. Next, to conserve the geometric deviation of the features, displacement 𝜹 𝑷 is projected in a specific direction which is the same as in the Laplacian mesh regularization method.

2D Example by Local Method

The example in Figure 3

Global Method

In the local methods proposed above, the mesh is adjusted with certain direction constraints.

Based on this principle, we propose another method to add direction constraints and solve the regularization problem globally.

Q 1 Q 2 Q 3 Q 4 P Q 5 Q 6 (a) (b)
In our combination and adjustment problem, where every vertex has at least one adjustment constraint, the method for adding constraints and solving equations efficiently is important. Finite Element Analysis (FEA) which has a standard solving process and can deal with complicated boundary problems, is a promising way to solve the issue.

FEA is an efficient tool in structural design and analysis. It decomposes complex structures into small basic elements and provides a standard solving process. This benefits its implementation in computer programs, and complex structural problems are solved as simpler ones with the same process. In calculating CFD problems, where mesh is moving during simulation, FEA is also used for mesh regularization [START_REF] Amirante | A moving mesh algorithm for aero-thermo-mechanical modelling in turbomachinery[END_REF][START_REF] Hsu | Mesh deformation based on fully stressed design: the method and 2-D examples[END_REF].

With FEA, the 3D bar element is used to calculate the adjustment of vertex, which is similar to the spring in spring analogy method. The penalty function approach is used to add adjustment direction constraints. The regularization problem defined by FEA is a linear problem, and the calculation only needs to be conducted once. The process of combination with FEA is shown below:

• Simulation of deviations on each feature independently. The deviation is not added to the vertex of the feature immediately, but stored as a displacement boundary constraint along the normal direction of the vertex.

• Generation of matrix for FEA based on the existing mesh model. Every edge of the mesh triangles is taken as a 3D bar element.

• Definition of the boundary displacement constraints corresponding to the deviations.

• Solving the FEA problem using the penalty function approach [START_REF] Zienkiewicz | The finite element method[END_REF] to generate the regularized skin model shape.

In the following, the basic concept of FEA and the penalty function approach are introduced.

Finite Element Analysis

Typically, for a static FEA, the stiffness equation is applied to the considered elements:

𝑲 𝒆 ⋅ 𝒒 𝒆 = 𝑷 𝒆 (34)
where 𝑲 𝑒 is the stiffness matrix, 𝑷 𝑒 is the nodal force vector, and 𝒒 𝑒 is the nodal displacement vector which needs to be solved for each element. They are defined in a local coordinate system associated with the corresponding element. To conduct the analysis and to solve the displacements of all the nodes in the structure, they should be transformed into a global coordinate system and assembled together. With the transformation matrix 𝑻 𝑒 , we obtain:

{ 𝑲 𝒆 = 𝑻 𝒆 𝑻 𝑲 𝒆 𝑻 𝒆 𝑷 𝒆 = 𝑻 𝒆 𝑻 𝑷 𝒆 𝒒 𝒆 = 𝑻 𝒆 𝒒 𝒆 (35)
𝑲 𝑒 , 𝒒 𝑒 and 𝑷 𝑒 are the stiffness matrix, displacement vector, force vector for elements in a global coordinate system respectively. Once the displacement vector is calculated, other parameters (like strain, stress) can be solved based on their relation. For more detailed information about FEA, such as how to establish the stiffness equation and solve it, we refer readers to [START_REF] Zeng | Fundamentals of Finite Element Analysis[END_REF][START_REF] Zienkiewicz | The finite element method[END_REF].

Penalty Function Approach

To solve the stiffness equation with specified boundary conditions, there are several methods, such as direct solving, the Lagrange multiplier method and the penalty function approach. As the penalty function approach keeps the structure of the original equations (matrix size, ranking and symmetry), it reduces the complexity of large-scale computing dramatically, and has been widely used in engineering applications. We apply this method to our mesh adjustment. The principle of the penalty function approach is shown below.

For a finite element model, its potential energy function Π is:

𝚷 = 𝟏 𝟐 𝒒 𝑻 𝑲𝒒 -𝑭 𝑻 𝒒 ( 36
)
Where 𝒒 is the displacement vector of the system, 𝑲 is the stiffness matrix, and 𝑭 is the external force vector added to the nodes. To introduce the displacement boundary conditions, the displacement relation between coordinates can be written as:

𝑪𝒒 = 𝒅 (37)
𝑪 is the matrix containing the displacement relations between the coordinates of the vertices, and 𝒅 is the vector of displacement boundary conditions (corresponding to the simulated form deviations along normal directions in our case). A very large number 𝑆 is introduced, modifying the potential energy function as:

𝚷 * = 𝟏 𝟐 𝒒 𝑻 𝑲𝒒 -𝑭 𝑻 𝒒 + 𝟏 𝟐 𝑆(𝑪𝒒 -𝒅) 𝟐 (38)
Since S is a very large number, Π * can take its minimum value only when (𝑪𝒒 -𝒅) ≈ 0, so that the displacement boundary condition is satisfied. For the minimum value, ∂Π * ∕ ∂𝒒 = 0, the resulting equation can be written as:

𝑲 * 𝒒 = 𝑭 * (39)
where 𝑲 * = 𝑲 + 𝑆𝑪 𝑇 𝑪, 𝑭 * = 𝑭 + 𝑆𝑪 𝑇 𝒅.

For identical finite element model and boundary conditions, the result of the calculation is influenced by the parameter 𝑆 in the penalty function. If 𝑆 is large enough, the boundary conditions play a critical role and the adjustment will be only along the tangential direction of the vertices. Thus, depending on the size and precision of the model, a relatively large 𝑆 is taken to guarantee the effectiveness of adjustment.

2D Example by FEA-Based Method

Figure 3-13 illustrates the FEA-based mesh regularization method in a simplified 2D situation.

The displacement boundary constraint is added to every vertex. The value of the displacement boundary constraint is manufacturing deviation, which has been simulated for every vertex, and the directions are vertex normals.

For a vertex with just one displacement boundary constraint, it can be adjusted along its local tangent direction. For a vertex on the edge, which is shown in the dashed circle region, this has two non-parallel displacement boundary constraints, and its position is fully constrained in 2D.

When dealing with the generation of skin model shapes, a similar method is extended to a 3D situation.

Simulation and Comparison

Model Used for Comparison

As the vertices and triangles on the boundary cause problems in combining features with form error, here we use the model with four typical edges in mechanical parts to evaluate the simulation result. This model is shown in Figure 3-14. The model is initially discretized and segmented into 7 features, with different colours given to distinguish the faces. The four edges, e1, e2, e3 and e4, represent the edges between acute angle, obtuse angle, right angle and edge with two tangent faces, respectively.

Simulation Results and Comparison

To assess the effectiveness of different methods, we compare the simulation result by: (1) the direct deviation addition method, (2) Laplacian mesh regularization and (3) FEA-based mesh regularization. The form error for each feature is first simulated by the random field method, which has been introduced in earlier sections. To visualize the shape error of the model, the deviation is amplified. It should be noted that all three methods use the same result from the random field simulation, but differ in how the deviations are combined and adjusted.

As in the comparison method, deviations are assigned directly to the discrete skin model shape as shown in Figure 3 the problems on the edges are eliminated, and the simulated manufacturing deviations are conserved. For this method, criteria are required (e.g. not contain mesh intersections, threshold value of smallest adjustment, etc.) to decide whether or not to continue the iterations.

The result for the FEA-based regularization method is given in Figure 3-15(d). As can be seen from the zoomed in areas, a more regular and smoother edge is generated compared to Figure 3-15(b, c), and only one linear problem needs to be solved. Considering its complexity and effectiveness, the FEA-based mesh regularization could be the best choice for our problem. Based on the discussion above, by applying the adjustment with the FEA method, a reliable combination process could be conducted and the skin model shape is generated. However, one thing that should be considered is the connection of the tangent area, as shown in e4 of Figure 3-14. When we separate the tangent connected features and assign manufacturing defects to them independently, dramatic changes may happen at the edge area, as in e4 in Figure 3-16(a). This is possible if the two tangent-connected features are machined in separate processes. If we machine them at the same time, a more continuous simulation result is preferred and we should regard the two connected features as one feature before simulation with the random field method.

Application with Mechanical Part

The model we use for method comparison is a simple one, but it does contain several different types of edge situation. To verify the effectiveness of the proposed simulation procedure and deviation combination method, a representative mechanical part with complex shapes is used to generate its skin model shape.

The model we used here is a cutter body for milling. The nominal shape of the model can be seen from

Conclusion

The objective of this chapter is to present different approaches to generate consistent skin model shapes, so they could be used in tolerancing, manufacturing and metrology. Geometrical problems during skin model shape simulation (self-intersection, increase of the deviation…) are given and analysed. Two solutions are proposed to solve these problems: the local method and the global method.

Local methods, like Laplacian mesh smoothing or spring analogy-based mesh adjustment, are used to regularize meshes at the edges. Due to local properties, iterative calculation and the need for stop criteria, they are not the best choice for us.

A global method based on FEA is then introduced to generate skin model shapes while avoiding the drawbacks of local methods. In FEA based method, the manufacturing deviations of each feature are considered as displacement boundary conditions. With the penalty function approach, the combination of manufacturing deviations and the adjustment of mesh (solve geometrical problems) could be done within only one global calculation without iterations.

Comparisons between local and global methods are conducted, and the FEA-based global method generate better result than local methods.

Based on the simulation procedure and methods we have proposed, consistent skin model shapes could be generated. However, the generation of skin model shapes is only the first step, and more research efforts are needed for further promising applications. In Chapter 4, assembly simulation is conducted using the generated skin model shapes. In Chapter 5, virtual metrology method is developed. The method is applied on single skin model shape and skin model shape assembly.

Skin Model Shape Assembly Simulation

In this chapter, the method for skin model shape assembly simulation is studied. By comparing assembly results using different types of model, the importance of assembly boundary conditions is emphasized. Based on the analysis, a Linear Complementarity Condition-based method is proposed. The effectiveness of the proposed method is validated by assembly examples.

Introduction

Due to the inaccuracies of manufacturing, deviations exist not only between nominal part and real part, but also between real parts from a single batch. In complex mechanical assemblies, the accumulation of geometrical deviations may cause problems in assemblability and product functioning. Therefore, tolerance analysis and deviation evaluation become part of the main tasks in quality control.

Nowadays, the wide use of Information Technology (IT) facilitates tolerance management, Computer Aided Tolerancing (CAT) tools are developed to help engineers with decisions-making.

One of these tools, the Small Displacement Torsor (SDT), is used to evaluate and represent manufacturing defects [START_REF] Bourdet | Controlling a complex surface with a 3 axis measuring machine[END_REF]. The Offset Zone [START_REF] Requicha | Toward a Theory of Geometric Tolerancing[END_REF] based method is introduced to represent tolerance zone geometry. Technologically and Topologically Related Surfaces (TTRS) model, Vectorial tolerancing model [START_REF] Gao | Generalized 3-D tolerance analysis of mechanical assemblies with small kinematic adjustments[END_REF][START_REF] Wirtz | Vectorial tolerancing: a basic element for quality control[END_REF], Jacobian-Torsor model [START_REF] Desrochers | Application of a Unified Jacobian Torsor Model for Tolerance Analysis[END_REF] and Gap Space model [START_REF] Zou | A gap-based approach to capture fitting conditions for mechanical assembly[END_REF] are proposed for tolerance analysis. Proportioned Assembly Clearance Volume (PACV) [START_REF] Teissandier | A computer aided tolerancing model: proportioned assembly clearance volume[END_REF] and Tolerance Map (T-Map ® ) [START_REF] Davidson | A New Mathematical Model for Geometric Tolerances as Applied to Round Faces[END_REF] are proposed to calculate the ranges of assembly deviation.

The above-mentioned tolerance analysis approaches use ideal features, which contain no form defects. To emphasize the influence of form defects and cover different precision levels of manufacturing defects, the skin model shape is used to conduct tolerance analysis [START_REF] Schleich | Approaches for the assembly simulation of skin model shapes[END_REF].

The use of skin model shapes for virtual assembly and tolerance analysis is intuitional, but it encounters several difficulties. For mating two ideal features (two planes for example), there is a single assembly result. For non-ideal models, however, like skin model shapes, different results could be generated. This is because the simulation methods are the simplification of real physical processes, and different simplification methods may lead to different results. In the work of Charpentier et al. [START_REF] Charpentier | A scientific point of view of a simple industrial tolerancing process[END_REF] the assembly of skin model shapes is analyzed theoretically, Lê et al. [START_REF] Lê | Experimental and theoretical investigations of mechanical joints with form defects[END_REF] and Ballu et al. [START_REF] Ballu | Experimental evaluation of convex difference surface for planar joint study[END_REF] compared the results between ideal feature-based simulation and real part assembly experiment. This research highlighted the differences between simulations and actual engineering situations. To simulate assembly and predict its precision, suitable models and methods should be chosen carefully to meet the accuracy requirements.

Assembly simulation can be classified into: (1) rigid body assembly, and (2) compliant body assembly. For rigid body assembly, Pierce et al. [START_REF] Pierce | A Method for Integrating Form Errors Into Geometric Tolerance Analysis[END_REF] proposed a method based on minimizing the distances between mating features. Schleich et al. [START_REF] Schleich | Approaches for the assembly simulation of skin model shapes[END_REF] compared several optimization methods that could be used for the assembly simulation. Armillotta [START_REF] Armillotta | Tolerance Analysis Considering form Errors in Planar Datum Features[END_REF] applied the optimization method for a planar feature assembly. In these optimization-based methods, the assembly processes are simplified to an objective and constraints. How process parameters (like assembly sequence and force) are included in the simulation is critical for a reliable result. To consider the assembly force during simulation, Samper et al. [START_REF] Samper | Modeling of 2D and 3D Assemblies Taking Into Account Form Errors of Plane Surfaces[END_REF] introduced a method based on difference surface and calculated the static equilibrium position. This method is also extended to the assembly of spherical joints [START_REF] Ledoux | Integrating form defects of mechanical joints into the tolerance studies[END_REF], cylindrical joints [START_REF] Homri | Tolerance analysis --Form defects modeling and simulation by modal decomposition and optimization[END_REF], and considering local deformations under loads [START_REF] Grandjean | On the role of form defects in assemblies subject to local deformations and mechanical loads[END_REF].

Based on this method, the simulation result could be closer to a real assembly than methods not considering assembly force. However it also has drawbacks in that only one pair of mating surfaces is considered for each mating step. In the work of Corrado et al. [START_REF] Corrado | Manufacturing signature and operating conditions in a variational model for tolerance analysis of rigid assemblies[END_REF][START_REF] Corrado | 3D Tolerance Analysis with Manufacturing Signature and Operating Conditions[END_REF], optimization-based assembly was improved by integrating the checking of force balance. This method could take several mating surfaces into account, and the external forces applied to the parts are considered.

However, due to the randomly generated rotation adjustment, the consistency of simulation results is not guaranteed.

The assembly of compliant parts is also complex. Besides the manufacturing deviations, deformations of parts during assembly should also be considered [START_REF]Geometrical product specifications (GPS)[END_REF]. To simulate these complicated processes, FEA-based methods are widely used. In the work of Liu and Hu [START_REF] Liu | Variation Simulation for Deformable Sheet Metal Assemblies Using Finite Element Methods[END_REF], and the work of Chang and Gossard [START_REF] Chang | Modeling the assembly of compliant, non-ideal parts[END_REF], the assembly process of sheet metal parts is broken down into four steps: place, clamp, fasten and release (PCFR). Based on this basic process cycle, extensive research has been done to improve simulation precision or efficiency [START_REF] Franciosa | Simulation of variational compliant assemblies with shape errors based on morphing mesh approach[END_REF][START_REF] Wagersten | Using Morphing Techniques in Early Variation Analysis[END_REF][START_REF] Huang | Mode-based Decomposition of Part Form Error by Discrete-Cosine-Transform with Implementation to Assembly and Stamping System with Compliant Parts[END_REF][START_REF] Das | Physics-driven Shape Variation Modelling at Early Design Stage[END_REF][START_REF] Huang | Simulation and integration of geometric and rigid body kinematics errors for assembly variation analysis[END_REF][START_REF] Hu | Simulation and analysis of assembly processes considering compliant, non-ideal parts and tooling variations[END_REF][START_REF] Falgarone | Variation Simulation During Assembly of Nonrigid Components. Realistic Assembly Simulation with ANATOLEFLEX Software[END_REF]. Moreover, to reach a high level of simulation precision, research has started to introduce the effect of local form deformations into a rigid body-based tolerance analysis method. As mentioned before, Grandjean et al. [START_REF] Grandjean | On the role of form defects in assemblies subject to local deformations and mechanical loads[END_REF] considered the local deformation of mating surfaces under loads. [START_REF] Pierre | Variational tolerancing analysis taking thermomechanical strains into account: Application to a high pressure turbine[END_REF]. Liu et al. [START_REF] Liu | Assembly Error Calculation with Consideration of Part Deformation[END_REF] considered deviations caused by temperature and load, and integrated different deviations by the Jacobian-torsor model. Gouyou et al. [START_REF] Gouyou | Tolerance analysis of overconstrained and flexible assemblies by polytopes and finite element computations: application to a flange[END_REF] combined polytopes with FEA for over-constrained assembly system analysis. Lustig et al. [START_REF] Lustig | Enhancement in coupling tolerance analysis and elastic deformations on the example of a Serial linear support system[END_REF][START_REF] Lustig | Linking Elastic Deformations and Fabrication Deviations in Computer-Aided-Analysis[END_REF] enhanced the function of commercial software by linking elastic deformations to the model. In this chapter, the rigid body assembly simulation is conducted using skin model shapes. The analysis of assembly boundary conditions is emphasized. A link between optimization objective and the assembly boundary conditions, which are the foundation of the optimization-based method, is established. A general assembly simulation method based on Linear Complementarity Condition (LCC) is proposed.

Pierre et al. introduced thermomechanical strains into tolerance analysis

The rest of the chapter is arranged as follows: Section 4.2 analyzes the influences of assembly boundary conditions. Based on this analysis, the LCC-based simulation is introduced in Section 4.3.

Lastly, in Section 4. [START_REF] Bourdet | Controlling a complex surface with a 3 axis measuring machine[END_REF], several examples are used to validate the proposed methods.

Assembly Considering Form Defects

The consideration of form defects during assembly imposes many more difficulties than when no form defects are considered. Thus, to model the influences of form defects on product assembly, there must first be a suitable investigation into the assembly configurations. In this section, the simulation models, assembly conditions and final assembly configurations are studied.

The study highlights the important elements and provides guidelines for developing a simulation method.

In the study, three different kinds of model will be used. Figure 4-1 illustrates the differences between these models by 2D "square" examples (whereas they are "cubes" in 3D). The most familiar one is the nominal model. In tolerance analysis, the relations between part defects and assembly precision cannot easily be modeled by functions. Thus, to provide a prediction of the assembly result, one solution is to create a virtual assembly of the parts with manufacturing defects. The manufacturing defects of the ideal surface-based deviation model already give rise to difficulties during assembly simulation, and the skin model shape is far more complex.

Experimental work by Ballu et al. [START_REF] Lê | Experimental and theoretical investigations of mechanical joints with form defects[END_REF][START_REF] Ballu | Experimental evaluation of convex difference surface for planar joint study[END_REF] points out the importance of assembly boundary conditions. For assembly simulation, two aspects are considered to solve the problem:

 Mating pairs and their forms, which correspond to the displacement boundary condition.

 Assembly force and torque, which correspond to the load boundary condition.

The following subsections are separated on the basis of the mating conditions. For each mating condition, the influences of different load conditions are analysed.

Mating of Two Planar Surfaces

To start from basic cases, we look at the mating of two parts with planar contact surfaces. The ideal surface-based deviation model contains translation and rotation deviations, while the mating planes are still ideal plane features. This means the plane could be expressed by a center point and a normal vector. Using normal vectors, the orientation of two parts after assembly could be determined. As in nominal models, the two planes can still mate without gaps, as can be seen from Figure 4-2(b).

Here, friction forces are ignored because:

 During assembly, the parts are assumed to have displacements to adjust their relative positions, which overcome the static frictions.

 Additional displacement constraints are added to constrain the DOFs.

Thus in the configuration in Figure 4-2(b), without friction, the top part would slide to the left because the planar surface is not normal to the force. As there is a degree of freedom, the

(a) (b) (c) (d)
tangential displacement is considered as constrained. This constraint replaces the friction force.

The concept of support, used in FEA, will be used in the following part of the paper to put these constraints in place.

The problem becomes difficult when considering skin model shapes in Figure 4-2(c). Due to the consideration of form defects, the two mating surfaces are no longer ideal plane features. The mating of two surfaces is also reduced to certain contact regions only. The assembly result depends on the local contacts between two surfaces, and the local contacts depend on:

 The form defects on mating surfaces.

 The position and direction of assembly forces.

As in the configuration in Figure 4-3(c), there could be tangential displacements due to local contacts which are not normal to the force. Thus, the tangential displacements are also constrained. indicates that when we are considering form defects in assembly, we should also consider assembly loads.

Assembly of Two Parts with Several Contacting Surfaces

The assembly of two parts is based on the mating of surfaces. Between two parts, more than one mating pair is usually used to fix the relative position in different directions. Figure 4-4(a)

shows the mating of two nominal parts. The two mating edges (mating planes in 3D) could be mated perfectly at the same time.

However, for ideal surface-based deviation models, due to the translation and rotation of nominal features, usually only one pair of edges or planes can be mated. The assembly configuration depends on the defects and loads on the models, some of which are shown in Figure c,d). The contacts between the surfaces of the parts may be surface contact, line contact or point contact. This makes it possible to find the configurations by searching and enumerating the boundary features, like points, edges and surfaces. Figure 4-4(c, d) shows the assembly of the same parts under different loads. The final positioning depends on both manufacturing defects and loads. Li and Roy [START_REF] Li | Relative positioning of toleranced polyhedral parts in an assembly[END_REF] examined the possible configurations using an ideal surface-based deviation model. In their work, they find extreme configurations as the limits. The searching method does not work for skin model shapes. Figure 4-5 shows three assembly configurations for the assembly of skin model shapes. As can be seen from the figure, it is possible to have contacts between surfaces at any point inside the surface, not limited to edges as for ideal surface-based deviation models. To express the form defects, discrete models are usually used [START_REF] Zhang | A discrete geometry framework for geometrical product specifications[END_REF][START_REF] Yan | Generation of consistent skin model shape based on FEA method[END_REF][START_REF] Zhang | Discrete shape modeling for skin model representation[END_REF]. The number of points, edges and facets is far greater than for ideal surface-based deviation models, making it unrealistic to search for them all.

Meanwhile, the contact positions are also influenced by assembly loads. As illustrated by Similar difficulties are also found with other kinds of mating pairs.

Assembly of Multiple Parts

When more assembly parts are considered, the problem is even more difficult. Not only is there an increase in mating surfaces and configurations, but the assembly loads can also be more complex. The nominal models are considered first. Without form defects, the three nominal models in Figure 4-7 are placed in their nominal position. In most research, assembly precision is studied considering only the geometrical deviation of ideal features, and not considering loads. This corresponds to the displacement boundary conditions. The examples above point out that the assembly forces (corresponding to load boundary conditions) also influence the final configurations of the assembly significantly. Thus it is not appropriate to ignore assembly loads. Based on the analysis, we proposed to introduce the load boundary conditions into the tolerance analysis.

In conclusion, to simulate an assembly that resembles the real situation, the following factors need to be considered:

 Manufacturing deviations introduced into the skin model shape.

 Relative displacements between skin model shapes with rigid body displacement.

 Displacement conditions:

o Non-interpenetration condition.

o Internal reaction forces.

o Displacement boundary condition.

 Load conditions:

o External forces and torques (load boundary conditions).

o Balance of internal and external loads.

 Multiple assembly part.

Methods to introduce these factors into the assembly simulation are explained in the next section.

LCC-Based Assembly Simulation

In the studies mentioned above, differences between models, loads and simulation results are explained. To make the simulation reflect a real physical assembly, researchers are trying to improve assembly simulation methods by introducing the influences of loads. Samper et al. [START_REF] Samper | Modeling of 2D and 3D Assemblies Taking Into Account Form Errors of Plane Surfaces[END_REF] introduced a method based on difference surfaces and calculated the static equilibrium position as an assembly result. Ledoux et al. [START_REF] Ledoux | Integrating form defects of mechanical joints into the tolerance studies[END_REF] applied this method to spherical joints. Corrado et al. [START_REF] Corrado | Manufacturing signature and operating conditions in a variational model for tolerance analysis of rigid assemblies[END_REF] improved an optimization-based simulation method by adding the balance checking process after assembly.

Apart from the methods already used for tolerance analysis, research in robotics and multibody dynamics is also proving interesting. In multibody dynamics, dynamic properties of parts are formulated based on methods like the Linear Complementarity Problem (LCP) [START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF][START_REF] Williams | RPI-MATLAB-Simulator: A Tool for Efficient Research and Practical Teaching in Multibody Dynamics[END_REF] or penalty functions [START_REF] Drumwright | A fast and stable penalty method for rigid body simulation[END_REF]. Contact dynamics based on the Discrete Element Method (DEM) are used to simulate the contact reaction between small particles [START_REF] Radjai | Contact dynamics as a nonsmooth discrete element method[END_REF][START_REF] Gilardi | Literature survey of contact dynamics modelling[END_REF]. Rigid body kinematics are used for computer graphics and game engines [START_REF] Bender | Interactive simulation of rigid body dynamics in computer graphics[END_REF]. Banerjee et al. [START_REF] Banerjee | Historical Origin and Recent Development on Normal Directional Impact Models for Rigid Body Contact Simulation: A Critical Review[END_REF] have provided a comprehensive review of rigid body contact simulation methods.

Inspired by the concept of Linear Complementarity Condition (LCC), this work provides a general simulation method considering assembly boundary conditions all together. However, these boundary conditions are usually nonlinear, which are difficult to model and solve. To handle this problem, the linearization-based method is proposed. In the first step, the factors proposed at the end of the last section are analyzed, and linearized mathematical constraints are established.

Based on these linearizations, the assembly is modeled using the LCC. The simulation process of the proposed method is explained at the end.

Integration of Manufacturing Defects

To introduce manufacturing defects into an assembly simulation, the corresponding skin model shapes are required. Discrete modeling of the part is widely used in CAE systems, where it has the advantage of modeling and visualization. Therefore in the assembly simulation, the discrete modeling method is adopted. This modeling uses vertices, edges and triangle patches. To explain our method, a simple 2D model is shown in Figure 4-12(a). As can be seen from the figure, two surfaces 𝑆 1 and 𝑆 2 are constructed using vertices and segments. The form of surface is controlled by the relative displacements of vertices on it.

Meanwhile, the sampling rate of vertices on the surface also influences what type of form can be simulated. Detailed simulation methods can be found in earlier articles [START_REF] Schleich | Skin Model Shapes: A new paradigm shift for geometric variations modelling in mechanical engineering[END_REF][START_REF] Yan | Generation of consistent skin model shape based on FEA method[END_REF][START_REF] Zhang | Discrete shape modeling for skin model representation[END_REF][START_REF] Homri | Tolerance analysis --Form defects modeling and simulation by modal decomposition and optimization[END_REF][START_REF] Samper | Form Defects Tolerancing by Natural Modes Analysis[END_REF].

Displacements between Skin Model Shapes

Rigid Body Displacement

The assembly of skin model shapes is a process where the models are repositioned to satisfy certain constraints. Therefore the expression of skin model shape displacement from its initial position to its final position is important. In our method, the skin model shapes are considered as rigid parts, and theories for rigid body displacement are used.

In Figure 4-12(a), the two surfaces are in their initial positions. The center points of the surfaces are 𝐶 1 and 𝐶 2 respectively. 𝐴 1 is a vertex on surface 𝑆 1 . After displacement, surface 𝑆 1 becomes 𝑆 1 ′ as shown in Figure 4-12(b). Its center point becomes 𝐶 1 ′ , and vertex 𝐴 1 becomes 𝐴 1 ′ .

Based on the rigid body displacement theory, the displacement of surface 𝑆 1 can be expressed by the displacement of any point on it, for example the displacement of center point 𝐶 1 .

To express the rigid body displacement, the Small Displacement Torsor (SDT) theory is used.

The concept of SDT was developed by Bourdet and Clément to fit geometrical surfaces to point clouds [START_REF] Bourdet | Controlling a complex surface with a 3 axis measuring machine[END_REF][START_REF] Bourdet | The concept of the small displacement torsor in metrology[END_REF]. In the field of metrology and tolerancing, the displacement is usually very small, and can be linearized. Using the SDT, the displacement of surface 𝑆 1 expressed on center point 𝐶 1 is:

(a) (b) 𝑆 2 𝑆 1 𝐶 1 𝐶 2 𝐴 1 𝐴 2 𝑑 𝒏 𝟐 𝒏 𝟏 𝑆 2 ′ 𝑆 1 ′ 𝐶 1 ′ 𝐴 1 ′ 𝐶 2 ′ 𝐴 2 ′ 𝑑 ′ 𝒏 𝟐 ′ 𝒏 𝟏 ′ [𝑺𝑫𝑻 𝟏 ] = [ 𝒓 𝟏 𝒕 𝟏,𝑪𝟏 ] 𝐶1 = [ 𝛼 1 𝑢 1,𝐶1 𝛽 1 𝑣 1,𝐶1 𝛾 1 𝑤 1,𝐶1 ] 𝐶1 (40)
The subscribe 𝐶1 indicates the displacement expressed at center point 𝐶 1 . 𝑺𝑫𝑻 𝟏 contains two components 𝒓 𝟏 and 𝒕 𝟏,𝑪𝟏 . 𝒓 𝟏 is the rotation of the rigid body, which is the same for all vertices on surface 𝑆 1 . Rotation 𝒓 𝟏 contains three rotation angles 𝛼, 𝛽 and 𝛾, which correspond to the rotations around the 𝑥, 𝑦 and 𝑧 axes respectively. 𝒕 𝟏,𝑪𝟏 is the translation of point 𝐶 1 . It contains 𝑢 1,𝐶1 , 𝑣 1,𝐶1 and 𝑤 1,𝐶1 which are translations of 𝐶 1 along the 𝑥, 𝑦 and 𝑧 axes.

For any vector attached to the surface, after rigid body displacement, its direction changes.

Take vector 𝒏 𝟏 as an example, after displacement 𝒏 𝟏 ′ could be calculated as:

𝒏 𝟏 ′ = 𝒏 𝟏 -𝒏 𝟏 × 𝒓 𝟏 (41)
The translation of any other vertex on surface 𝑆 1 , for example vertex 𝐴 1 , could also be expressed using 𝑺𝑫𝑻 𝟏 on center point 𝐶 1 . Assume the translation of vertex 𝐴 1 is 𝒕 𝟏,𝑨𝟏 , then:

𝒕 𝟏,𝑨𝟏 = 𝒕 𝟏,𝑪𝟏 + 𝐴 1 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ × 𝒓 𝟏 (42)
where 𝐴 1 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ is the vector pointing from 𝐴 1 to 𝐶 1 . The displacement of surface 𝑆 2 could be expressed by Equations ( 40), ( 41) and ( 42) in a similar way. The use of SDT benefits the linearization of distance between surfaces, which will be described in the following section.

Distance between Surfaces

Based on the rigid body assumption, there should be no interpenetration between part surfaces. This can be expressed by the constraints on surface distance.

For a vertex 𝐴 1 on surface 𝑆 . Then, using Equations ( 41) and ( 42), the distance after displacement is 𝑑 ′ :

𝑑 ′ = 𝐴 2 ′ 𝐴 1 ′ ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ • 𝒏 2 ′ = (𝐴 2 𝐴 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ + 𝒕 𝟏,𝐀𝟏 -𝒕 𝟐,𝐀𝟐 ) • (𝒏 2 -𝒏 2 × 𝒓 2 ) = (𝐴 2 𝐴 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ + 𝒕 𝟏,𝑪𝟏 + 𝐴 1 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ × 𝒓 𝟏 -𝒕 𝟐,𝑪𝟐 -𝐴 2 𝐶 2 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ × 𝒓 𝟐 ) • (𝒏 2 -𝒏 2 × 𝒓 2 ) = 𝐴 2 𝐴 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ • 𝒏 2 + 𝒕 𝟏,𝑪𝟏 • 𝒏 2 + 𝐴 1 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ × 𝒓 𝟏 • 𝒏 2 -𝒕 𝟐,𝑪𝟐 • 𝒏 2 -𝐴 2 𝐶 2 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ × 𝒓 𝟐 • 𝒏 2 -𝐴 2 𝐴 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ • (𝒏 2 × 𝒓 2 ) -(𝒕 𝟏,𝑪𝟏 + 𝐴 1 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ × 𝒓 𝟏 -𝒕 𝟐,𝑪𝟐 -𝐴 2 𝐶 2 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ × 𝒓 𝟐 ) • (𝒏 2 × 𝒓 2 ) (43)
In tolerance analysis, the rotation and translation are usually very small, and their product is second order small, which can be ignored in a linear approximation. Therefore the last term in Equation ( 43) is not considered, and after simplification, the distance 𝑑 ′ is:

𝑑 ′ = 𝒏 2 • 𝒕 𝟏,𝑪𝟏 + 𝒏 2 × 𝐴 1 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ • 𝒓 𝟏 -𝒏 2 • 𝒕 𝟐,𝑪𝟐 + 𝐴 1 𝐶 2 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ × 𝒏 2 • 𝒓 𝟐 + 𝑑 = [ 𝒏 2 × 𝐴 1 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ 𝒏 2 𝐴 1 𝐶 2 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ × 𝒏 2 -𝒏 2 ] • [ 𝒓 𝟏 𝒕 𝟏,𝑪𝟏 𝒓 𝟐 𝒕 𝟐,𝑪𝟐 ] 𝑇 + 𝑑 = 𝒂 • 𝑺𝑫𝑻𝒔 + 𝑑 (44)
where 𝒂 is the vector containing information based on the initial positions of skin model shapes. 𝑺𝑫𝑻𝒔 is the vector of SDT for all the skin model shapes. SDT is expressed at the center point of the skin model shape. For the discrete skin model shapes, there are many vertices on the surfaces and many corresponding distances. Equation ( 44) is organized in matrix form as:

𝒅 ′ = 𝑨 • 𝑺𝑫𝑻𝒔 + 𝒅 (45)
The local distance 𝒅 ′ is calculated using vector products.

Displacement Condition

Non-interpenetration and Internal Reaction Forces

In Equations ( 44) and ( 45), if the value of 𝑑 ′ is less than zero, this means the vertex 𝐴 1 is inside surface 𝑆 2 . Therefore, to avoid interpenetration, a constraint on distance is defined:

𝒅 ′ ≥ 𝟎 (46)
Due to the non-interpenetration condition, there are only two configurations between vertices: gap or contact. This is illustrated in Figure 4-13. The gap configuration is shown in Figure 4-13(a), where the distance between two vertices 𝐴 1 and 𝐴 2 is greater than zero. The contact configuration can be seen in Figure 4-13(b), where 𝐴 1 and 𝐴 2 are in the same position and the distance is equal to zero.

𝑆 2 𝑆 1 𝐶 1 𝐶 2 𝐴 1 𝐴 2 𝑆 1 𝑆 2 𝐶 1 𝐴 1 𝐶 2 𝐴 2 𝒗𝒓𝒇 𝟏 𝒏 𝟏 𝒏 𝟐 𝒗𝒓𝒇 𝟐
When the two surfaces are in contact, there could be reaction forces applied to the contact vertices. For example, in Figure 4-13(b), reaction forces 𝒗𝒓𝒇 𝟏 and 𝒗𝒓𝒇 𝟐 are applied to vertices 𝐴 1 and 𝐴 2 respectively. They are of the same magnitude, but in opposing directions. Let 𝒏 𝟏 and 𝒏 𝟐 be their directions, 𝑟𝑓 1 and 𝑟𝑓 2 be their magnitudes, then: 𝒗𝒓𝒇 𝟏 = 𝑟𝑓 1 • 𝒏 𝟏 and 𝒗𝒓𝒇 𝟐 = 𝑟𝑓 2 • 𝒏 𝟐 .

The signs of 𝑟𝑓 1 and 𝑟𝑓 2 are positive if the force directions are pointed into their material (inside the part). Because the adhesion force between parts is not considered, the value of the reaction forces should be positive. If we consider that the non-contact vertices have the value of reaction forces equal to zero, then every vertex has a reaction force greater than or equal to zero:

𝒓𝒇 ≥ 𝟎 (47)
where 𝒓𝒇 is the vector containing reaction force magnitudes for all vertices.

Displacement Boundary Condition

In order to constrain the displacements of parts, the concept of support is introduced, derived from supports in FEA. Supports are used either to constrain possible sliding of parts due to the absence of friction (example of  Reaction forces of support: the forces provided by the support to carry out the constraints on displacement. In the case in Figure 4-2(b), the reaction force of support corresponds to the static friction force stopping the part from sliding.

Depending on the specific assembly cases, forces may or may not be considered. Any boundary condition that contains constraints on displacement and reaction forces could be modeled using support, such as the assembly fixture and manual operations. Therefore, flexible use of support could simplify the modeling of the assembly simulation.

For the contact between parts and support, there is only one configuration. The part should always be in contact with the support, as can be seen from Once the support is introduced into the simulation, to maintain the contact, constraints on part displacement are added at the same time. In Equation ( 40), the displacement of the part is expressed by six torsor variables. To constrain the displacement, the following methods can be used:

 Set one or several torsor variable values equal to zero.

 Establish linear constraint between torsor variables.

The linear constraints between torsor variables take the form:

{ 𝒗 𝟏 • 𝒓 𝟏 = 0 𝒗 𝟐 • 𝒕 𝟏,𝑪𝟏 = 𝟎 (48)
Where 𝒗 𝟏 is the vector of rotation axis that will be constrained, and 𝒓 𝟏 is the corresponding rotation torsor variable. 𝒗 𝟐 is the direction vector of translation that will be constrained, and 𝒕 𝟏,𝑪𝟏 is the corresponding translation torsor variable.

In application, if the torsor variable equals zero, the corresponding column in the matrix is deleted. If a linear constraint between torsor variables is established, an additional row is added to the matrix to represent this constraint.

Balance of Forces and Moments

The importance of assembly loads has been well explained in previous sections. They are therefore introduced as assembly load boundary conditions.

The assembly loads applied to parts are generally considered to be force and torque. The above examples consider only one force load applied to one part, whereas more forces or torques are possible. Therefore during the assembly process, there are external forces 𝒇, external torques

𝑆 2 𝐶 2 𝐴 2 𝒗𝒓𝒇 𝟏 𝒗𝒓𝒇 𝟐 Support 𝑆 2 𝐶 2 𝐴 2 𝒗𝒓𝒇 𝟐 Support 𝑆 2 𝐶 2 𝐴 2 𝒗𝒓𝒇 𝟏 Support (a) (b) (c)
𝒕𝒒, internal or support reaction forces 𝒓𝒇. For each skin model shape, after assembly they should find a balanced position under internal and external loads. This balance can be expressed as:

{ ∑ 𝑟𝑓 𝑖,𝑗 • 𝒏 𝒊,𝒋 + ∑ 𝑓 𝑖,𝑘 • 𝒏 𝒊,𝐤 = 𝟎 ∑ 𝑳 𝒊,𝒋 × (𝑟𝑓 𝑖,𝑗 • 𝒏 𝒊,𝒋 ) + ∑ 𝑳 𝒊,𝒋 × (𝑓 𝑖,𝑘 • 𝒏 𝒊,𝐤 ) + ∑ 𝑡𝑞 𝑖,𝑞 • 𝒏 𝒊,𝐪 = 𝟎 (49)
Where:

 𝑟𝑓 𝑖,𝑗 is the 𝒋-th reaction force applied to the 𝒊-th skin model shape, 𝒏 𝒊,𝒋 is its direction. The reaction forces contain forces between contacting parts, and the forces from supports.

 𝑓 𝑖,𝑘 is the 𝒌-th external force applied to the 𝒊-th skin model shape, 𝒏 𝒊,𝒌 is its direction. 𝑡𝑞 𝑖,𝑞 is the 𝒒-th external torque, and 𝒏 𝒊,𝒒 is its axis.

 The balance of moment also considers the moments caused by external forces and internal reaction forces.

Multiple Assembly Parts

The assembly usually uses two or more skin model shapes at the same time. Generally, one of the skin model shapes should be fixed. This depends on two aspects: the effect of support and the assembly sequence.

The effect of support has been explained earlier. The displacements should be constrained and reaction forces of the support should be considered.

For multi-stage assembly, the assembly sequence influences the precision of the final product.

To consider the assembly sequence, the following simulation method is applied:

 For each assembly station, two or more subassemblies are joined together. One of them is considered as fixed. The displacements of other subassemblies can be adjusted.

 After assembly, the relative displacement between these subassemblies is fixed. They are then treated as a new subassembly in the next assembly station.

The assembly should therefore be decomposed into subassemblies by the user according to the assembly sequence.

Formulation of the Assembly Problem

In the above analysis, the factors that influence assembly have been transformed into linear constraints, which can be seen from Equations ( 46), ( 47) and [START_REF] Raffin | Curvilinear displacement of free-form-based deformation[END_REF]. Thereafter, the assembly simulation of skin model shapes consists in actually finding the displacement 𝑺𝑫𝑻𝒔 and reaction forces 𝒓𝒇 that satisfy these constraints.

Notice the two configurations in Figure 4-13, if the parts are not in contact, then distance 𝑑 ′ is greater than zero, and there is no reaction force. If they are in contact, 𝑑 ′ equals zero and the reaction force is positive. For the contact between part and support in Figure 4-14, the distance 𝑑 ′ is always equal to zero and the reaction forces can be positive. These cases could be concluded as:

{ 𝑑 ′ ≥ 0 𝑟𝑓 = 0 and { 𝑑 ′ = 0 𝑟𝑓 ≥ 0 (50)
Moreover, for all the potential contact point pairs, Equation ( 50) can be expressed in the form of vector product:

𝒓𝒇 𝑻 𝒅 ′ = 𝟎 (51)
Equation ( 51) is called the LCC (Linear Complementarity Condition), which has been widely used in the simulation of multi-body dynamics [START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF][START_REF] Williams | RPI-MATLAB-Simulator: A Tool for Efficient Research and Practical Teaching in Multibody Dynamics[END_REF]. Based on the linearized expression of distance 𝒅 ′ in Equation ( 45), distance 𝒅 ′ in Equation ( 51) is replaced by 𝑺𝑫𝑻𝒔:

𝒓𝒇 𝑇 𝒅 ′ = 1 2 [ 𝒅 ′ 𝒓𝒇 ] 𝑇 [ 𝒓𝒇 𝒅 ′ ] = 1 2 [ 𝒅 ′ 𝒓𝒇 ] 𝑇 [ 𝟎 𝑰 𝑰 𝟎 ] [ 𝒅 ′ 𝒓𝒇 ] = 1 2 ([ 𝑨 𝟎 𝟎 𝑰 ] [ 𝑺𝑫𝑻𝒔 𝒓𝒇 ] + [ 𝒅 𝟎 ]) 𝑇 [ 𝟎 𝑰 𝑰 𝟎 ] ([ 𝑨 𝟎 𝟎 𝑰 ] [ 𝑺𝑫𝑻𝒔 𝒓𝒇 ] + [ 𝒅 𝟎 ]) = 1 2 [ 𝑺𝑫𝑻𝒔 𝒓𝒇 ] 𝑇 [ 𝟎 𝑨 𝑇 𝑨 𝟎 ] [ 𝑺𝑫𝑻𝒔 𝒓𝒇 ] + [ 𝟎 𝒅 𝑇 ] [ 𝑺𝑫𝑻𝒔 𝒓𝒇 ] (52)
Let 𝑿 = [ 𝑺𝑫𝑻𝒔 𝒓𝒇 ], then LCC in Equation ( 51) is simplified as:

1 2 𝑿 𝑻 [ 𝟎 𝑨 𝑇 𝑨 𝟎 ] 𝑿 + [ 𝟎 𝒅 𝑇 ]𝑿 = 𝟎 (53)
The LCC indicates the existence of configurations in which the skin model shapes could be assembled. If the skin model shapes could be assembled, then 𝑿 satisfies all these constraints in Equations ( 46), ( 47), ( 49) and ( 53):

{ 𝒓𝒇 ≥ 0 𝑨 ⋅ 𝑺𝑫𝑻𝒔 + 𝒅 ≥ 0 ∑ 𝑟𝑓 𝑖,𝑗 • 𝒏 𝒊,𝒋 + ∑ 𝑓 𝑖,𝑘 • 𝒏 𝒊,𝐤 = 𝟎 ∑ 𝑳 𝒊,𝒋 × (𝑟𝑓 𝑖,𝑗 • 𝒏 𝒊,𝒋 ) + ∑ 𝑳 𝒊,𝒋 × (𝑓 𝑖,𝑘 • 𝒏 𝒊,𝐤 ) + ∑ 𝑡𝑞 𝑖,𝑞 • 𝒏 𝒊,𝐪 = 𝟎 1 2 𝑿 𝑻 [ 𝟎 𝑨 𝑇 𝑨 𝟎 ] 𝑿 + [ 𝟎 𝒅 𝑇 ]𝑿 = 𝟎 (54)
On the other hand, if the skin model shapes cannot be assembled, any 𝑿 that satisfies Equations ( 46), ( 47), [START_REF] Raffin | Curvilinear displacement of free-form-based deformation[END_REF] will result in the left side of Equation ( 53) having a positive value. Thus the LCC is not satisfied.

To this end, minimizing the left side of the LCC corresponds to finding the assembly configuration 𝑿. The assembly problem is transformed into a quadratic optimization problem, and linear constraints in Equations ( 46), ( 47) and ( 49) should be satisfied:

𝑓𝑖𝑛𝑑 𝑿 = [ 𝑺𝑫𝑻𝒔 𝒓𝒇 ] 𝑚𝑖𝑛 ( 1 2 𝑿 𝑻 [ 𝟎 𝑨 𝑇 𝑨 𝟎 ] 𝑿 + [ 𝟎 𝒅 𝑇 ]𝑿) 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 { 𝒓𝒇 ≥ 0 𝑨 ⋅ 𝑺𝑫𝑻𝒔 + 𝒅 ≥ 0 ∑ 𝑟𝑓 𝑖,𝑗 • 𝒏 𝒊,𝒋 + ∑ 𝑓 𝑖,𝑘 • 𝒏 𝒊,𝐤 = 𝟎 ∑ 𝑳 𝒊,𝒋 × (𝑟𝑓 𝑖,𝑗 • 𝒏 𝒊,𝒋 ) + ∑ 𝑳 𝒊,𝒋 × (𝑓 𝑖,𝑘 • 𝒏 𝒊,𝐤 ) + ∑ 𝑡𝑞 𝑖,𝑞 • 𝒏 𝒊,𝐪 = 𝟎 (55)
The solving of the optimization problem uses quadratic optimization algorithms in mathematical packages, such as the interior point method.

Overview of Simulation Process

Using the formulated quadratic optimization problem, the assembly simulation is conducted according to the following steps:

Application and Analysis

Generating Skin Model Shapes

In the application of the proposed simulation method, skin model shapes with different levels of manufacturing defects are generated first. In the work of Zhang et al. [START_REF] Zhang | Discrete shape modeling for skin model representation[END_REF], a statistical shape analysis-based method was used to simulate skin model shapes. Schleich et al. [START_REF] Schleich | Skin Model Shapes: A new paradigm shift for geometric variations modelling in mechanical engineering[END_REF] used the random field method to simulate the form defects. There are also other methods available [START_REF] Formosa | Modal expression of form defects[END_REF][START_REF] Wilma | Manufacturing Signature for Tolerance Analysis[END_REF][START_REF] Cao | Geometrical Simulation of Multiscale Toleranced Surface With Consideration of the Tolerancing Principle[END_REF][START_REF] Homri | Tolerance analysis --Form defects modeling and simulation by modal decomposition and optimization[END_REF][START_REF] Corrado | Manufacturing signature and operating conditions in a variational model for tolerance analysis of rigid assemblies[END_REF] to simulate manufacturing defects. The method used here is based on the work in Chapter 0, which is explained by Figure 2-1.

Mating Two Squares with Manufacturing Defects

To validate the simulation result of the proposed method, a simple example is used where the influences of assembly forces can be observed. The load balance of the mating is also checked.

This example corresponds to the case shown in Figure 4 To validate whether the simulation result is feasible or not, we can fix the contact points and calculate reaction forces at these points from the balance constraints in Equation [START_REF] Raffin | Curvilinear displacement of free-form-based deformation[END_REF]. Because there is no adhesion force, the calculated forces should be positive. Considering the top surface, the calculation results are shown in Table 4-

Assembly Examples

The example in Figure 4-16 uses surfaces with large form defects, where it is easy to check the simulation result visually. In the following simulations of 3D parts, realistic defects are considered.

The assembly of a cube in the corner is shown in Figure 4 Different from the typical 3-2-1 positioning, for each contact plane pair in the figure, they have two contact points. This emphasizes the fact that depending on detailed form defects and assembly loads, the positioning of the part could be different from the simplified assumptions.

In this assembly of a cube, the mating and load conditions are simple. Thus the case in Figure 4-11 is used to show the consideration of system equilibrium. 3D mesh models are shown in Figure 4-18. The model contains a cube with a hole, an L-shaped part and a pin. The pin passes through the L-shaped part and is inserted into the hole in the cube, and the cube is fixed. In simulations where loads are not considered, the position of the pin is floating inside the hole, which provides infinite solutions. Among these solutions, only certain configurations are close to the real assembly.

To determine these configurations, load boundary conditions are necessary. For each part, gravity is applied at center point. An assembly force is applied to the bottom of the pin along its axis direction.

First, a large value is given to the assembly force. After assembly, the three parts are in contact at points which are highlighted in Figure 4-19(a, b). In a second case, the assembly force is decreased, and the assembly result is shown in Figure 4-19(c, d). The difference is obvious. At first when the force is large enough, mating surfaces between these parts are tightly in contact. When the force is decreased, to get balance, the pin and L-shaped part start to rotate. In the final position, the pin has two contact regions with the hole in the cube, which helps to balance the system. For the case in Figure 4-16(a), if we change the position of the external assembly force slightly, for instance, by moving it to the position of the blue arrow in Figure 4-20(a), the three contact points do not change. This means it is in a stable configuration.

In Figure 4-20(b), the contact points for each surface are located very close together.

Therefore, although the assembly is balanced under current load conditions, it is not stable. A slight change in the load value or direction may lead to a change in assembly configuration.

In the simulation result in Figure 4-19(c) and (d), the L-shaped part has a rotation around the axis of its hole. This is due to the mesh size of the cylinder surfaces, and the balance of the force and moment. The section view in Figure 4-21 is used to explain the phenomenon.

In Figure 4-21, both the hole and pin are discrete and expressed by segments. The pin is fixed and the hole could translate and rotate under the force applied to its center point. Besides, the configuration is not sufficiently stable.

Simulation Analysis

In this section, the proposed assembly simulation method is applied to different types of model. The importance of manufacturing defects and assembly boundary conditions (load and displacement) are emphasized by these cases:

 From the surface mating case in Figure 4-16, the influence of assembly load on contact configuration is highlighted. The magnitude of reaction forces and contact regions are checked to validate the simulation result.  The assembly simulation is a generalization and extension of the surface mating. For the case in Figure 4-17, the mating condition and load condition are more complex.

However, the non-interpenetration constraints and balance of loads are still tenable.

The calculation amount is increased, but the simulation method and process remain the same.

 Assembly of three parts in Figure 4-19 introduced cylindrical mating pairs. Different assembly results are generated in consideration of the assembly loads.  The hand saw example shows the ability of the proposed method to handle complex simulations. Moreover, a statistical assembly simulation could be conducted. A batch of skin model shapes are generated in advance, then they are selected at random and deviations are evaluated after assembly.

In addition, the stability of assembly simulation results is analyzed. These unstable results are decided by form defects, position and magnitude of assembly loads, and the precision of meshes.

Conclusion

In this chapter, we studied tolerance analysis taking manufacturing defects into consideration.

First, the issue of assembly with form defects is discussed. Differences in the assembly of nominal model, ideal surface-based deviation model and skin model shape are highlighted. The analysis shows that assembly boundary conditions (load and displacement) must be considered when form defects are introduced.

Next, the assembly simulation method is proposed. The distance between potential contact point and plane is linearized by SDT. The balance of the system under assembly loads is modelled by LCC. The proposed method considers form defects and assembly loads at the same time, while the dynamic process is ignored to simplify the calculation.

Simulation of different assembly examples is conducted in the last section. From simple to complex, the effectiveness of the proposed method is proved. Meanwhile, the contact regions and reaction forces are also checked to guarantee that the simulation result is feasible.

The simulation results highlight two challenges for the proposed methods. The first is due to the random generation of the manufacturing defects: assemblies may be impossible if the defects are too large. As a result, assemblies must be checked before simulation itself. However, before this, the automatic generation of the defects must be under control and skin model shapes must verify assembly specifications.

The second challenge concerns computation time, knowing that time increases with the number of parts, meshes, and the number of assemblies. In fact, about 70% of the calculation time is spent searching for the nearest potential contact points. To reduce the search time, the algorithm and data structure could be improved. To reduce the time for quadratic optimization, the problem could be reduced by considering only the closest potential contact points.

The current work presents a general way to conduct tolerance analysis with skin model shapes.

Deviation Evaluation

In this chapter, measurement uncertainties and GeoSpelling are introduced. The Small Displacement Torsor-based deviation evaluation method is explained and improved. Using this method, virtual metrology can be conducted on a single skin model shape or a skin model shape assembly. Additionally, the deviation evaluation method is customized and integrated into an online Virtual Laboratory for e-learning.

Introduction

In the previous chapters, skin model shape generation and assembly simulation have been studied. This chapter presents the evaluation of characteristics on the isolated skin model shape or the skin model shape assembly. With this evaluation, virtual metrology and evaluation of the functional characteristics of an assembly can be carried out. Virtual metrology could be used during the design of the product, to prepare the metrology process before any part is manufactured. It is also integrated into an online Virtual Laboratory (VL) for e-learning.

The evaluation of geometrical deviation is not limited to deciding whether the manufactured part can be accepted or not. Indeed, it provides a basis for the error source diagnosis of multistage manufacturing, and helps to improve the manufacturing processes. To give information about the precision of an assembly at the last stage of assembly simulation, an evaluation is made.

Depending on the evaluation objective and the available equipment, different measurement strategies could be planned. Specialized instruments like gauges are used to check whether the part deviation is acceptable or not. Using general measurement instruments, like a caliper and micrometer, the deviation value can be evaluated. With the advances in measurement instruments, more general tools (e.g. Coordinate Measurement Machine (CMM), 3D camera, white light scanner, laser scanner, 3D Computed Tomography (CT), etc.) are available and used in industry [START_REF] Schmitt | Advances in Large-Scale Metrology -Review and future trends[END_REF]. To accompany these tools, software-based virtual gauge, statistical analysis package and error compensation have been developed. For these general measurement instruments, the standardization of data exchange format and precision of algorithms inside software play critical roles.

In this chapter, the method used for deviation evaluation is introduced. The concept of method uncertainty will be explained first. To avoid method uncertainty, the GeoSpelling-based tolerance specification approach is used. The application examples of GeoSpelling are given using the skin model shapes previously generated. Within the framework of GeoSpelling, an SDT-based algorithm is used to carry out some of the operations and evaluate the deviations. The algorithm is also modified and integrated into a Virtual Laboratory (VL) for e-learning.

Method Uncertainty and GeoSpelling

Uncertainties are inherent in manufacturing and measurement processes and they influence all stages of the product lifecycle. For measurement, it is easy and intuitive to realize that there are uncertainties due to the measurement equipment. The solution for this type of uncertainty relies on improving equipment precision. However, this does not mean that the problem is solved.

There are more kinds of uncertainties, which may have significant influences on the measurement result. In the standard ISO 17450-2 [START_REF]Geometric Product Specification-General concepts-Part 2: Basic tenets, specifications, operators, uncertainties and ambiguities[END_REF], the ambiguity of specification and method are addressed.

The impact of ambiguities between function description and specification, method and implementation, are emphasized and analyzed. In the work of Mathieu and Ballu [START_REF] Mathieu | GEOSPELLING: a common language for Specification and Verification to express Method Uncertainty[END_REF], a model is used to express the relation between these uncertainties and design/verification processes. The detailed measurement and evaluation methods corresponding to specifications are not defined, and these methods may vary depending on the experience of the metrologist and the available equipment. Lastly, the uncertainties introduced by the measurement equipment in implementation also contribute to the uncertainty of verification.

It can be clearly seen that besides the uncertainty caused by measurement equipment, there are many specification-related uncertainties in upstream processes. These accumulate from the early stages, and may lead to a bias in manufacturing and verification objective. Some current research is trying to solve the problem of specification uncertainty. In the work of Sarigecili et al. [START_REF] Sarigecili | Interpreting the semantics of GD&T specifications of a product for tolerance analysis[END_REF], the semantics of GD&T are interpreted and mapped to an ontology-based model. Since the interpretation system is designed by related experts, there may be no Designers are still puzzling over the definition of specification.

In fact, the problem of specification uncertainty and method uncertainty was noticed a long time ago [START_REF] Ballu | The processing of measured points in coordinates metrology in agreement with the definition of standardized specifications[END_REF]. To avoid these uncertainties from the very beginning, a coherent and comprehensive language was developed to express the semantics of specifications in various stages. This language is called "GeoSpelling" and was initially proposed by Ballu and Mathieu [START_REF] Ballu | Analysis of dimensional and geometrical specifications: standards and models[END_REF].

GeoSpelling is based on a limited number of operations. The design of GeoSpelling considered the semantic coherence, the demand for digitalization and automation, the readability for both humans and machines, and the learning cost for engineers. Due to these characteristics, the concepts developed in GeoSpelling have been integrated into ISO standard 17450 [START_REF]Geometric Product Specification-General concepts-Part 1: Model for geometrical specification and verification n[END_REF][START_REF]Geometric Product Specification-General concepts-Part 2: Basic tenets, specifications, operators, uncertainties and ambiguities[END_REF].

Figure 5-2 Operations of GeoSpelling [START_REF] Ballu | Formal Language for GeoSpelling[END_REF].

In the work of Mathieu and Ballu [START_REF] Mathieu | A Model for a Coherent and Complete Tolerancing Process[END_REF][START_REF] Mathieu | GEOSPELLING: a common language for Specification and Verification to express Method Uncertainty[END_REF], the problem of specification uncertainty is discussed in detail, and cases are studied to prove the effectiveness of a GeoSpelling-based expression approach. In the work of Dantan et al. [START_REF] Dantan | Geometrical product specifications -model for product life cycle[END_REF], GeoSpelling is used to express specifications covering product lifecycle. Both Dantan et al. [START_REF] Dantan | Geometrical Specification Model for Gear -Expression, Metrology and Analysis[END_REF] and Zhang et al. [START_REF] Zhang | Geometric Product Specification of Gears: The GeoSpelling Perspective[END_REF] have applied GeoSpelling to specification expression and tolerance verification. Later, in the work of Ballu et al.

Examples of GeoSpelling

The graphical standardized GD&T system on a 2D/3D model is currently the common geometrical specification expression method used in technical documentation. To have a comparison between this traditional method and a GeoSpelling-based method, we will use them to express the same specifications on the model.

The example is the position tolerance of two cylindrical surfaces. First, the graphical standardized GD&T expression is illustrated. The nominal model and specifications can be seen in or should be, defined in these standards. However, the standards are not complete and not always relevant. We first explain the specification in natural language, see below:

 The non-ideal surface S1 is identified by partition from the skin model shape. S1 corresponds to the nominal plane NPL1.

 The ideal plane PL1 is associated to surface S1:

o Constraints: minimum material distance between PL1 and S1 is greater than zero (PL1 is outside the material).

o Objective: minimize the maximum material distance between PL1 and S1.

 The non-ideal surface S2 is identified by partition from the skin model shape. S2 corresponds to the nominal cylinder NCY2.

 The ideal cylinder CY2 is associated to surface S2:

o Constraints: minimum material distance between S2 and CY2 is greater than zero (CY2 is outside the material). The axis of CY2 is CY2a. The angle between CY2a and PL1 is 90 deg.

o Objective: maximize the diameter of ideal cylinder CY2, which is Dia(CY2).

 The non-ideal surfaces S3 and S4 are identified by partition from the skin model shape, which correspond to nominal NCY3 and NCY4 respectively.

 The extracted median lines for S3 and S4 are calculated, and named EML3 and EML4.

 A straight line SL5 corresponding to the axis of cylinder tolerance zone is constructed:

Deviation Evaluation Based on SDT

To evaluate the deviation of measurement data (such as point cloud or mesh), an ideal feature is associated to the non-ideal surfaces. This association process is generally transformed into an optimization problem, which minimizes distance between ideal feature and non-ideal surfaces, or minimizes other characteristics. During association, non-linear constraints should be considered to represent 3D spatial relations between features.

To take these non-linear constraints into consideration and solve the problem efficiently, the SDT-based method was developed [START_REF] Bourdet | The concept of the small displacement torsor in metrology[END_REF]. The constraints and optimization objective are linearized using the SDT, thus the problem could be solved efficiently by linear programming algorithms. In the rest of this section, the linearization method and then formulation of the problem will be introduced.

Elementary Surfaces and Invariances

In coordinate metrology, the measured surfaces are saved as point clouds. Engineers need to associate ideal features (like plane, cylinder, and sphere) to these point data, and the deviation values from points to the corresponding ideal features are evaluated.

Based on the TTRS model [START_REF] Desrochers | A dimensioning and tolerancing assistance model for CAD/CAM systems[END_REF], the ideal features are generalized and classified into seven classes of elementary surfaces, as shown in Table 5-3. These surfaces are the most commonly used ones in mechanical parts, and mating them formed mechanical joints, which related to the functional requirements.

The classification of elementary surfaces is also based on the degree of invariance [START_REF]Geometric Product Specification-General concepts-Part 1: Model for geometrical specification and verification n[END_REF]. For any ideal surface, its degree of invariance corresponds to the rotations or translations that keep the surface unchanged. For example, the ideal cylindrical surface will not change when it rotates or translates along its axis. However, for general surfaces, any rotation or translation could change it. The degrees of invariance for elementary surfaces are described in Table 5-3.

For each surface belonging to the seven types of elementary surfaces, their position and orientation can be determined by the Minimum Geometric Datum Elements (MGDE) associated to them. Point, straight line and plane are considered as "basic MGDEs" because other MGDEs could be generated by combining them. For basic MGDEs, they have degrees of invariance that are easy to recognize. For the combination of basic MGDEs, certain rules were proposed by Desrochers et al. [START_REF] Desrochers | A dimensioning and tolerancing assistance model for CAD/CAM systems[END_REF] to establish the appropriate new MGDE and its degrees of invariance. In the second column of Table 5-3, the MGDE for each surface is shown in orange, and the third column of the table indicates which MGDE is used.

Express Degree of Invariance by SDT

In associating an ideal surface to a point cloud, we usually fix the position of the point cloud, and rotate or translate the ideal surface to minimize the distance between them. Due to the degrees of invariance of the ideal surfaces, some rotation or translation may not be helpful for minimization. For example, in the case of fitting a cylinder to a point cloud, rotation and translation along the cylinder axis will not influence the distance between cylinder and points.

Meanwhile, this means that there is an infinity of solutions which minimize the distance and thus give the same value for this distance. To avoid these unnecessary rotations and translations, constraints are added to the displacements along/around directions of degrees of invariance.

𝒏 𝟐𝟏 ′ is the same vector after translation and rotation. Based on the small displacement assumption, we have:

𝒏 𝟐𝟏 ′ = 𝒏 𝟐𝟏 -𝒏 𝟐𝟏 × 𝒓 𝟐 (61)
Then distance 𝑑 ′ can be written as:

𝑑 ′ = 𝐶 2 ′ 𝐶 1 ′ ⃗⃗⃗⃗⃗⃗⃗⃗⃗ • 𝒏 𝟐𝟏 ′ = (𝐶 2 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ + 𝒕 𝟏,𝑪𝟏 -𝒕 𝟐,𝑪𝟐 ) • (𝒏 𝟐𝟏 -𝒏 𝟐𝟏 × 𝒓 𝟐 ) = 𝐶 2 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ • 𝒏 𝟐𝟏 + 𝒏 𝟐𝟏 • 𝒕 𝟏,𝑪𝟏 -𝒏 𝟐𝟏 • 𝒕 𝟐,𝑪𝟐 -𝐶 2 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ × 𝒏 𝟐𝟏 • 𝒓 𝟐 (62)
b) When the distance between point and straight line is close to zero, e.g. 𝑑 ′ = 0, then:

𝟎 = 𝐶 2 ′ 𝐶 1 ′ ⃗⃗⃗⃗⃗⃗⃗⃗⃗ × 𝒏 𝟐 ′ = 𝐶 2 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ × 𝒏 𝟐 + 𝒕 𝟏,𝑪𝟏 × 𝒏 𝟐 -𝒕 𝟐,𝑪𝟐 × 𝒏 𝟐 -𝒏 𝟐 (𝐶 2 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ • 𝒓 𝟐 ) + 𝒓 𝟐 (𝐶 2 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ • 𝒏 𝟐 ) (63)

Formulating and Solving the Problem

During the association process, the ideal features are repositioned to fit the measured points.

The displacement of ideal features is expressed by small displacement torsors 𝑺𝑫𝑻𝒔, which are the variables we need to find.

The association of an ideal feature to a non-ideal surface is transformed into an optimization problem. The optimization objective is generally to:

 Minimize/maximize the distance between ideal feature and non-ideal surface. The distance could be unsigned distance, signed material distance, or quadratic distance.

 Minimize/maximize the inherent characteristic of ideal feature. For example, maximize the radius of a cylindrical feature to inscribe the cylindrical surface mesh.

Using the linearization method introduced in earlier sections, the objective function is expressed as the linear function 𝑶𝒃𝒋(𝑺𝑫𝑻𝒔). Its value depends on the displacement of ideal features.

There are several kinds of constraint between features during association:

 Constraints between ideal feature and non-ideal surface. These constraints are distance constraints, such that the signed material distance should be greater than zero.

 Constraints between ideal features. These constraints could be distance constraints or angle constraints. For example, the two ideal planes are parallel (angle equals zero), and the distance between them is 30mm.

These constraints are also linearized using SDT, and organized to have the general form:

𝑪𝒏𝒔𝒕(𝑺𝑫𝑻𝒔) = 𝟎 (64)
Therefore, the generalized linear programming problem is:

𝑓𝑖𝑛𝑑 𝑺𝑫𝑻𝒔 𝑚𝑖𝑛 𝑶𝒃𝒋(𝑺𝑫𝑻𝒔) 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 𝑪𝒏𝒔𝒕(𝑺𝑫𝑻𝒔) = 𝟎 (65)
For most of the cases displacement, objective and constraints are already linearized, and the problem could be solved using optimization algorithms efficiently (e.g. simplex, interior point, etc.).

For some special cases, however, non-linear expressions are required. Due to the linearization of the problem, the algorithm will find the local solution depending on initial positioning. Meanwhile, an iterative optimization and repositioning process is required to achieve a certain precision level.

Improvements on Evaluation Method

Based on the evaluation method described above, improvements are proposed to enhance its functions. With the first improvement it is possible to consider linear relations between characteristics, while with the second improvement the invariance of feature group can be calculated. These are explained in the following.

Add Linear Equations as Association Constraints

As introduced above, in our algorithm, constraints on characteristics (distance and angle) are linearized by SDT, therefore the association problem can be solved efficiently. However, these constraints are on a single characteristic, whereas there could be constraints between several characteristics. A basic example could be the equality between two distances to express symmetry.

In general, the most frequently used constraints between several characteristics can be expressed by a linear combination of basic characteristics. Two examples are shown below.

{ 𝑐 1 𝑑 1 ′ + 𝑐 2 𝑑 2 ′ = 𝐶 𝑐 1 𝛼 1 ′ + 𝑐 2 𝛼 2 ′ = 𝜃 (66)
It is difficult to associate feature group to measured points due to the complex constraints between features and datum. To handle this problem, several methods have been proposed.

Choley et al. [START_REF] Choley | A new variational association process for the verification of geometrical specifications[END_REF] associated complex geometry to measurement data with respect to datum system. Lehtihet et al. [START_REF] Lehtihet | On the composite position tolerance for patterns of holes[END_REF] studied the composite position tolerance of hole patterns using nonlinear formulation.

In coordinate metrology, two approaches are considered to evaluate the deviation of a feature group. The first approach considers that the members of the feature group (single features) may have relative displacements. The Theoretically Exact Dimensions (TED) [130] are explicitly expressed by constraints (in dimension and orientation) between members.

This could be explained by the example in Figure 5 Due to the constraints, the axes of the two cylinders are parallel, and the distance between them is the nominal distance 𝑑. For the second approach, the initial positions of features are extracted from the nominal model, and all the feature members are associated at the same time. The feature members have the same displacements thus the TED from the nominal model are conserved.

As can be seen from For both approaches, if the algorithm works well, they should generate the same association result. In applying the first approach, explicit constraints and situation features are needed to establish the relations between feature members. For a complex feature group, this procedure could be tedious and error-prone. Most of the time when we are doing measurements, the To use the second approach, the degree of invariance for the feature group is required. As explained before, the constraints on the degree of invariance could avoid unnecessary calculations during association. This invariance of feature group could be deduced in two ways:  Based on a table provided by Clément et al. [START_REF] Clément | Theory and Practice of 3-D Tolerancing for Assembly[END_REF],

Feature Group Association

 Based on the invariance of each feature member.

The methods that could be used to associate feature groups are shown in Figure 5678. The feature group invariance calculation method is flexible, and will be explained next.

Using SDT, feature displacement is linearized. The set of displacements of a feature according to one degree of invariance is represented by a line in a 6-dimensional real space ℝ 6 . The direction of this line can be described by an invariance direction vectors 𝒗. For the 𝑖-th single feature inside the feature group, its invariance directions according to all the invariance degrees form a basis:

𝐵 𝑖 = {𝒗 𝒊,𝟏 , 𝒗 𝒊,𝟐 , ⋯ , 𝒗 𝒊,𝒏 }, 𝑛 ≤ 6 (73)
where 𝒗 𝒊,𝒏 is the 𝑛-th invariance direction of the 𝑖-th feature. The basis 𝐵 𝑖 spans a subspace of ℝ 6 , which is represented as 𝕀 𝑖 .

Similarly, the sets of positional displacements of a feature are also lines in ℝ 6 , and have positional direction vectors 𝒗 ′ . The positional directions of the 𝑖-th feature form the basis 𝐵 𝑖 ′ :

𝐵 𝑖 ′ = {𝒗 𝒊,𝟏 ′ , 𝒗 𝒊,𝟐 ′ , ⋯ , 𝒗 𝒊,𝒎 ′ }, 𝑚 ≤ 6, 𝑚 + 𝑛 = 6 (74)
where 𝒗 𝒊,𝒎 ′ is the 𝑚-th positional direction of the 𝑖-th feature. The basis 𝐵 𝑖 ′ spans a subspace ℙ 𝑖 . The two subspaces are dual spaces [START_REF] Gerbino | How to investigate constraints and motions in assemblies by screw theory[END_REF][START_REF] Adams | Application of screw theory to constraint analysis of assemblies of rigid parts[END_REF], and they have the relation that:

𝕀 𝑖 = 𝐶𝑜𝑚𝑝𝑙𝑒𝑚𝑒𝑛𝑡(ℙ 𝑖 ) (75)
It should be noted that the invariance and positional directions of single features inside a feature group are expressed by SDT at a single point.

The above expressions are for single features. For a feature group, it also has invariance and positional directions and corresponding bases. These bases span the subspaces 𝕀 𝐹𝐺 and ℙ 𝐹𝐺 . 𝕀 𝐹𝐺 is the intersection of the invariance spaces of feature members, and ℙ 𝐹𝐺 is the Minkowski sum of the positional spaces of feature members:

{ 𝕀 𝐹𝐺 = ⋂ 𝕀 𝑖 𝑞 𝑖=1 ℙ 𝐹𝐺 = ⋃ ℙ 𝑖 𝑞 𝑖=1 (76)
Where 𝑞 is the number of feature members. Based on Equations ( 75) and ( 76), the invariance space of a feature group can be deduced from the positional spaces of feature members:

𝕀 𝐹𝐺 = 𝐶𝑜𝑚𝑝𝑙𝑒𝑚𝑒𝑛𝑡(ℙ 𝐹𝐺 ) = 𝐶𝑜𝑚𝑝𝑙𝑒𝑚𝑒𝑛𝑡(⋃ ℙ 𝑖 𝑞 𝑖=1 ) (77)
To calculate feature group invariance space 𝕀 𝐹𝐺 and its basis 𝐵 𝐹𝐺 , matrices are constructed corresponding to invariance and positioning. The invariance matrix 𝐼 𝑖 is size 𝑛 × 6, where each row inside the matrix corresponds to an invariance direction vector. The positional matrix 𝑃 𝑖 is size 𝑚 × 6, where each row represents a positional direction vector:

{ 𝐼 𝑖 = [𝒗
The sum operation in Equation ( 77) is conducted by connecting positional matrix 𝑃 𝑖 :

𝑃 𝐹𝐺 = [𝑃 1 , 𝑃 2 , ⋯ , 𝑃 𝑞 ] 𝑇 (79)
Using SDT, constraints on displacements in positional directions are written in the form of a linear equation group:

𝑃 𝐹𝐺 [ 𝑡
where 𝑡 𝑖,𝑃𝐹𝐺 in the column vector are torsor variables of the feature group, and the subscript index 𝑃𝐹𝐺 indicates the point at which the torsor is expressed. 𝟎 is 6 × 1 column zero vector.

The row vectors in 𝐼 𝐹𝐺 are the vectors which form the basis of the null space of 𝑃 𝐹𝐺 . To find invariance matrix 𝐼 𝐹𝐺 , the 𝑃 𝐹𝐺 is first transformed into row reduced echelon form by methods like SVD or Gaussian Pivot. After transformation, 𝑃 𝐹𝐺 * has the form:

𝑃 𝐹𝐺 * = [ 1 0 0 0 0 0 0 1 0 0 0 0 ⋮ 0 0 0 0 0 0 ] (81)
The linear equation group in Equation ( 80) is then: The user interface in the virtual metrology laboratory is shown in Figure 56789. Region 1 is the tool bar, which has three tabs including measurement instruments, objects and history. The measurement instruments are currently the digital caliper, micrometer and indicator with gauges.

Multiple functions of the caliper are integrated, measurement of inner size, outer size, shoulder and depth.

The Object tab provides parts which can be measured in the laboratory. These parts include skin model shapes. For teaching purposes, the deviations of the skin model shape can be amplified.

The skin model shapes contain manufacturing defects, therefore the measurement results vary, depending on the location of the surface measurement. This will help students to understand the uncertainties of measurement. Using skin model shapes with amplified deviations enables them to visualize the profile of form error. The History tab is a list of former operations, and these can be modified or repeated by users.

For example, comparisons could be made by changing certain operations and statistical analysis of the uncertainties can be carried out by repeating the operations. In future, a batch of parts could be measured by repeating these operations on different skin model shapes.

Region 2 is the information bar. Course aims, measurement projects, and real time response information are shown here. Data management is in the right-hand side of this bar.

Region 3 is the main operation area. It is a 3D interactive environment for the virtual laboratory, with objects and instruments imported from the tool bar (region 1). Movements of objects and measurements are conducted here.

A measurement example is shown in

Measurement of Outer Size

The function to measure outer size is designed for the caliper and the micrometer. It is used to measure the outer size between two plane surfaces, or the diameter of a cylindrical pin surface.

The principle of the algorithm is explained using the caliper example as shown in Both measurements, using top end or bottom end, have similar processes. To illustrate the method, the top end example in Figure 5-14(a) will be used.

First, the caliper is put close to the final measurement position, which is shown in Figure 5-15(a). The distances between 𝑃𝐿 1 and part surface 𝑆 1 are evaluated by calculating the distances between sampling points and part surface.

Next, surface PL1 is associated to the part surface. The maximum distance between them is minimized without part interpenetration. After association, the position of the main part of the caliper is fixed, the slider can move only along the guide way.

The distances between 𝑃𝐿 2 and part surface 𝑆 2 are calculated avoiding interpenetration. Then the slider is moved to minimize the distance 𝑑 𝑚 between 𝑃𝐿 1 and 𝑃𝐿 2 . The magnitude of 𝑑 𝑚 is the depth. The measurement of depth is expressed using GeoSpelling in Depth measurement using bottom end surfaces uses a similar process. However, since the size of the measurement surface is small, the caliper's initial position should be chosen carefully.

The linearization and solving of the problem use methods that are the same as the measurement cases introduced earlier.

Conclusion

The evaluation of geometrical defects is critical for quality control and process improvement.

Besides the precision limitation of measurement instruments, the method used to express specification also contains ambiguities.

In this chapter, the concept of method uncertainty is introduced first and its influences on product lifecycle explained. To solve this problem, the GeoSpelling-based specification expression was studied, and cases given as examples.

To carry out the GeoSpelling operations, the SDT-based deviation evaluation method was reviewed. The method was developed to calculate the invariance of a feature group automatically.

In addition, a virtual laboratory was developed for e-learning. It takes manufacturing uncertainties into account by introducing the skin model shapes. These defects can be visualized and evaluated. The SDT-based algorithm was modified for various pieces of measurement equipment to be able to conduct different measurement tasks.

In future work, measurement equipment like the Coordinate Measurement Machine (CMM)

and optical devices will be introduced into the VL. The corresponding measurement algorithms will be developed and integrated into the function package.

Conclusion

The management of geometrical deviation requires a consistent and complete model, avoiding ambiguities and covering the whole product lifecycle. To handle these problems, the concept of GeoSpelling and skin model were proposed, and they were retained as the basis of ISO standards.

GeoSpelling is a language based on limited operations. It was designed with coherent semantics in order to avoid specification related uncertainties (such as correlation uncertainty, specification uncertainty and method uncertainty). The design of GeoSpelling also considers the readability and learning cost for engineers. Moreover, it is a language which meets the requirements of digitalization and automation during application.

The manufactured part contains various types of deviation, but in current computer-aided tool systems, only the nominal model is integrated. To consider these deviations and uncertainties in product lifecycle management, the concept of the skin model was proposed. However, this model is continuous and infinite. To be integrated into computer-aided tools, its discrete instance, called skin model shape, is generated. The skin model shape can be used in tolerance analysis, metrology process improvement, etc. Additionally, visualization of the skin model shape helps engineers to understand the manufacturing defects and the semantics of GeoSpelling.

Contributions

Based on the framework of GeoSpelling and the skin model, this thesis concentrates on several implementation topics. The main contributions of the work are listed below:  A survey of deviation simulation methods for discrete skin model shapes has been conducted.

The simulation methods have been classified into three categories: random noise methods, mesh morphing methods, and mode-based methods. These methods have been compared, giving consideration to their characteristics, such as multi-scale, surface complexity, measurement data integration, parametric control and calculation complexity. More importantly, analysis has been provided from the point of view of engineering applications, which helps engineers or researchers to choose the most suitable solution.

 Solutions have been proposed to handle self-intersection or distortion on surface edges.

The earlier studies deviation simulation methods aim to simulate manufacturing defects on single surfaces. A complete skin model shape is required for product management, or for visualization purposes. However, due to manufacturing defects on the surfaces, once they are combined together, there can be problems on surface edges (such as self-intersection or distortion). To solve this problem, two solutions have been proposed. The first adds deviations to the nominal model directly, then the mesh is smoothed along local tangential directions. The second considers deviations on mesh as displacement boundary constraints, and solves the adjustment problem globally using FEA.

 A rigid body assembly simulation method using skin model shapes has been proposed.

One application for skin model shapes is assembly simulation and tolerance analysis.

The difference between nominal model, ideal surface-based deviation model and skin model shape have been explained. For the skin model shape, the influences of load and displacement boundary conditions cannot be ignored. Thus an assembly simulation method considering assembly loads has been developed. The method calculates the balanced situation under external loads using LCP and SDT based method.

 Verification has been conducted based on skin model shape.

The influences of specification uncertainty during the design, manufacturing and measurement stages have been introduced. To avoid specification-related uncertainties, the GeoSpelling-based approach has been studied. The semantics of specifications have been analyzed, and expressed using GeoSpelling. The operations used by GeoSpelling (such as partition, extraction and association) have been performed using SDT-based algorithms. Case studies have been conducted using the generated skin model shape.

 A library of algorithms for measurement has been developed for integration into an online virtual metrology laboratory. The laboratory will be used for teaching and self-learning. Models inside the laboratory are skin model shapes, so that manufacturing defects and measurement uncertainties can be introduced. The defects on skin model shape could be amplified to help users to understand the influence of form defects. The measurement uses SDT based algorithms.

Future Work

In this thesis, work has been done on the application of skin model shapes and GeoSpelling.

Efforts have been made to integrate skin model shapes into CAE systems. However, this work is not sufficient, and there are still many promising research topics to be conducted in the future.

Following the work in this thesis, applying skin model shapes to CAT systems, several critical studies are needed:  Simulate defects close to deviation in real manufacturing.

To reflect the deviations on the actual part, the manufacturing defects on skin model shapes should be close to reality. This ensures that the subsequent application of skin model shapes, like assembly simulation, will provide reliable results.

The generation of skin model shapes must be based on manufacturing process simulation and deviation measurement. The knowledge acquired from simulation and measurement enrich the generation process defining the type of geometrical defect and defect value to be considered.

 Simulate defects respecting the specifications.

Currently, the simulation of skin model shapes is conducted using random control parameters, or measurement data. If tolerances are defined on the CAD model, the skin model shape may not respect these tolerances. However, in assembly simulation and tolerance analysis, the models are assumed to satisfy corresponding specifications. The differences between design specification and skin model shape cause uncertainties, thus the simulation result is not reliable.

In earlier work, the control of deviation value on simple surfaces was studied.

However, the method is not general enough in consideration of surface complexity and specification complexity, and a more generic and reliable method is required. The problem could be solved in two ways. One method could be to analyze the specifications first, then generate control parameters that guarantee that the specifications have been fulfilled. Or a skin model shape could be generated which violate the specification. Then, based on the specification, the deviation value is scaled.

 Improve assembly simulation methods.

In this thesis, the importance of form defects, load boundary conditions and displacement boundary conditions are emphasized. An assembly simulation method is proposed to include these factors in the analysis. With the proposed methods, the assembly problem could be solved efficiently. However, there is still room for improvement.

The most time-consuming part of the method is the search for potential contact point pairs. A more efficient algorithm could be found in the field of Computational Geometric (CG). For the LCC-based method, the time required for solving quadratic optimization problems increases with the number of vertices. When it can be seen that only closest potential contact point pairs may be in contact in the next iteration, a point selection strategy could be introduced to reduce the number of point pairs.  Examine and improve GeoSpelling for metrology.

In current studies, GeoSpelling was applied to simple and common cases to illustrate its use. The most frequently used operations, such as partition and association, are studied. The formal language for operations, like extraction, collection, filtration and construction, still need to be studied. Meanwhile, the ISO standards for geometrical specification are still improving, and there is a need to express all these new added specifications using the GeoSpelling approach. In future, statistical characteristics and non-rigid workpieces should also be considered.

Apart from these specific problems, to implement skin model shapes and GeoSpelling that cover the whole product lifecycle, several problems should be considered for long-term development:

 The multi-scale representation of the skin model shape.

With the development of manufacturing technologies, the precision requirement on products is increasing gradually. Meanwhile, for the same part, the specification can be defined at various precision levels. This leads to the consideration of multi-scale solutions. Additionally, the storage and management of multi-scale skin model shape data throughout the product lifecycle may also be considered.

 Consideration of internal defects for skin model shape.

The skin model shape is generally considered as the external surface of the part, however there could be internal defects, like porosity for a cast part. With the application of additive manufacturing and topological optimization, parts with complex internal structure are widely generated. The defects on these internal structures influence the mechanical characteristics of the part. Therefore they should also be taken into consideration.

 Complex mechanical behavior in simulation.

In current research, the physical behavior of mechanisms is simplified for different purposes. To handle complex situations and generate more realistic simulation results, factors like friction force, local and global deformation etc. should also be considered. To consider these factors, simulation methods are not limited to FEA.

Related methods (such as boundary element method, point based simulation method, etc.) could also be considered.

 Extending the application of GeoSpelling and covering the whole product lifecycle.

GeoSpelling was initially designed for metrology to reduce specification uncertainties. However, similar uncertainties exist in various fields and stages in the product lifecycle. To improve the PLM, the concept of GeoSpelling could be extended to different applications, such as assembly process description. Currently, the assembly process is defined on ideal parts, not on parts with form defects. Boundary conditions must be introduced into GeoSpelling. Thus, there are many more possibilities still to be explored.

Figure 2 - 1

 21 Figure 2-1 Process for the generation of the Skin Model Shapes. ....................................................................

Figure 2 - 2

 22 Figure 2-2 Classification of simulation methods. ..............................................................................................

Figure 2 - 3

 23 Figure 2-3 Discrete rectangular surface model. ................................................................................................

Figure 2 - 4

 24 Figure 2-4 Random Gaussian deviation in 1D and 3D. ....................................................................................

Figure 2 - 5

 25 Figure 2-5 Result of 1D Gaussian simulation. ..................................................................................................

Figure 2 - 6

 26 Figure 2-6 Control points with sphere and ellipsoid influence hull. ...............................................................

Figure 2 - 7

 27 Figure 2-7 Simulation result of random mesh morphing. ...............................................................................

Figure 2 - 8

 28 Figure 2-8 Form error simulation on cylinder surfaces. ..................................................................................

Figure 2 - 9

 29 Figure 2-9 Form error simulation by second order shape morphing on rectangular surfaces. ......................

Figure 2 -

 2 Figure 2-10 Examples of Zernike polynomials. ................................................................................................

Figure 2 -

 2 Figure 2-11 Modes (u,v) of DCT used for simulation of form error.................................................................

Figure 2 -

 2 Figure 2-12 Simulation result by DCT. .............................................................................................................

Figure 2 -

 2 Figure 2-13 Embedding a line mesh by its eigenvectors of Laplacian matrix..................................................

Figure 2 -

 2 Figure 2-14 Eigenvectors (modes) generated from graph Laplacian. .............................................................

Figure 2 -

 2 Figure 2-15 Simulation of the plane with manufacturing deviation by graph Laplacian. ...............................

Figure 2 -

 2 Figure 2-16 Vibration modes used for deviation simulation. ..........................................................................

Figure 2 -

 2 Figure 2-17 Simulation result by natural vibration modes. .............................................................................

Figure 2 -

 2 Figure 2-18 Deformation patterns generated by eigenvectors corresponding to the 5 largest eigenvalues.

Figure 2 -

 2 Figure 2-19 Random field simulation result for different correlation lengths. ...............................................

Figure 2 -

 2 Figure 2-20 Decomposition modes by PCA based method. ............................................................................

Figure 3 - 1

 31 Figure 3-1 Examples of connected features. ...................................................................................................

Figure 3 - 2

 32 Figure 3-2 Estimation of vertex normal. ..........................................................................................................

Figure 3 - 3

 33 Figure 3-3 Two connected features with and without form error. .................................................................

Figure 3 - 4

 34 Figure 3-4 3D view of the intersection and gap between features with form error. ......................................

Figure 3 - 5

 35 Figure 3-5 Principle of deviation addition on nominal model directly. ...........................................................

Figure 3 - 6

 36 Figure 3-6 Combination result by deviation addition on vertices of nominal model directly. .......................

Figure 3 - 7

 37 Figure 3-7 Influence on the precision of features. ..........................................................................................

Figure 3 - 8

 38 Figure 3-8 Over-stretch and intersection of mesh when deviation is larger than the mesh size. ..................

Figure 3 - 9

 39 Figure 3-9 Principle of the umbrella operator. ................................................................................................

Figure 3 -

 3 Figure 3-10 Projection of umbrella operator during mesh regularization process. ........................................

Figure 3 -

 3 Figure 3-11 Principle of Spring Analogy. ..........................................................................................................

Figure 3 -

 3 Figure 3-12 2D mesh regularization example. .................................................................................................

Figure 3 -

 3 Figure 3-13 Example of FEA-based mesh regularization method in 2D. .........................................................

Figure 3 -

 3 Figure 3-14 Model used for simulation and comparison. ...............................................................................

Figure 3 -

 3 Figure 3-15 (a) Direct deviation addition with colour, (b-d) direct deviation addition, Laplacian mesh regularization and FEA-based regularization. .......................................................................................

Figure 3 -

 3 Figure 3-16 Changes at tangent area e4. ........................................................................................................

Figure 3 -

 3 Figure 3-17 Body of milling tool (a) Nominal model (b) Model after segmentation.......................................

Figure 3 -

 3 Figure 3-18 Skin model shapes with simulated deviations. ............................................................................

Figure 4 - 1

 41 Figure 4-1 Three kinds of model. .....................................................................................................................

Figure 4 - 2

 42 Figure 4-2 Mating of two parts with two plane surfaces. ...............................................................................

Figure 4 - 3

 43 Figure 4-3 Influences of form defects and force. ............................................................................................

Figure 4 - 4

 44 Figure 4-4 Assembly of two parts. ...................................................................................................................

Figure 4 - 5

 45 Figure 4-5 Three assembly configurations for skin model shapes. .................................................................

Figure 4 - 6

 46 Figure 4-6 Pin and hole assembly with gaps. ..................................................................................................

Figure 4 - 7

 47 Figure 4-7 Assembly of 3 nominal models. .....................................................................................................

Figure 4 - 8

 48 Figure 4-8 Three configurations for ideal surface-based deviation model assembly. ....................................

Figure 4 - 9

 49 Figure 4-9 Influences of assembly forces on configurations. ..........................................................................

Figure 4 -

 4 Figure 4-10 Pin, hole and washer assembly. ...................................................................................................

Figure 4 -

 4 Figure 4-11 Pin, hole and L-shaped part assembly. .........................................................................................

Figure 4 -

 4 Figure 4-12 Discrete models and relative displacements. ..............................................................................

Figure 4 -

 4 Figure 4-13 Two configurations between surfaces. ........................................................................................

Figure 4 -

 4 Figure 4-14 Different types of supports and reaction forces. .........................................................................

Figure 4 -

 4 Figure 4-15 Two square surfaces with manufacturing defects. ......................................................................

Figure 4 -

 4 Figure 4-16 Assembly forces change the contact positions. ...........................................................................

Figure 4 -

 4 Figure 4-17 Assembly of a cube. ......................................................................................................................

Figure 4 -

 4 Figure 4-18 Assembly of a cube, a pin and an L-shaped part. .........................................................................

Figure 4 -

 4 Figure 4-19 Variation of assembly configuration. ...........................................................................................

Figure 4 -

 4 Figure 4-20 Stability of the assembly configuration. .......................................................................................

Figure 4 -

 4 Figure 4-21 Balance of the discrete hole and pin under external force. .........................................................

Figure 4 -

 4 Figure 4-22 Section view of the hand saw model. ..........................................................................................

Figure 4 -

 4 Figure 4-23 Example of skin model shape for assembly simulation. ..............................................................

Figure 4 -

 4 Figure 4-24 Assembly simulation result and contact regions. ........................................................................

Figure 4 -

 4 Figure 4-25 Translation and rotation of saw body. .........................................................................................

Figure 5 - 1

 51 Figure 5-1 Method uncertainty. ......................................................................................................................

Figure 5 - 2

 52 Figure 5-2 Operations of GeoSpelling. ............................................................................................................

Figure 5 - 3

 53 Figure 5-3 Position tolerance example. ...........................................................................................................

Figure 5 - 4

 54 Figure 5-4 Generated skin model shape. .........................................................................................................

Figure 5 - 5

 55 Figure 5-5 Result of evaluation. .......................................................................................................................

Figure 5 - 6

 56 Figure 5-6 Example of median line extraction. ................................................................................................

Figure 5 - 7

 57 Figure 5-7 Two feature group association approaches. ................................................................................

Figure 5 - 8

 58 Figure 5-8 Methods to associate feature group. ...........................................................................................

Figure 5 - 9

 59 Figure 5-9 Interface of the virtual metrology laboratory. .............................................................................

Figure 5 -

 5 Figure 5-10 Measurement example using Vernier caliper. ...........................................................................

Figure 5 -

 5 Figure 5-11 Real time measurement using digital caliper. ............................................................................

Figure 5 -

 5 Figure 5-12 Principle of outer size measurement. ........................................................................................

Figure 5 -

 5 Figure 5-13 Principle of inner size measurement. .........................................................................................

Figure 5 -

 5 Figure 5-14 Measurement of depth. .............................................................................................................

Figure 5 -TablesTable 2 - 1

 521 Figure 5-15 Measurement of depth using top end surfaces. ........................................................................

  The development and implementation of GeoSpelling and skin model are global studies which cover various fields. Here, we concentrate on developing the following sub-objectives:  Generating skin model shapes containing manufacturing defects;  Taking form defect into consideration in tolerance analysis;  Applying skin model shapes in virtual metrology. The work is presented in the following sequence: In Chapter 2, the manufacturing defect simulation methods are studied. They are classified and compared based on their principles and characteristics. Simulation examples are given to show their applications. These methods are also compared based on various criteria.



  Deviation combination. When the geometric deviations are simulated, they are combined and added to the original nominal model.

Figure 2 - 1

 21 Figure 2-1 Process for the generation of the Skin Model Shapes: (a) Nominal model of the part, (b) Discretization and segmentation of the nominal model, (c) Simulate manufacturing deviations on each feature, (d) Combine all deviation data and nominal model into a complete skin model shape.

Figure 2 -

 2 Figure 2-1 shows the simulation result at each step. In Figure 2-1(b), the nominal model is discretized and segmented. Then for each feature, form defects are generated as shown in Figure 2-1(c). After form defect generation, the deviation values are combined with the nominal model to generate a complete skin model shape in Figure 2-1(d).

Figure 2 - 2

 22 Figure 2-2 Classification of simulation methods.

Figure 2 - 3

 23 Figure 2-3 Discrete rectangular surface model. Before introducing detailed simulation methods, we provide a general equation for the simulation methods. The skin model shape 𝑆𝐾 containing manufacturing defects is the vector of vertex coordinates, which could be reduced to deviations along the normal direction for some methods. Its corresponding nominal model (or mean model for statistical based method) is 𝑆𝐾 ̅̅̅̅ .

Figure 2 - 4

 24 Figure 2-4 Random Gaussian deviation in 1D and 3D.

Figure 2 - 5

 25 Figure 2-5 Result of 1D Gaussian simulation.

Figure 2 - 6

 26 Figure 2-6 Control points with sphere and ellipsoid influence hull.

Figure 2 - 7

 27 Figure 2-7 Simulation result of random mesh morphing.

Figure 2 - 8

 28 Figure 2-8 Form error simulation on cylinder surfaces.

Figure 2 - 9

 29 Figure 2-9 Form error simulation by second order shape morphing on rectangular surfaces.

Figure 2 -

 2 Figure 2-10 shows the shapes of the first 10 standard Zernike polynomials.

Figure 2 -

 2 Figure 2-10 Examples of Zernike polynomials [57].

The

  Figure 2-11, and one example of a generated surface with defect is shown in Figure 2-12.

Figure 2 -

 2 Figure 2-11 Modes (u,v) of DCT used for simulation of form error.

Figure 2 -

 2 Figure 2-12 Simulation result by DCT.

meshes in Figure 2 -

 2 [START_REF] Ballu | A new Design Method based on Functions and Tolerance Specifications for Product Modelling[END_REF] are embedded such that the 𝑦 coordinates are values of an eigenvector while the 𝑥 coordinates are the order number of the coordinates in the vector. Initially, their 𝑦 coordinates are set to zero, which is shown in Figure2-13(a). Eigenvectors in Figure2-13(b-d) correspond to the 2 nd , 4 th and 10 th eigenvalues respectively. As can be seen from the figure, as the eigenvalue increases, there are more oscillations.

Figure 2 -

 2 14. A surface generated by eigenvectors of a Laplacian matrix is shown in Figure 2-15. The coefficients are generated randomly.

Figure 2 -

 2 Figure 2-13 Embedding a line mesh by its eigenvectors of Laplacian matrix: (a) Initial position (b) 2 nd eigenvector (c) 4 th eigenvector (d) 10 th eigenvector.

Figure 2 -

 2 Figure 2-14 Eigenvectors (modes) generated from graph Laplacian.Simulation by the combination of weighted eigenvectors is a general way to simulate the manufacturing deviations. Unlike DCT or Zernike Polynomials, the calculation of modes is not limited by the shape of the nominal geometry.

Figure 2 -

 2 Figure 2-15 Simulation of the plane with manufacturing deviation by graph Laplacian.

  considered are numbers 1, 3, 7, 8, 9, 34 as shown in Figure 2-16. By specifying random coefficients following normal distribution, a surface with defects like those in Figure 2-17 can be generated.

Figure 2 -

 2 Figure 2-16 Vibration modes used for deviation simulation. As shown in Figure 2-17, modes in different frequencies have quite different patterns; this makes it possible to simulate various form errors from these modes. The simulation result has some obvious features of Mode 34. This is because we treated modes with different frequencies equally and weight factors were given randomly. Normally, these artifacts may not appear in a real application, as coefficients are generated from the real measured data or based on experience.

Figure 2 -

 2 Figure 2-17 Simulation result by natural vibration modes.

Figure 2 -

 2 18, some eigenvectors are visualized by their mesh color to show the change of values. The surfaces have certain deviation patterns like other mode based methods.

Figure 2 -

 2 Figure 2-18 Deformation patterns generated by eigenvectors corresponding to the 5 largest eigenvalues. For the example, we have considered the 20 largest eigenvectors, and coefficients following normal distribution are generated randomly. Correlation lengths 𝑙 𝜌 equal to 5, 10 and 20 are considered to show the differences.Figure 2-19 shows examples of generated surfaces.

Figure 2 -

 2 Figure 2-19 Random field simulation result for different correlation lengths.

Figure 2 -

 2 Figure 2-20 Decomposition modes by PCA based method[START_REF] Ballu | Comparison of mode decomposition methods tested on simulated surfaces[END_REF].Modes generated by PCA have similar properties to other spectral modes. As we do not have measurement data, an example from[START_REF] Ballu | Comparison of mode decomposition methods tested on simulated surfaces[END_REF] is used. The author decomposed simulated deviations

Figure 3 -

 3 2(b) is closer to reality. Thus in what follows in this paper, the vertex normal is estimated for single features after segmentation, as shown in Figure 3-2(b).

Figure 3 - 1

 31 Figure 3-1 Examples of connected features: (a) two connected mesh planes and the section plane in 3D, (b) 2D section view.

Figure 3 - 2

 32 Figure 3-2 Estimation of vertex normal, (a) before feature segmentation, (b) after feature segmentation.

Figure 3 - 3 shows

 33 Figure 3-3(a), where dashed lines represent a section view of two connected planes without form error, circles or squares in the dashed lines represent vertices on the faces, and the two faces are connected at the corner vertex. Assuming that the deviation of vertices outside the model is positive (+) and inside the model is negative (-), there are four types of configuration for the corner vertex, which can be seen in Figure 3-3(b). Solid lines, circles and squares represent the section view of features with form error. One can imagine intersections between neighbouring surfaces, but Figure 3-3(b) illustrates cases without intersection.

Figure 3 - 3

 33 Figure 3-3 Two connected features with and without form error.

Figure 3 - 4

 34 Figure 3-4 3D view of the intersection and gap between features with form error.

Figure 3 -

 3 [START_REF] Requicha | Toward a Theory of Geometric Tolerancing[END_REF] shows an example of deviation addition on an edge vertex.

  (a) Nominal features (b) Features with form errorsBy applying this combination method, Figure3-6 shows the combination results of Figure3-3(b) and the results look appropriate for these configurations.

Figure 3 - 5

 35 Figure 3-5 Principle of deviation addition on nominal model directly.

Figure 3 - 6

 36 Figure 3-6 Combination result by deviation addition on vertices of nominal model directly.

Figure 3 - 7

 37 Figure 3-7 Influence on the precision of features of (a) obtuse dihedral angle (b) acute dihedral angle.

Figure 3 -

 3 8(b) shows the combination result with amplified deviation. These overstretching and intersection problems could also happen in precision simulation, when the mesh size is smaller than or close to deviation value.

Figure 3 - 8

 38 Figure 3-8 Over-stretch and intersection of mesh when deviation is larger than the mesh size.

Figure 3 - 9

 39 Figure 3-9 Principle of the umbrella operator.

Figure 3 -

 3 Figure 3-11 Principle of Spring Analogy.

  -12 shows mesh regularization by spring analogy. Dashed lines indicate the nominal model while solid lines indicate the model with deviations.

Figure 3 -

 3 12(a) shows the model before regularization, and Figure 3-12(b) shows the model after regularization. As can be seen from Figure 3-12(b), vertices are adjusted only along their local tangent direction. The local mesh regularization method works in the same way for other situations, such as obtuse or acute dihedral angles.

Figure 3 -

 3 Figure 3-12 2D mesh regularization example (a) before regularization, (b) after regularization.

Figure 3 -

 3 Figure 3-13 Example of FEA-based mesh regularization method in 2D.

Figure 3 -

 3 Figure 3-14 Model used for simulation and comparison.

  -15(a) and (b). In Figure 3-15(a), the color denotes the simulation result from the random field method, while in Figure 3-15(b) the problems on the edges of the mesh are emphasized. It can be clearly seen that in Figure 3-15(b), the edges of the mesh have different degrees of stretch, compression and intersection.

Figure 3 -

 3 Figure 3-15(c) shows the result for the Laplacian mesh regularization. This is a local method and 200 iterations are conducted to smooth the mesh. With the proposed regularization strategy,

Figure 3 -

 3 Figure 3-15 (a) Direct deviation addition with colour, (b-d) direct deviation addition, Laplacian mesh regularization and FEA-based regularization.

Figure 3 -

 3 16(b) shows the result when the tangent features are considered together in simulation. Deviations of the tangent features are changing continuously. The decision to separate the two tangent features or not depends on the manufacturing process. Thus, to guarantee that the skin model shape represents realistic parts, it is better to keep the simulation process in accordance with the manufacturing process.

Figure 3 -

 3 Figure 3-16 Changes at tangent area e4 (a) Separated features, not continuous deviation, (b) Single feature, continuous deviation.

Figure 3 -

 3 17 (a). Its surfaces are segmented into single features, which are distinguished by different colours in Figure 3-17 (b). Modal based simulation method is used to simulate the manufacturing deviations. For each feature, four modes are chosen randomly and the maximum sum of the deviation could be 0.4.

Figure 3 -

 3 Figure 3-17 Body of milling tool (a) Nominal model (b) Model after segmentation.

Figure 3 -

 3 Figure 3-18 Skin model shapes with simulated deviations (a-c) Adding deviations on nominal model directly (d-e) Adding deviations on nominal model by FEA based method. The simulated deviations are amplified 15 times. In Figure 3-18 (a-c), the simulated manufacturing deviations are added on the nominal model directly.Figure 3-18 (a) aims at

Figure 3 -

 3 18 (a) aims at showing the global deviation shape in colour, while Figure 3-18 (b-c) emphases local mesh problems. As shown in Figure 3-18 (b), triangle meshes close to the corner and edge are heavily compressed. In the upper box of Figure 3-18 (c), due to the deviation values on vertices are larger than the mesh size of their neighbouring triangles, self-intersection happens and shapes like saw tooth are generated on the edge. In the lower box of Figure 3-18 (c), some triangles are

  compressed while others are stretched. Beside these four problems have been pointed out, there exist more problems if we check the model in detail.

Figure 3 -

 3 Figure 3-18 (d-e) shows the result by adding manufacturing deviations with FEA based global adjustment method. Same places are zoomed in to show the effectiveness of the method. It can be seen clearly that after adjustment, meshes close to the edge are more equally distributed, while the shapes of the geometry (original edges and the simulated manufacturing deviations) are conserved.

  Figure 4-1(a) and (b) shows the nominal model in both 3D and 2D (each edge in 2D corresponding to a plane surface in 3D). The nominal model represents the design intent and there are no manufacturing defects. Thus the surface (or edge) is a plane (or a straight line).

Figure 4 -

 4 1(c) shows one kind of deviation model. Here, the manufacturing defects are abstracted as the displacements of ideal planes (or edges). In the work of Li and Roy[START_REF] Li | Relative positioning of toleranced polyhedral parts in an assembly[END_REF], different simulation configurations are studied using this kind of deviation model. In our work, this kind of deviation model will be called the ideal surface-based model.Concerning the expression of manufacturing defects, skin model shape simulation goes further than the ideal surface-based deviation model. Based on engineering requirements, different levels of manufacturing defects can be simulated.

Figure 4 -

 4 1(d) shows the skin model shape with form defects. These form defects are considered as the combination of various kind of deviations, including the translation and rotation of nominal features, as in the ideal surface-based deviation model.

Figure 4 - 1

 41 Figure 4-1 Three kinds of model: (a) 3D nominal model (b) 2D nominal model (c) ideal surfacebased deviation model (d) skin model shape.

Figure 4 -

 4 Figure 4-2 shows the assembly, and the red arrow indicates the assembly force. The nominal model should be perfectly mated, as can be seen in Figure 4-2(a).

Figure 4 - 2

 42 Figure 4-2 Mating of two parts with two plane surfaces: (a) nominal models (b) polyhedral deviation models (c) skin model shapes.

Figure 4 - 3

 43 Figure 4-3 Influences of form defects and force: (a) force applied to center point (b) force applied close to left side.

  understanding, a close view of the mating surfaces is shown in Figure4-3. The two lines in different colors represent two mating surfaces in 2D. In Figure4-3(a), the assembly force is applied close to the center of the surface, and one mating configuration is found. When moving the assembly force to the left side of the surface, the contact regions are also shifting. The shifting of contact regions is determined by the form defects and assembly force at the same time. This

Figure 4 - 4

 44 Figure 4-4 Assembly of two parts: (a) nominal models (b) assembly of ideal surface-based deviation models, point on the corner is in contact with the surface (c, d) two possible cases where the loads applied to the model change.

Figure 4 - 5 ,

 45 Figure 4-5, with the changing of position and direction of assembly loads, the assembly configurations also change.

Figure 4 - 5

 45 Figure 4-5 Three assembly configurations for skin model shapes.

Figure 4 -

 4 6(a) shows the section view of the pin and hole assembly using nominal models. For ideal surface-based deviation models in Figure4-6(b), the possible contact points can be found by searching for the boundary features, while for skin model shapes in Figure4-6(c), it will be very difficult to search for all potential contact points on the surfaces.

Figure 4 - 6

 46 Figure 4-6 Pin and hole assembly with gaps: (a) nominal models (b) ideal surface-based deviation models (c) skin model shapes.

Figure 4 - 7

 47 Figure 4-7 Assembly of 3 nominal models.For ideal surface-based deviation models, due to the gaps caused by manufacturing defects, there are more configurations when loads change. Figure4-8 shows three different configurations.

Figure 4 -

 4 8 shows three different configurations. When considering an assembly of multiple parts, the forces may be applied to one part only (Figure 4-8 (a)) or to several parts (Figure 4-8 (b) and (c)). Occasionally, some forces may be neglected, reducing the forces applied to several parts to forces on one part only, thus simplifying the problem. Depending on the position, orientation and magnitude of the forces, different configurations appear.

Figure 4 - 8

 48 Figure 4-8 Three configurations for ideal surface-based deviation model assembly, different results are generated under different assembly loads. The influences of the loads are explained with skin model shapes.Figure 4-9(a) shows the case

  Figure 4-9(a) shows the case when the dominant force is applied to only one part, and other forces are negligible. Based on this boundary condition, the parts are positioned as shown in Figure 4-9(a). Generally, more than one

  force must be considered. This case is explained by Figure4-9(b), where the gravity of parts and assembly forces are considered.

Figure 4 - 9 Figure 4 - 10 .

 49410 Figure 4-9 Influences of assembly forces on configurations: (a) force applied to only one part (b) force applied to multiple parts.The influence of forces is also illustrated by the pin, hole and washer assembly, as shown in Figure4-10. If no force is applied to the washer, its position floats around the pin without penetration. When force is applied to the washer, there is displacement and contact with the pin.

Figure 4 -

 4 Figure 4-10(b) shows one possible configuration.

Figure 4 -

 4 Figure 4-10 Pin, hole and washer assembly: (a) no force applied to the washer (b) force applied to the washer. The example in Figure 4-10 could be modified as shown in Figure 4-11. In this example, the washer is extended to become an L-shaped part. Let us consider a force applied to the lower part of the L-shaped part. When sufficient force is applied to the pin, as in Figure 4-11(a), there are contacts between the L-shaped part and the base part. If the force on the pin is reduced, the Lshaped part may rotate, as shown in Figure 4-11(b).

Figure 4 -

 4 Figure 4-11 Pin, hole and L-shaped part assembly: (a) sufficient force applied to pin (b) force applied to pin is reduced.

Figure 4 -

 4 Figure 4-12 Discrete models and relative displacements: (a) models before displacement (b) models after displacement.

Figure 4 -

 4 Figure 4-13 Two configurations between surfaces: (a) gap (b) contact.

Figure 4 -

 4 2(b)), or to model contact with adjacent parts (application example in section 4.4.5). Support contains constraints on part displacement, and reaction forces from the support: Constraints on displacement: the displacement of the part along/around a certain direction is constrained. In the case in Figure4-2(b), the translation of the part in a tangential contact direction is constrained.

Figure 4 - 14 .

 414 The contact position is defined by the user in light of the assembly situation. The support can provide a displacement boundary constraint either along a local normal direction, or along a local tangential direction, or a combination of the two.

Figure 4 -

 4 Figure 4-14 Different types of supports and reaction forces: (a) support and reaction force along local normal direction (b) support and reaction force along local tangential direction (c) combination of normal and tangential directions.

a)

  Choose the fixed part(s); b) Define pairs of mating surfaces; c) Define load boundary conditions; d) Formulate the problem, and solve the optimization problem; e) If the calculated 𝑺𝑫𝑻𝒔 is smaller than a predefined value, stop the simulation; if not, go to step d) and repeat the optimization. In the next section, detailed examples and comparisons are shown to validate the effectiveness and efficiency of the proposed method.

- 3 .

 3 The two square surfaces with form defects are shown in Figure4-15. The bottom surface is fixed and the assembly force is applied to the top surface. Vertices on the surfaces have normal vectors along the 𝑍 direction, which is also their reaction force direction. By changing the position of assembly forces, different configurations are observed. As can be seen from Figure4-16, in the first case the assembly force is applied at the center point (Figure4-16(a)), and in the second case the force is applied close to one corner (Figure4-16(b)). The magnitude of the force is equal to 20N. The contact points between two surfaces (red dots) also change according to the position of the assembly forces. This means that the orientation and position of the mating surfaces change.

1 .

 1 AF indicates the external assembly force applied to the surface, 𝑟𝑓 1 , 𝑟𝑓 2 and 𝑟𝑓 3 corresponding to the reaction forces on each contact point respectively. As can be seen from the Table, all reaction forces are positive, which means the solution is feasible and unique.

Figure 4 -

 4 Figure 4-15 Two square surfaces with manufacturing defects.

Figure 4 -

 4 Figure 4-16 Assembly forces change the contact positions: (a) force is applied to the center (b) force is applied close to one corner.

- 17 .

 17 Several forces are applied to the cube, their positions and directions are shown in Figure 4-17(a). The magnitude of these three forces is equal to 10N. With these assembly forces, the mating surfaces are in contact at the points highlighted in Figure 4-17(b).

Figure 4 -

 4 Figure 4-17 Assembly of a cube: (a) forces applied to the cube (b) contact points.

Figure 4 -

 4 Figure 4-18 Assembly of a cube, a pin and an L-shaped part.

Figure 4 - 4 . 4 . 4

 4444 Figure 4-19 Variation of assembly configuration: (a, b) assembly force is sufficient (c, d) assembly force is decreased.4.4.4 Stability of the Assembly ConfigurationIn former simulation examples, although the configurations of the simulation results are balanced, they may not be stable. A turbulence in the assembly boundary condition may totally change the simulation result. We therefore analyze the stability of the simulation result in the following.

Figure 4 -

 4 Figure 4-20 Stability of the assembly configuration (a) mating of two surfaces (b) assembly of a cube.

Figure 4 -

 4 21(a) shows the initial relative position between pin and hole. Due to the form error, the force applied to the center point of the hole does not pass through the center point of the pin. A moment is caused by this external force and the hole rotates. The balanced position is shown in Figure4-21(b). As can be seen from the figure, the balanced position is not only determined by the external force, but the discrete expression of the surface may also influence the configuration.

Figure 4 -

 4 Figure 4-21 Balance of the discrete hole and pin under external force: (a) initial configuration (b) balanced configuration.4.4.5 Assembly of Mechanical ProductUsing these three examples, the characters of the proposed method are illustrated. However, they are too simple and not close to real engineering cases. Thus in this section a mechanical

Figure 4 -

 4 Figure 4-22 Section view of the hand saw model: (1) body (2) eccentric wheel (3) rod (4) slider and saw (5) wood material. The behavior of the wood in relation to the hand saw is modeled using the support and the reaction force of the support on the body and the saw. The rotation of the saw body around 𝑋 and

Figure 4 -

 4 Figure 4-23 Example of skin model shape for assembly simulation.

First, a batch

  of skin model shapes (20 for each part, a total of 80) are simulated. Next, we select them at random and conduct an assembly simulation based on the boundary conditions described above. The motor torque is estimated as 𝑻 𝑴 = 1000 (𝑁 • 𝑚𝑚). Hand force and torque are 𝑭 𝑯 = 100 (𝑁) and 𝑻 𝑯 = 500 (𝑁 • 𝑚𝑚) respectively. Gravity applied to the saw body is 2 (𝑁), and gravity applied to other parts is 1 (𝑁). The cutting force 𝑭 𝑪 and supporting forces 𝑭 𝑺𝟏~𝟒 are unknown external forces. 200 assemblies are generated and one of them is shown in Figure 4-24. The red dots indicate the contact regions between mating surfaces. For each assembly, the body are computed, and the histograms of these displacements for the 200 assemblies are presented (Figure 4-25).

Figure 4 -

 4 Figure 4-24 Assembly simulation result and contact regions.

Figure 4 - 25

 425 Figure 4-25 Translation and rotation of saw body: (a) translation along "X" axis, (b) translation along "Y" axis, (c) rotation around "Z" axis.

Figure 5 - 1

 51 Figure 5-1 Method uncertainty[START_REF] Mathieu | GEOSPELLING: a common language for Specification and Verification to express Method Uncertainty[END_REF].The classification of uncertainties is presented in Figure5-1. The correlation uncertainty is the difference between design intent and actual specification on the document. Different designers may have a different understanding of the symbolic GD&T notes, and translate the same functional requirement to different specifications. This uncertainty at the initial design stage influences all the following processes. Meanwhile, the specification itself contains uncertainties.

  system. But the correlation uncertainty and method uncertainty are ignored.

  [START_REF] Ballu | Formal Language for GeoSpelling[END_REF], a more formal version of GeoSpelling is proposed. Basic instructions like logical expressions and loops are introduced with examples.GeoSpelling is based on operations which could be applied to both ideal and non-ideal features. These operations are: partition, extraction, filtration, association, collection and 5-2 shows examples of each operation. In the application of GeoSpelling, we follow the formal definition proposed in[START_REF] Ballu | Formal Language for GeoSpelling[END_REF].

Figure 5 - 3 .

 53 Figure 5-3. Only datum features, toleranced features and tolerance value can be found on the drawing. There is no additional information to support the work of metrology. The appropriate interpretation of the specification implies the reading of the GPS standards. All the information is,

Figure 5 - 3

 53 Figure 5-3 Position tolerance example.To associate datum features and evaluate the deviation value, different measurement strategies may generate different results. To avoid this method uncertainty, in the following the specifications are expressed with the GeoSpelling approach. The skin model shapes containing manufacturing deviations are used to visualize the form deviations, which can be seen in Figure

Figure 5 - 4

 54 Figure 5-4 Generated skin model shape: (a) the complete model, (b) only related surfaces.

Figure 5 - 5 o

 55 Figure 5-5 Result of evaluation.Because the extraction of a median line from a nominally cylindrical surface is frequently used, a function called "MedianLine" is used[START_REF] Ballu | Formal Language for GeoSpelling[END_REF]. The process is as below:

Figure 5 - 6

 56 Figure 5-6 Example of median line extraction.

Table 5 - 4

 54 Linearized constraints on degree of invariance. and the straight line are 𝑃𝑇 1 (𝐶

- 7 .

 7 In Figure 5-7(a), two point clouds represent the measurement data of two cylinders. The two ideal cylinders (rectangles in 2D) are initially positioned close to the point clouds. The axes of the cylinders (dashed lines) are not parallel, and the distance between them differs from the nominal model. During the association process, angle and distance constraints are considered, and the result is shown in Figure 5-7(c).

Figure 5 - 7

 57 Figure 5-7 Two feature group association approaches: (a) two features could have relative displacements, (b) two features positioned according to TED, and they have the same displacements, (c) association results are the same for both approaches.

Figure 5 -

 5 7(b), because the two cylinders are extracted from the nominal model, their axes are parallel, and the distance between them is equal to the nominal distance 𝑑.During the association process, the two cylinders will have the same displacements, thus the axes will maintain the distance and angle between them. Their final position and orientation are shown in Figure5-7(c).

Figure 5 - 8

 58 Figure 5-8 Methods to associate feature group.

  tasks with no programming knowledge, or they can perform complex functions by programming with the Lua language[135].In developing of the VL, DIGINEXT developed the platform, data base, server, 3D environment and basic operations, and I developed the measurement function package for instruments like calipers, micrometers, and dial indicators. The functions inside the package are based on the method introduced in this chapter, and they are customized for each instrument.

Figure 5 - 9

 59 Figure 5-9 Interface of the virtual metrology laboratory.

Figure 5 - 10 .

 510 The caliper and skin model shape are first placed in a suitable position, which is shown in the operation area of the figure. The measurement can then be conducted by pressing a key. In the case of Figure5-10, the caliper will close until it comes into collision with the part. The orientation of the caliper is adjusted automatically to be perpendicular to the surface. The measurement result is displayed in Region 4. This result is also recorded and all the results of the session can be exported as a CSV file. The following sections explain how the previous developments on metrology have been used in the VL to execute the measurement operation.

Figure 5 -

 5 Figure 5-10 Measurement example using Vernier caliper.

4 5. 6 . 2

 462 Figure 5-10. Next, specific functions are used to mate the part surface and the measuring instrument. The measurement is carried out automatically. In the rest of this section, the measurement functions for instruments are introduced.

Figure 5 -

 5 Figure 5-11 Real time measurement using digital caliper.

Figure 5 - 12 .

 512 The two opposite plane surfaces of the caliper used for measurement (𝑃𝐿 1 and 𝑃𝐿 2 ) are shown in the figure. The distance between these two measuring surfaces is 𝑑 𝑚 . Points are sampled from the measuring surfaces. A part is placed between these two measuring surfaces, where the distance between surfaces 𝑆 1 and 𝑆 2 is to be measured. For the 𝑖-th point, a ray which5.6.2.3 Measurement of DepthUsing the caliper, two methods are provided to measure depth. In Figure5-14(a), depth is measured using the top surfaces of the caliper (shoulder measurement), while in the bottom surfaces are used. The top end has large contact surfaces which are stable in measurement, and the bottom end can be used to measure depth of hole.

Figure 5 -

 5 Figure 5-14 Measurement of depth: (a) using top end, (b) using bottom end.

Figure 5 -

 5 Figure 5-15 Measurement of depth using top end surfaces: (a) initial position, (b) mating the plane surface PL1, (c) minimizing the distance between PL1 and PL2.

Table 5 -

 5 7 Depth measurement expressed by GeoSpelling. ..........................................................................

	XIV

  ). Ces spécifications symboliques rencontrent des limites et ne permettent pas toujours d'éviter les ambiguïtés et d'exprimer de nouveaux besoins fonctionnels et de fabrication. Pour Le but du « skin » modèle est de représenter la géométrie des pièces réelles. Dans les systèmes de CAO, le produit est conçu et visualisé dans les conditions nominales, ce qui satisfait aux contraintes de conception mais ne contient aucun défaut de forme. En revanche, le « skin » modèle intègre les défauts de fabrication. C'est donc un lien entre la partie physique réelle et la partie nominale virtuelle. Pour les besoins de calcul et de visualisation, la notion d'instance de « skin » modèle a été présentée. Elle prend la forme d'un maillage représentatif des défauts géométriques d'une pièce virtuelle. à base de bruit aléatoire, des écarts aléatoires normaux à la surface (et parfois selon des directions tangentielles) sont générés au niveau des sommets du maillage de la surface suivant une distribution aléatoire qui peut être par exemple gaussienne. Cette méthode permet de simuler des défauts de surface locaux. Le procédé de morphing consiste à déformer le maillage de la surface soit selon une surface quadratique soit selon une déformation aléatoire.Une déformée quadratique, peut simuler efficacement un défaut de forme dû aux défauts cinématiques des machines-outils. La méthode basée sur des déformations aléatoires a l'avantage de permettre de contrôler directement la forme et la valeur de l'écart, ce qui rend le processus de génération flexible. Dans le cas des méthodes modales, les surfaces sont générées par combinaison linéaire de modes de déviation. Selon la technique utilisée pour le calcul des modes, ces méthodes sont classées en méthodes trigonométrique et en méthodes spectrales. analogie avec les ressorts, sont utilisées pour régulariser les maillages au niveau des arêtes. Ces méthodes étant basées sur les propriétés locales, elles nécessitent un calcul itératif avec le choix toujours délicat de critères d'arrêt. Après étude, ces méthodes apportent une solution mais ne représentent pas le meilleur choix pour l'objectif visé.Une méthode globale basée sur les Eléments Finis (EF) est ensuite introduite pour générer des formes de « skin » modèle tout en évitant les inconvénients des méthodes locales. Dans la méthode basée sur les EF, les écarts de fabrication de chaque entité sont considérés comme des conditions aux limites de déplacement. Avec la méthode de pénalisation, la combinaison des écarts de fabrication et de l'ajustement du maillage sont effectué dans un seul calcul global, sans la simulation mettent en évidence deux défis pour les méthodes proposées. Le premier est relatif à la génération aléatoire des défauts de fabrication : les assemblages peuvent être impossibles si les défauts sont trop importants. Par conséquent, l'assemblabilité des pièces doit être vérifiée avant la simulation elle-même. De plus, initialement, la génération automatique des défauts de chaque pièce doit être maîtrisée et les formes de l'instance du « skin » modèle doivent vérifier les spécifications d'assemblage.Le second défi concerne le temps de calcul qui augmente rapidement avec le nombre de pièces, la taille des maillages et le nombre d'assemblages. En fait, environ 70% du temps de calcul est consacré à la recherche des points de contact potentiels les plus proches. Pour réduire le temps de recherche, l'algorithme et la structure de données pourraient être améliorés. Pour

	De plus, un laboratoire virtuel est en développement à l'université de Bordeaux pour
	l'apprentissage en ligne. Il prend en compte les incertitudes de fabrication en introduisant les
	instances de « skin » modèle. Ces défauts peuvent être visualisés et évalués. L'algorithme basé sur
	les torseurs de petit déplacement a été adapté pour divers équipements de mesure (pied à
	coulisse, micromètre, comparateur) afin de pouvoir effectuer différentes tâches de métrologie
	conventionnelle.
	généraliser et communiquer les spécifications géométriques sans ambiguïtés, Ballu et Mathieu ont
	proposé un langage (GeoSpelling) et un modèle conceptuel (« Skin » modèle). Ceux-ci ont été
	sélectionnés comme base de la norme IS0 17450-1.

GeoSpelling est un langage basé sur un nombre limité d'opérations. La conception de GeoSpelling a pris en compte la cohérence sémantique, la demande de numérisation et d'automatisation, la lisibilité pour les humains et les machines, et le coût d'apprentissage pour les ingénieurs.

À ce stade, la génération d'instances numériques du « skin » modèle et les activités connexes pour le tolérancement ne répondent pas encore aux besoins de l'industrie et un travail important de recherche doit être poursuivi. L'objectif de cette thèse est de compléter les recherches antérieures et de développer de nouvelles approches pour mieux répondre aux demandes.

Au chapitre 2, les méthodes de génération de défauts virtuels de fabrication sur des surfaces isolées sont étudiées et comparées.

Dans la plupart des systèmes commerciaux de simulation des tolérances, les écarts sont générés par la modification des caractéristiques de position et d'orientation des surfaces nominales, qui ne contiennent pas de détails sur la forme des surfaces fabriquées. Cette approximation ne satisfait plus la demande croissante de précision du produit et n'est pas suffisante pour couvrir tout le cycle de vie. Certains chercheurs ont commencé à introduire des pièces virtuelles dans la simulation, avec des écarts de fabrication détaillés, comme les défauts de forme. Ceux sont encore des prototypes de laboratoire.

Le but de ce chapitre est de présenter des méthodes de génération appropriées pour créer des instances de « skin » modèle, qui sont des représentations discrètes de la pièce réelle fabriquée. Nous avons classé les méthodes de génération actuelles en trois catégories, les méthodes de bruit aléatoire, les méthodes de morphing et les méthodes modales. Cette classification n'est pas basée sur leurs origines, mais sur leurs approches mathématiques et numériques. Par exemple, les modes vibratoires et les champs aléatoires sont des méthodes ayant différentes origines, mais les deux sont basées sur un calcul de vecteurs propres. L'objectif est d'explorer les différentes méthodes disponibles, d'en observer les similarités et les différences et de montrer le type de résultat généré par chacune. XXI Pour la méthode Les méthodes de génération existantes couvrent une large gamme de types de surfaces et peuvent satisfaire la plupart des applications. Parmi ces méthodes, il n'est pas possible de définir de « meilleure » méthode répondant à tous les besoins, nous avons donc défini des critères pour les comparer de différents points de vue. Ces critères portent sur la prise en charge multi-échelles (forme, rugosité), la complexité des surfaces (plan, cylindre, sphère, cône, surface quelconques…), l'intégration de données de mesure (à partir d'un lot de pièces fabriquées), l'ajustement paramétrique (possibilité de modifier la surface à partir de paramètres) et la complexité de calcul. Cette analyse et la comparaison l'accompagnant doivent aider les chercheurs et les ingénieurs à mieux comprendre ces méthodes avant leur application et à éviter de réinventer des méthodes similaires. Au chapitre 3, sont présentés des méthodes pour générer des instances de « skin » modèle complètes et cohérentes. Avec les défauts de forme simulés au Chapitre 2 au niveau des surfaces individuelles, l'objectif de ce chapitre est de présenter différentes approches pour générer des instances de « skin » modèles cohérentes au niveau global de la pièce, de sorte qu'elles puissent être utilisées en tolérancement, en fabrication et en métrologie. La volonté étant de générer les surfaces individuellement afin de pouvoir contrôler les défauts indépendamment sur chaque surface, des problèmes géométriques surviennent lors de la combinaison de ces surfaces pour former la surface complète de la pièce (auto-intersection, amplification des défauts ...). Ces limites sont présentées et analysées. Deux solutions sont proposées pour obtenir une surface cohérente : une méthode locale et une méthode globale. XXII Les méthodes locales, telles que le lissage laplacien des mailles ou l'ajustement du maillage basé sur l'itérations. Des comparaisons entre les méthodes locales et globales sont effectuées, et la méthode globale basée sur les EF génère plus rapidement un maillage plus régulier que les méthodes locales. Au chapitre 4, la simulation d'assemblage est effectuée avec les instances de « skin » modèle générées. Tout d'abord, la question de l'assemblage en considérant les défauts de forme est discutée. Les différences entre l'assemblage de modèles nominaux, de modèles à base de surfaces idéales (avec défauts d'orientation et de position) et des instances du « skin » modèle (avec défauts de forme) sont mises en évidence. L'analyse montre que les conditions aux limites de l'assemblage (charge et déplacement) doivent être prises en compte lorsque des défauts de forme sont introduits. Ensuite, une méthode de simulation d'assemblage est proposée. La distance entre le point de contact potentiel et le plan est linéarisée à l'aide du torseur de petits déplacements. Les conditions de contact entre surfaces sont traduites sous la forme d'une condition de complémentarité linéaire entre les contacts et les forces de réaction. L'équilibre du système soumis aux charges d'assemblage est écrit sous la forme des équations d'équilibre statique d'un système mécanique. La méthode proposée considère les défauts de forme et les charges d'assemblage simultanément, tandis que l'aspect dynamique est ignoré pour simplifier le calcul. La simulation de différents exemples d'assemblage est effectuée dans la dernière section. De simple à complexe, l'efficacité de la méthode proposée est démontrée. Les régions de contact et les forces de réaction sont vérifiées pour garantir que le résultat de la simulation est réalisable. XXIII Les résultats de réduire le temps de résolution du calcul d'optimisation quadratique, le problème pourrait être réduit en ne considérant que les points de contact potentiels les plus proches. Ce chapitre propose donc une méthode générale pour mener une analyse de tolérance avec des instances de « skin » modèles (pièces virtuelles avec défauts de forme). La simulation prend en compte les conditions aux limites physiques en formulant le problème de contact sous la forme d'une condition de complémentarité linéaire et en utilisant les torseurs de petits déplacements pour linéariser les déplacements. Dans le chapitre 5, une méthode d'évaluation des défauts géométriques d'une instance de « skin » modèle isolé ou d'un assemblage est développée. La méthode est introduite dans le développement d'un laboratoire virtuel pour l'enseignement. L'évaluation des défauts géométriques est essentielle pour le contrôle qualité et l'amélioration des processus. Outre la limitation de précision des instruments de mesure, la méthode utilisée pour déterminer le défaut lié à une spécification introduit également des incertitudes. Dans ce chapitre, le concept d'incertitude de méthode est d'abord introduit et son influence sur la mesure est expliquée. Pour résoudre ce problème, l'expression des spécifications basée sur GeoSpelling a été étudiée, et des cas d'étude sont donnés en exemples. Pour mener à bien les opérations de GeoSpelling, la méthode d'évaluation des déviations basée sur les petits déplacements a été reprise. Dans les recherches antérieures, la méthode est appliquée sur des entités isolées. Pour associer un groupe d'éléments, deux approches sont proposées. La première approche considère le groupe d'éléments comme la combinaison de plusieurs entités isolées avec des contraintes de position et d'orientation entre chaque entité. La deuxième approche considère le groupe d'éléments comme une entité unique, ce qui signifie que XXIV la position relative entre les entités est définie directement par le modèle nominal de la CAO. Pour utiliser la seconde approche, la connaissance des degrés d'invariance du groupe d'éléments est requise. Une méthode basée sur les torseurs de petits déplacements et les sous espaces supplémentaires des déplacements a été développée pour calculer automatiquement l'invariance d'un groupe d'éléments.

Le dernier chapitre présente les apports de la thèse et les perspectives de recherche envisagées. Ainsi, il faudra intégrer les données issues de mesures lors de la génération de défauts de forme. Les spécifications doivent également être prises en compte lors de la génération de l'instance de « skin » modèle. La précision de la méthode de simulation d'assemblage proposée devrait être évaluée et l'efficacité du calcul numérique doit être améliorée. Pour les développements à long terme, l'expression des spécifications avec GeoSpelling doit être améliorée et les conditions d'assemblage doivent être intégrées. Enfin, les déformations sous charges doivent être prises en compte au niveau des pièces et au niveau des contacts lors de la simulation de l'assemblage. Les travaux présentés, ouvrent donc la voie à de nombreuses applications et futurs développements.

Table 2 -

 2 1.

Table 2 -

 2 1 Comparison of different simulation methods.

	Method	Classification Multi-scale	Surface Complexity	Measurement Data Integration	Parametric Control	Calculation Complexity
	1D Gaussian	Non-modes	★	★★★	★	★	★★★
	Second Order Shape Morphing	Non-modes	★	★★	★★	★★★	★★★
	Random Shape Morphing	Non-modes	★★★	★★★	★	★★★	★★
	Trigonometric Functions	Modes	★★★	★	★★★	★★	★★★
	Spectral Methods	Modes	★★★	★★★	★★★	★★	★★

  [START_REF] Charpentier | A scientific point of view of a simple industrial tolerancing process[END_REF] , it is projected to the closest segment on 𝑆 2 . The projection direction is the normal direction of this segment (𝒏 𝟐 ). A foot point 𝐴 2 is found, and the distance 𝑑 between 𝐴 1 and 𝐴 2 is the local distance between surfaces 𝑆 1 and 𝑆 2 . The displacements of the surfaces are expressed by SDT at their center point as [ 𝒓 𝟏 𝒕 𝟏,𝑪𝟏 ]

𝐶1

and [ 𝒓 𝟐 𝒕 𝟐,𝑪𝟐 ] 𝐶2

Table 4 -

 4 1 Assembly and reaction forces for cases in Figure 4-16.

		AF	rf1	rf2	rf3
	Case (a)	20 N	9.05 N	5.64 N	5.32 N
	Case (b)	20 N	3.63 N	12.73 N	3.64 N

  Simulation parameters are linked with physical boundary conditions by formulating the problem with SDT and LCC. However, although assembly forces are considered, the simulation is still handling rigid bodies. Thus in the next step, local and global part deformations could be introduced. With local deformation the interference fit can be simulated, which has a wide range of application in machinery. With global deformation it is possible to consider deformations caused by assembly loads.

  [START_REF] Charpentier | A scientific point of view of a simple industrial tolerancing process[END_REF] ) and 𝑆𝐿 2 (𝐶 2 , 𝒏 𝟐 ), after displacement [ 𝒓 𝟏 𝒕 𝟏,𝑪𝟏 ] and [ 𝒓 𝟐 𝒕 𝟐,𝑪𝟐 ], they will be 𝑃𝑇 1 (𝐶 1 ′ ) and 𝑆𝐿 2 (𝐶 2 ′ , 𝒏 𝟐 ′ ).a) When the ratio between displacement and distance is small, we calculate vector 𝒏 𝟐𝟏 which is parallel to the plane defined by 𝐶 1 , 𝐶 2 and 𝒏 𝟐 , and orthogonal to 𝒏 𝟐 :

	𝒏 𝟐𝟏 =	𝐶 2 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ -(𝐶 2 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ⋅𝒏 𝟐 )𝒏 𝟐 ‖𝐶 2 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ -(𝐶 2 𝐶 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ⋅𝒏 𝟐 )𝒏 𝟐 ‖

  𝒊,𝟏 , 𝒗 𝒊,𝟐 , ⋯ , 𝒗 𝒊,𝒏 ]

	𝑇
	𝑃 𝑖 = [𝒗 𝒊,𝟏 ′ , 𝒗 𝒊,𝟐 ′ , ⋯ , 𝒗 𝒊,𝒎 ′ ]	𝑇

Table 5

 5

-7.

Table 5 -

 5

	% Associate PL1 to S1
	D1min=Dmatmin(PL1,S1);
	D1max=Dmatmax(PL1,S1);
	PL1=Associate(PL1,D1min>=0,D1max);
	% Minimize distance between PL1 and PL2
	D2min=Dmatmin(PL2,S2);
	dm=Distance(PL1,PL2);
	A12=Angle(PL1,PL2);
	PL2=Associate(PL2,[D2min>=0,A12=0],dm);
	% Evaluation of the depth
	DEPTH=Eval(dm);

7 Depth measurement expressed by GeoSpelling. % PL1,PL2: measurement plane surfaces on caliper % S1,S2: non-ideal surfaces on the part

o Constraints: the distance between SL5 and CY2a is zero (SL5 passing through CY2a), and the angle between SL5 and CY2a is 90° (SL5 is orthogonal to CY2a).

o Objective: minimize the maximum distance between SL5 and the extracted median lines (EML3 and EML4).

 The maximum distance between SL5 and EML3, EML4 is evaluated. This maximum value should be smaller than half of the position tolerance value (0.1mm/2=0.05mm).

As explained above, there is plenty of information hidden behind the graphical GD&T expression. Metrologists will interpret the specification based on personal understanding, and metrology process planning is also influenced by available instruments.

To avoid ambiguities, and to be readable for both engineer and computer, these operations could be expressed by GeoSpelling in the form of a programming language. The formal expression of these operations is shown in Table 5-1. The evaluation result is shown in Figure 5-5.  The median line (ML) is identified by the collection of all the center points of CE. For the expression in Table 5-2, function "Solve" and loop instruction "foreach" are used.

Table 5-1 Formal expression of position.

Corresponding expression using GeoSpelling is shown in

SET_PL is considered as a set of an infinite number of planes because SCY surface is considered as continuous. The expression in Table 5-2 is a definitional expression. To be applicable to metrology, one must consider a finite number of planes. The solutions are sampled from the infinity solution set. Figure 56shows an example of median line extraction.

In the above sections, method uncertainty and GeoSpelling are introduced. Application examples using skin model shapes are also given to illustrate its usage. GeoSpelling defined the operations for deviation evaluation. By composition of the operations, multiple types of geometrical specifications can be considered, beyond the standardized ones.

Elementary Surfaces Type

Elementary Surfaces Example

Degrees of Invariance MGDE Elements

General surface Using the concept of SDT, the constraints on displacements can be written as linear equations.

Assume the displacement of a surface is expressed by [𝑺𝑫𝑻] = [𝒓 𝒕 𝑪 ] 𝐶 at its center point 𝐶. The constraints on degrees of invariance are shown in Table 5-4, where the helical surface is not included as more detailed parameters about the helix are required. The general surface is also not included as no constraint is needed.

Linearization of Characteristics

To take non-linear relations into consideration, in the work of Ballu [START_REF] Ballu | Identification de modèles géométriques composés pour la spécification et la mesure par coordonnées des caractéristiques fonctionnelles des pièces mécaniques[END_REF], the non-linear distance and angle constraints between features are linearized using SDT. In the following, two examples corresponding to the linearization of the orientation characteristic and the distance characteristic are given. For more detailed cases, please refer to related works [START_REF] Bourdet | The concept of the small displacement torsor in metrology[END_REF][START_REF] Ballu | Identification de modèles géométriques composés pour la spécification et la mesure par coordonnées des caractéristiques fonctionnelles des pièces mécaniques[END_REF][START_REF] Mathieu | Virtual Gauge with Internal Mobilities for the Verification of Functional Specifications[END_REF].

Linearization of Orientation Characteristics

The orientation of features can be described by the angle between them. This angle is calculated based on the vectors associated to the two features. Assume the two features have normalized vectors 𝒏 𝟏 and 𝒏 𝟐 associated to them. After rotation 𝒓 𝟏 and 𝒓 𝟐 , they become 𝒏 𝟏 ′ and 𝒏 𝟐 ′ . The angles between vectors before and after rotation are 𝛼 and 𝛼 ′ . Linearization depends on the configuration of the vectors, 3 cases are considered: a) When (𝛼 -𝛼 ′ ) is small compared with 𝛼 and its supplementary angle (𝜋 -𝛼), then:

Because the displacement is very small, based on Equation ( 41), 𝛼 ′ is linearized as:

b) When 𝛼 = 0, then:

c) When 𝛼 = 𝜋, then:

Linearization of Position Characteristics

The relative position between features can be expressed by distances between them. These distances are abstracted as distances between points, straight lines and planes. The linearization of distance between point and plane has been used in assembly simulation in an earlier chapter.

Here, we take the distance between point and straight line as an example.

Where 𝑑 ′ and 𝛼 ′ are distance and angle, 𝑐 1 and 𝑐 2 are coefficients, 𝐶 and 𝜃 are constants. To be able to solve this kind of constraint, the constraints between characteristics are analyzed.

Linear Equation Constraints for Distance Characteristics

Based on the linearization work in Section 5.4.3.2, it is easy to combine distance characteristics. The linear combination of distances is shown below:

Because each distance can be expressed by SDT in the form:

Where 𝒂 and 𝒃 are coefficient matrices, 𝑒 is a value dependent on feature initial position.

Then we can transform Equation ( 67) into a linear expression:

To find solutions that satisfy the equation constraints, the related features should be associated at the same time.

Linear Equation Constraints for Angle Characteristics

For angle characteristics, a general example is:

Based on Equation ( 57), the angle characteristic can be expressed in the form:

Where 𝒂 is coefficient matrix, 𝒓 is rotation vector of the SDT of features, and 𝑒 is a value dependent on initial relative orientation between features. Then Equation ( 70) is transformed into:

Then the solution that satisfies the constraint between orientation characteristics can be found.

Calculate Degree of Invariance for Feature Group

The degree of invariance and corresponding expressions by SDT have been established for a single feature. In mechanical design, to perform certain functions, a group of features (pattern of features) is also commonly used. Problems in measuring and tolerancing feature groups are complex due to the combination of different types of feature making the system over-constrained.

The torsor variables 𝑡 𝑖,𝑃𝐹𝐺 ′ are different from 𝑡 𝑖,𝑃𝐹𝐺 due to the translation of the coordinate basis. The simplified matrix 𝑃 𝐹𝐺 * represents the displacement space of the feature group, and can be decomposed into two parts:

The left side of the matrix, with the columns that have non-zero coefficients, represents the positional directions of the feature group. The right side of the matrix, which has all coefficients are zero, represents the invariance directions. Therefore, the null space is identified by the partitioning of the matrix. The directions of row vectors in 𝐼 𝐹𝐺 * (invariance matrix 𝐼 𝐹𝐺 after transformation) are identified from 𝑡 ′ . Once 𝐼 𝐹𝐺 * is calculated, an inverse transformation should be conducted to get the invariance matrix 𝐼 𝐹𝐺 in the original coordinate basis. Each row in 𝐼 𝐹𝐺 is a base vector in 𝐵 𝐹𝐺 , which spans 𝕀 𝐹𝐺 .

Using this calculation method, the invariance of the feature group can be identified automatically.

Development of Virtual Laboratory

The developed method was introduced into a Virtual Laboratory (VL) for e-learning. Using skin model shapes, the manufacturing defects of the part can be considered, and these defects are evaluated by SDT-based algorithm. In the following, the development objectives of the virtual laboratory are described. The modifications to the SDT-based algorithm to fit different types of measurement equipment are explained. is perpendicular to the surface is used to find the distance 𝑑 𝑖 between the measuring surface and the part. 𝑑 𝑖 is the material distance which is expressed using SDT in Equation [START_REF] Hasegawa | Calculation of the fractal dimensions of machined surface profiles[END_REF].

Project and Objectives

During assembly, the position and orientation of the caliper and slider are adjusted to minimize the distance 𝑑 𝑚 . The distance 𝑑 𝑖 is constrained to be positive to avoid interpenetration between the caliper and the part. The measurement result is the distance 𝑑 𝑚 . The measurement process can be expressed by the GeoSpelling approach, as shown in Table 5-5. Using the linearization method introduced earlier, the measurement process can be executed.

Solving the problem relies on linear programming algorithms, such as the simplex method and the interior point method.

Measurement of Inner Size

The function to measure inner size is designed for the caliper. It is used to measure the inner distance between two plane surfaces, or the diameter of a cylindrical hole surface. The principle of the algorithm is explained using the caliper example shown in Figure 5-13.

The definition of 𝑑 𝑖 and 𝑑 𝑚 is similar to the example in Figure 5

Title: Assembly simulation and evaluation based on generation of virtual workpieces with form defect

Abstract:

The geometry of a real manufactured part differs from the virtual workpieces designed in Computer Aided Design (CAD) systems. This difference is due to the accumulation of unavoidable manufacturing deviations. The objective of this work is to implement virtual workpieces with form defects (Skin Model Shape) in engineering applications to meet the industry's increasing demands in product geometry quality management. Various aspects are covered here, in particular form defect generation, assembly simulation and virtual metrology. Methods to generate form defects on simple surfaces are reviewed and classified. Due to form defects, the combination of simple surfaces to generate a whole part led to inconsistency on the edges. A global FEA-based method and a local mesh smoothing based method are used to overcome this issue.

To predict the deviation of functional characteristics, assembly simulation is conducted using skin model shapes. An approach is developed based on the Linear Complementarity Condition and the Small Displacement Torsor to take into account assembly boundary conditions, such as displacements and loads. Methods to evaluate deviation values on skin model shapes are also studied. Product specifications are expressed with GeoSpelling, and evaluated using the Small Displacement Torsor method. The developed methods are integrated into an online Virtual Laboratory for elearning.

The above-mentioned studies complement and extend the tolerance management methods based on GeoSpelling and skin models.