Nicolas Gobillot

Fabrice Guet

David Doose

Christophe Grand

Charles Lesire

Luca Santinelli

Measurement-based real-time analysis of robotic software architectures

I. INTRODUCTION

Works that develop or evaluate advanced functions for autonomous robots flourish. They propose advanced control techniques, environment reconstruction, navigation in complex, unknown, dynamic environments, interaction with human-beings, etc. In order to have all these techniques be used on actual robots in interaction with their environment, we must be able to guarantee their safe and correct behavior. This consideration has already been a major concern of roboticists for several years. Timing properties are crucial for the safety of the system as they could impact system execution and its functional performance. It is for instance proved that theoretically sound control techniques may become hazardous when facing timing disturbances [START_REF] Cervin | The jitter margin and its application in the design of real-time control systems[END_REF], [START_REF] Henriksson | On dynamic real-time scheduling of model predictive controllers[END_REF]. In this work we refer to a robotic system as a real-time system, where real-time systems are computing systems which must react within precise time constraints to events in the environment. Real-time systems are often called predictable or deterministic systems, since they demand for predictable and always reproducible behavior.

A. Real-time applications in robotics

Proving the real-time execution of applications has not been of major concern for roboticists, as they have mainly focused on algorithmic performance. Most of the works have considered that missing some deadlines would have no significant impact on the application. They have only considered making some measures of computation time of the functions [START_REF] Pizzoli | REMODE: Probabilistic, Monocular Dense Reconstruction in Real Time[END_REF], [START_REF] Sanfourche | Evo: A Realtime Embedded Stereo Odometry for MAV Applications[END_REF], [START_REF] Holzer | Adaptive Neighborhood Selection for Real-Time Surface Normal Estimation from Organized Point Cloud Data Using Integral Images[END_REF], [START_REF] Steinbrücker | Volumetric 3D Mapping in Real-Time on a CPU[END_REF] or response time over the network for distributed algorithms [START_REF] Wang | Hierarchical Auction-Based Mechanism for Real-Time Resource Retrieval in Cloud Mobile Robotic System[END_REF] for few runs, in order to evaluate timing properties of the systems. If the system is not strictly real-time, such measurements can suffice for having confidence that the algorithms can be embedded.

Authors are with ONERA -The French Aerospace Lab, Toulouse, France firstname.lastname@onera.fr

The lack of interest from roboticists for proving realtime constraints also comes from the applicability of realtime analyses. Indeed, robotics system are characterized by their modularity, the complexity of the algorithms and the diverse tasks executions models. On the other hand, realtime analysis methods need predictability and determinism for being applied. A way to obtain that is the use of the Worst-Case Execution Time (WCET), which is defined as the worst-possible execution time for a task. The WCET is a single value which is assumed as the execution time of the task in all its execution instances. The safety of real-time systems is assured once the timing constraints are respected.

B. Related works

Best practices towards real-time analysis of robotic systems have relied on model-based approaches. It is indeed necessary to have a model of the software running on the system in order to be able to check that real-time constraints will be ensured, depending on the scheduling policy of the system, and on the tasks properties (periods, deadlines, priorities). In [START_REF] Biggs | Modelling and Analysis of a Redundant Mobile Robot Architecture Using AADL[END_REF], the authors make an AADL model of their architecture, including hardware models, and analyse the latency of several components. However, they do not state how they get the worst case computation time of some functions. Moreover, their model is made a posteriori, then making the result questionable with respect to what is really implemented on the system. The SmartSoft framework [START_REF] Schlegel | Design Abstraction and Processes in Robotics: From Code-Driven to Model-Driven Engineering[END_REF] provides tools for describing components and how they are mapped to real-time tasks of the system. They eventually use the Cheddar tool [START_REF] Singhoff | Cheddar : a Flexible Real Time Scheduling Framework[END_REF] to simulate the scheduling of tasks in order to evaluate the schedulability of the architecture. However, they do no state how they get the WCET of tasks nor give any hint on how they could be computed.

C. Contribution

In this paper, we present an approach for evaluating the real-time execution of a robotic application that tends to be less pessimistic than classical real-time approaches and then more applicable to robotic systems.

This approach combines existing works into a global process for real-time analysis. The quality improvement of the analysis relies on two major points:

The approach is presented in Section II with a brief description of existing works that have been integrated. In the same section it is proposed a discussion on the way execution traces are obtained and applied. Section III then presents a case study in which this process is applied to the software architecture of a mobile robot performing an exploration mission.

II. REAL-TIME VALIDATION PROCESS

The real-time validation process presented in this paper is summarized on Fig. 1. First, it is necessary to make a model of the several components of the software architecture, and then to specify how they will be deployed on the system. The analysis needs to know the WCET of the elementary functions executed by the architecture which depend on the target platform. To obtain these WCET, the process is to generate some instrumented code and execute it on the actual target platform. From the traces obtained during this execution, it is possible to have a statistically sound evaluation of the WCET. From the models and these WCET, it is now possible to compute the Worst-Case Response Time (WCRT) of the several components in order to conclude on the architecture schedulability.

A. Component and Architecture Modeling

The modeling is done using a robotic, component-based Domain Specific Language, MAUVE [START_REF] Gobillot | A Modeling Framework for Software Architecture Specification and Validation[END_REF]. A MAUVE model of a software architecture is made of 4 nested levels:

• codels, which are elementary codes, i.e. functions implemented in C or C++ language and which implement the algorithmic part of the architecture; only codels signature is defined in the MAUVE DSL;

• components, composed of a shell describing their interface (input and output data ports, provided or required operations) and a core describing their behavior by a finite state machine; in each state, it is possible to call a codel or to access to ports (reading/writing) or operations;

• architectures, in which components are instantiated and connected one to the other;

• deployments, which indicate how the components of an architecture are mapped to real-time tasks, by defining their periods, priorities, deadlines and affinities (which CPU core they will run on). Listing 1 presents a snapshot of the Mauve model of a path planning component. Its shell is defined by a parameter (the path resolution), input ports (current pose of the robot, goal pose, current map), and an output port (the computed path). The core is defined by a state machine containing the Idle state (waiting for a new goal to arrive) and the Planning state. In that state, the component calls the astar codel. As the codels contain the algorithmic part of the components, the execution time of the tasks mainly rely on the execution time of the codels. It's also important to notice the model of the architecture is used in both the code generation and the schedulability analysis. Consequently, both the real execution and the analysis are based on the same formal specification.

B. Code generation and tracing

In order to do the schedulability analysis, we need to obtain timed execution traces of our components. From the MAUVE models of components and architecture, embedded executable code is generated using the Orocos [START_REF] Soetens | Realtime Hybrid Task-Based Control for Robots and Machine Tools[END_REF] middleware. Each instance of a component is a real-time task, where the behavior follows the state-machine description, and that calls manually defined codels. Architectures and deployments are plugged into Orocos deployment scripts which instantiate components, connect them together and define their activities.

Along with the component code, the code behavior is traced using the LTTng framework [START_REF] Desnoyers | The LTTng tracer: A Low Impact Performance and Behavior Monitoring for GNU/Linux[END_REF]. This framework has been chosen as it is known to have a low impact on the execution behavior 1 . These traces are generated within the components code so that we can log each important part of the architecture. The resulting trace provides:

• call time of codels, and the caller component name;

• return time of the same codels;

• start and end times of each component execution cycle;

• entry and exit times of each state of each component. From the timing data of these traces it is possible to extract the execution time of each codel execution by removing the interactions (preemptions or delays) of the more prioritary tasks. At the end of this step we obtain for each codel a sequence of execution times. It is then possible to get measures from several runs, where codels are executed, but not necessarily with the same architecture or in the same condition.

C. WCET estimation using the Extreme Value Theory

The static timing analysis [START_REF] Wilhelm | The worstcase execution-time problem -overview of methods and survey of tools[END_REF] requires a precise model of the hardware as well as a formal representation of the tasks to estimate WCETs. Unfortunately, some hardware used in robotics systems, such as commercial off-the-shelf platforms, are not well documented and cannot be formalized precisely. Moreover, some complex algorithm like SLAM or planning leads to pessimistic WCET representations.

Probabilistic timing analysis is emerging as alternative to the static timing analysis. As it is a probabilistic approach, it is the notion of probabilistic Worst-Case Execution Time (pWCET) that has to be considered. The pWCET is the worst distribution of execution times, where each value C i is associated with a probability p i ; p i is the probability to exceed C i . Probabilistic approaches, also called Measurement-Based Probabilistic Timing Analysis (MBPTA), only need measurements of execution time for inferring the pWCETs estimates. These approaches combine measurements and probabilistic analysis for WCET estimation [START_REF] Cucu-Grosjean | Measurement-Based Probabilistic Timing Analysis for Multi-path Programs[END_REF].

MBPTA approaches consist of two main steps: i) the collection of measurements of task execution time in real execution conditions as traces, and ii) the probabilistic analysis for inferring the pWCET estimation. The measurements are important for extracting observable features such as average behaviors and trends that can appear while executing tasks. In order to obtain the pWCET from the measurements, the probabilistic analysis is based on the Extreme Value Theory (EVT), allowing to infer the rare events (where the worstcases should be), very costly to measure. Doing this, the 1 http://lttng.org/features/#high-performance analysis provides safe pWCET estimates by mathematically inferring larger execution times than those measured. By safe estimates we mean that the pWCET is larger than or equal to the exact (and unknown) pWCET distribution, and to any possible empirical execution time distribution from the measurements.

Theoretically, the EVT is established for independent and identically distributed (i.i.d.) measurements [START_REF] Cucu-Grosjean | Measurement-Based Probabilistic Timing Analysis for Multi-path Programs[END_REF], but in practice the less constraining hypotheses of stationarity and extreme independent measure(ment)s are enough to apply the EVT [START_REF] Santinelli | On the Sustainability of the Extreme Value Theory for WCET Estimation[END_REF], [START_REF] Guet | On the reliability of the probabilistic worst-case execution time estimates[END_REF]. Thus, the EVT can be reliably applied to realistic cases (non time-randomized systems) with certain degrees of dependences between executions, such as the estimation of the WCET of codels in our component-based architecture. The safety of the pWCET estimates with EVT comes from the respect of the hypotheses, [START_REF] Santinelli | On the Sustainability of the Extreme Value Theory for WCET Estimation[END_REF], [START_REF] Guet | On the reliability of the probabilistic worst-case execution time estimates[END_REF], [START_REF] Cucu-Grosjean | Measurement-Based Probabilistic Timing Analysis for Multi-path Programs[END_REF].

The MBPTA is able to provide two main results: i) the WCET of each codels (WCET estimate) and ii) the WCET estimate reliability. While the WCET of the codels is used for schedulability analysis, the WCET estimate reliability is used to characterize the applicability of the MBPTA.

a) WCET estimates: A pWCET estimate is inferred from a trace of measurements providing tuples associating an execution time C to a probability p giving the confidence in not exceeding C. For instance, the confidence could be chosen to p = 10 -7 , allowing to have sufficient confidence (1 -10 -7 ) on the WCET model C.

b) WCET estimate reliability: The MBPTA tool we have developed is conceived as a logical workflow checking the applicability of the EVT with specific tests [START_REF] Guet | On the reliability of the probabilistic worst-case execution time estimates[END_REF]. For an input trace of execution times it is possible to provide a pWCET estimate and an associated reliability with regard to the EVT applicability hypotheses. In particular, the reliability qualifies the traces of measurements (e.g. lack of information, verification of certain statistical properties,) and defines the quality of the pWCET estimates in representing the worst-case execution behavior of the codels. The reliability accounts for the safe inference of the pWCET with the EVT indicating that resulting WCET estimates have a sense (if they do not diverge). Based on this reliability, we can decide whether more timing samples must be gathered and analysed or not.

D. Schedulability analysis

In order to check the schedulability of the software architecture, we need to compute the Worst Case Response Time (WCRT) of the architecture's components. The Response time correspond to the time between the triggering of the task and the end of its execution, including tasks preemptions and delays. A component is schedulable if its WCRT is lower than its deadline. Obviously, a complete system is schedulable if all the components are schedulable.

Classical schedulability analysis [START_REF] Klein | A practitioner's handbook for real-time analysis[END_REF] deals with monolithic task models, where each component of our architecture would have a unique WCET. Recent WCRT approaches [START_REF] Baruah | Dynamic-and static-priority scheduling of recurring realtime tasks[END_REF], [START_REF] Stigge | The Digraph real-time task model[END_REF], [START_REF] Gobillot | Periodic statemachine aware real-time analysis[END_REF] have used more detailed task models to improve the quality of the WCRT analysis. [START_REF] Gobillot | Periodic statemachine aware real-time analysis[END_REF] models the tasks as Periodic State Machines (PSM), where each PSM fires a transition at each execution cycle, and each transition has a specific WCET. We have then developed a transformation from the MAUVE models of components to PSM in order to have the WCRT analysis take into account the state machines of our components. This transformation is only used for the analysis, the executed code being generated from the initial state machine description.

III. EXPLORATION BY A MOBILE ROBOT

The real-time analysis process has been used on a software architecture of a real robot performing an indoor navigation mission. The aim of this experiment is to autonomously and safely join some waypoints in a partially known environment while updating a 2D map and avoiding (dynamic) obstacles. The robot used in these experiments is a two-wheeled Pioneer 3DX mobile robot from Adept Mobile Robots. It is equipped with an Hokuyo UTM-30LX laser range finder.

This experiment requires the following functionalities that need to be executed conjointly. First, a localization and mapping algorithm (SLAM) is used to improve the robot odometry and build a map based on occupancy grid model. Functions dealing with motion control and obstacle avoidance are also implemented. For each target waypoint, the navigation and the guidance functions compute the path to reach this target point and manage the robot motion orders required to follow this path. The last functions are associated to the system hardware, particularly the drivers of the laser range finder (Hokuyo) and the mobile robot (P3DX) which allows to control the robot motions and provides the odometry estimate based on encoders measurement.

A. Components and architecture models

Each system functionality is implemented in different components using the MAUVE DSL. Figure 2 gives an overview of the system architecture. The low-level functions (drivers, in orange) are in the following components: hokuyo is the laser driver that gives periodically a laser scanning; p3dx driver is the robot driver that takes linear and angular velocities as input and provides, as output, the robot pose computed from the encoders. The robot is controlled by a set of components (guidance, control and safety switch, in blue) in order to follow the path computed by the planning components (in green). A command from the operator can also arrive from the teleop component.

The navigation and SLAM functions are dispatched onto several components (in green): one is a wrapper of the GMapping algorithm [START_REF] Grisetti | Improved Techniques for Grid Mapping with Rao-Blackwellized Particle Filters[END_REF]; pose computes a pose from the robot odometry and the SLAM correction; navigation provides functions to first compute a path to join a target waypoint, based on the classical A* algorithm, and second send successive points of the path to the guidance component.

The safety switch has two functions: first it is a safety component that reduces or sets to zero the robot's velocity control if an obstacle is detected in front of it, and second can switch between the teleop command and the control command.

B. pWCET estimation

From the models of the above-mentioned architecture, we have generated the components' code and deployed these components to have some execution traces. The target platform we use on our Pioneer robot is a Intel i5 laptop with a 4-core CPU at 2.7GHz, running Ubuntu 14.04 configured to use the real-time SCHED FIFO Linux scheduler.

The pWCET estimation uses measurements of the execution time of elementary functions (codels). The estimation is independent of the complete software architecture running on the system, and consequently it is possible to make several runs with different architectures, and then combine the measurements for the pWCET estimation. The runs on our platform have then be split into:

• automatic control runs, in which we only deployed the control-related components; we manually sent some guidance goals to the guidance components, for several path lengths and on several environments (typically with/without obstacles) to have a good variability of conditions;

• SLAM runs, in which we only deployed the SLAMrelated components; we manually drove the robot while building a map, performing different loops over the environment, with several paths;

• navigation runs, in which we only deployed the navigation components with a static global map from the previous run, without moving the robot: the navigation component has been asked to compute paths with random objectives.

We incrementally performed all these runs by checking the pWCET estimate reliability criteria. Once the criteria are good enough, we stopped making more runs, and the computed pWCET estimates were considered reliable enough for the sequel of the schedulability analysis. An example of the pWCET estimate is given in Fig. 3 and4. Figure 3 shows the measured execution times for the avoid collision codel of component guidance. Figure 4 shows the result of the pWCET estimation. The inverse cumulative probability distribution of measures is represented by red circles. The pWCET estimates is given by the green curve. This curve gives estimation of the WCET of the codel for the desired probability threshold.

C. Schedulability analysis

The schedulability analysis then uses the components and architecture models, as well as an estimation of the WCET of all the codels called in the components. The deployment specification is given in Tab. I. We have set the period of almost all components to 100ms., while timeconsuming components, namely gmapping and navigation are scheduled at a period of 1s. All deadlines are equal to periods. The priorities have been defined so that the drivers are the components with the highest priority, in order to not miss data. Then the components that manage the safety and integrity of the platform have the highest priorities. Finally, the lowest priorities have been set to the less important components (with respect to robot safety). The affinity, i.e. the core on which the components' tasks will run, are set to reserve one core to the gmapping component, which is highly time consuming, while other components run on a specific core. In order to apply the schedulability analysis method, we first need to estimate WCET values. From the pWCET computations, we have to select the probability threshold p. Tables II and III present the several values of the WCRT of each component, depending on the probability threshold p. We remind the reader that p reprensents the confidence we can have that the actual WCET will not exceed the estimated value. We have made p vary from 10 -3 to 10 -9 , which are common values used in safety assessment of critical systems. We can first notice that the WCRT are logically increasing when we reduce p, as we want to be more and more confident in the WCET estimates. Moreover, for components running on core 1, we can notice that all components are schedulable down to confidence 10 -8 . For a confidence of 10 -9 , only the teleop and navigation components are not schedulable. A similar behavior can be observed for the unique component running on core 2: gmapping is schedulable for all chosen values of p, but for 10 -9 . A discussion on the impacts of the results on the architecture design is discussed hereafter.

D. From WCRT to architecture design

The WCRT analysis proves that the real-time constraints of all components will be respected at execution with a confidence at least of 10 -8 . In the case we need a higher confidence (i.e. p ≤ 10 -9 ), the WCRT analysis results will help us designing our architecture. For instance, regarding the gmapping component, two solutions are possible without changing the component code:

• increase the component period (e.g. to 2.5 seconds),

to have the component be schedulable with a high confidence, despite that the processing will be slower;

• make the component aperiodic, reacting to incoming data (scans and poses) when they arrive and when the component has time to process them; most of the time the processing will be faster, but a freeze of the processing is still possible. As a conclusion, the WCRT analysis, will first give some confidence in the fulfillment or real-time constraints for timecritical components, and will also provide some information useful for designing an embedded robotic architecture.

IV. CONCLUSION

In this paper, we have presented a methodology for analyzing the real-time properties of a component-based software architecture. This methodology relies on MAUVE models of the architecture from which instrumented embedded code is generated. We are then able to estimate the pWCET of components from actual measurements. From these pWCET and the models, we applied a WCRT analysis that provides schedulability results with a given confidence. We have also discussed on a real application of the presented process can first give some confidence on the critical parts of the architecture, and insights for the design of the overall software architecture.

From the experimental results we have conducted, several future works have appeared. First, we must have a better control on the real-time execution of components. The realtime skills provided by the existing real-time middlewares used in robotics are not sufficient. To do this, we plan to either directly generate real-time tasks for the operating system, or to develop a light RT middleware that would allow to better control the real-time execution. Second, the pWCET estimation can also provide an estimate of the time taken by the glue (i.e. the middleware and some calls like reading/writing within ports), from the measurements that contain both codel calls and components triggers. This glue estimation is available but has not been taken into account in the WCRT analysis yet. Finally, we plan to enhance the probabilistic WCRT computation by taking into account preemptions and delays probabilities, in order to have even more precise schedulability results.

V. ACKNOWLEDGEMETS

Fig. 1 :

 1 Fig.1: Real-time analysis process Each step of this process is further explained in this section.

Listing 1 : 8 9 15 var p a t h : P a t h 16 var n e w g o a l : b o o l = f a l s e 17 18 s t a t e m a c h i n e { 19 i n i t i a l s t a t e I d l e { 20 run = { 21 i

 18151618192021 Mauve code of a path planning component 1 s h e l l P a t h P l a n n i n g S h e l l { 2 p r o p e r t y r e s o l u t i o n : d o u b l e = 0 . 3 3 i n p u t p o r t map : O c c u p a n c y G r i d 4 i n p u t p o r t p o s e : P o s e S t a m p e d 5 i n p u t p o r t g o a l : P o s e S t a m p e d 6 o u t p u t p o r t p a t h : P a t h 7 } c o d e l a s t a r ( i n p : PoseStamped , i n g o a l : PoseStamped , i n map : O c c u p a n c y G r i d ) : P a t h 10 11 c o r e P a t h P l a n n i n g C o r e ( P a t h P l a n n i n g S h e l l ) { 12 var p o s e : P o s e S t a m p e d 13 var g o a l : P o s e S t a m p e d 14 var map : O c c u p a n c y G r i d f ( read ( g o a l , g o a l ) == NewData ) t h e n { 22 n e w g o a l = n e w p l a n ( g o a l , p o s e ) ; n s i t i o n i f ( n e w g o a l ) s e l e c t P l a n n i n g 26 } 27 s t a t e P l a n n i n g { 28 e n t r y { 29 read ( pose , p o s e ) ; read ( map , map ) ; 30 p a t h = a s t a r ( p o s e , g o a l , map ) ; 31 w r i t e ( p a t h , p a t h ) ; 32 } 33 t r a n s i t i o n s e l e c t I d l e 34 }}}

Fig. 2 :

 2 Fig. 2: Navigation component-based architecture. Ovals represent components. Edges represent port connections between components.

  Fig. 3: Execution time measures for codel avoid collision. Each circle represents a measure.

Fig. 4 :

 4 Fig. 4: pWCET estimate for codel avoid collision

TABLE I :

 I Deployment specification

	component	period (ms) deadline (ms) priority affinity
	p3dx driver hokuyo safety switch pose guidance control teleop navigation gmapping	100 100 100 100 100 100 100 1000 1000	100 100 100 100 100 100 100 1000 1000	10 9 8 7 6 5 2 2 3	1 1 1 1 1 1 4 1 2

TABLE II :

 II WCRT results for components running on core 1. Values are in milliseconds. Values in bold are exceeding the component deadline.

	component	1e-3 1e-4 1e-5 1e-6 1e-7 1e-8	1e-9
	p3dx driver hokuyo safety switch pose guidance control teleop navigation	8 28 31 32 35 38 45 764	10 32 35 36 39 42 49 817	11 33 36 37 40 43 50 841	13 36 41 42 45 48 55 876	17 40 46 47 50 53 60 906	29 52 61 62 65 68 75 984 1146 55 78 90 91 94 98 499

TABLE III :

 III WCRT results for components running on core 2. Values are in milliseconds. Values in bold are exceeding the component deadline.

	component 1e-3 1e-4 1e-5 1e-6 1e-7 1e-8	1e-9
	gmapping	703	711	715	722	764	994 2258

This work is partially supported by the CPSELabs project funded by European Community's Horizon 2020 Programme under grant agreement no 644400.


Jean Maihle

Jean Marc Linares
email: jean-marc.linares@univmed.fr

Jean Michel Sprauel

The statistical gauge in geometrical verification Part I. Field of probability of the presence of matter

Keywords:

a b s t r a c t Currently, most measurement and geometrical verification processes are based on the characterization of distances or angles between geometrical elements. In a previous paper, a method was developed for the verification of ISO specifications using the statistical information included in a set of points. This method is based on the development of measurement sequences and hence a sequence of geometrical constructions. The main problem in this verification process is to propagate uncertainties. With a correct expression of a virtual gauge, it is possible to avoid this propagation; indeed by using virtual gauge manipulation, initial best-fitted surfaces are never touched. The aim of this paper is to present a new approach using this statistical information to describe the measured part and develop a virtual gauge system for verification. This will be done in two parts. The first part will focus on the statistical modelling of the measured surface; the second part will describe a verification method based on a virtual, i.e. numerical, gauge.

Introduction

well as product quality. Moreover, in the past few years, companies have increasingly turned to outsourcing on an international level. Thus, mastering the assembly and lifespan of a mechanical system has become more and more complex. One of the most important parameters for this kind of problem is the inevitable geometric variation of a part when compared with its geometric specifications. These dispersions are due to the natural variability of the manufacturing process. At the present time, ISO geometric standards provide efficient concepts to specify parts. Thanks to current standards, variations can be limited thus respecting functional requirements. These specifications determine a geometric zone in which all real surfaces of the part are acceptable. Standards also provide methods which describe procedures for verifying these specifications. However, with the progress made in manufacturing processes and increasing demands for quality, tolerance intervals in geometric specifications are continually shrinking. Tolerance intervals are becoming near equals to the possibilities of the means of measurement. This situation has brought about an increase in the uncertainties of measurement. To cut back on verification uncertainties, two methods are possible: improving the quality of the In a global industrial context, companies have to cope with measurement process or the part verification process. stiff competition not only concerning costs, but also in rapidity as Studies done on improving measurement processes have enabled researchers to create efficient methods of calibration and to compensate for geometric errors stemming from the means of measurement. But using these methods does not completely diminish measurement uncertainties. They do not take into account instantaneous random deviations or a form defect which cannot be compensated for. New technologies, whether optical or mechanical, can noticeably increase the rapidity in acquiring a measurement. Increasing the number of points measured allows us to use more and more complex and more and more accurate surface featuring techniques. Added to the new generations of three-dimensional modellers, these techniques open up new possibilities such as:

• best-fitting complex shapes, • using statistical tools like Bayesian inference, and • accurately estimating tangent surfaces.

In addition to raw measurement uncertainties, the geometrical verification process amplifies the error bars of modelled surfaces. Nowadays, the typical method of geometric verification uses vectorial calculation. It is based on distance determinations: point/point, point/line, and point/plane. The disadvantage of this kind of calculation is the fact that uncertainties can be generated while creating geometrical constructions. With current calculation capabilities, it is possible to make calculations which are less penalising

ARTICLE IN PRESS

Nomenclature a E E E E E --→ E E E E E E E E v E E

Nbp

number of measured point prob(OM) probability of a point M to be inside the matter pa level of risk of a SMB P random vector representing mean estimated geometric parameters P i ith coordinate of the random vector P q(M) orthogonal projection of a point M r, s, t, , curvilinear coordinate r 0 , s 0 , y 0 , 0 , 0 fixed curvilinear coordinate representing a given point of a surface R radius S p,x , S p,y , S p,z parametric equation of the surface S, respectively, projected onto x, y, z E orientation vector of the random vector P This paper proposes a verification process using virtual calibration with a statistical criterion. The goal is to provide the metrologist with an indication of the risk he is taking when he claims that a part is geometrically standard. The novelty of this method is to include measurement uncertainties in the verification process. This indication must also make it possible to question the measurement strategy used, the sampling strategy and the number of points probed.

This paper is divided into two distinct parts. The first part will put forward a statistical modelling method of real surfaces. This model is presented in the form of a field of probability of the presence of matter (FPPM). It enables us to take any point in space and calculate the probability of its being within the boundaries of the object in question.

The second part deals with a method of geometrical verification of a part using virtual gauges. Virtual gauges will be modelled by a metric tensor and, thanks to the FPPM, will enable n parameter number us to have a geometrical check on the product using a statistical criterion.

State of the art

Geometrical verification using coordinate measuring machines (CMM) can be summed up in four different actions (Fig. 1):

• Measuring the real part.

• Determining the features of the real surfaces of the part.

• Interpreting geometrical specifications.

• Testing conformance.

A great variety of research work has been carried out on the measurement of geometrical quantities as well as on verification. These studies cover many fields which can be organized into four groups: best-fitting methods, uncertainty calculations, tolerance zone descriptions, and conformance testing. Fig. 2 proposes different channels of research for each group.

In the first part of this paper, we will present different studies done in the fields concerning how to determine the features of measured surfaces as well as fields concerning how to estimate uncertainty.

ARTICLE IN PRESS

Best-fitting

The data supplied by CMM is given in the form of point coordinates. To estimate the surface of the real part, best-fitting must be carried out. Best-fitting is the operation in which a given feature representing the real surface is associated with a set of measured points. This association is completed by optimizing a criterion like least squares, infinite norm or likelihood functions. Since the beginning of CMM use, research work has been done on numerous bestfi tting methods [START_REF] Bourdet | The concept of small displacement torsor in metrology[END_REF][START_REF] Cox | Least square adjustment in the presence of discrepant data[END_REF][START_REF] Takamasu | Basic concept of feature-based metrology[END_REF]. In fact, the minimized function is never linear, but to simplify the optimization problem, the objective function is often linearized. A method based on a likelihood function proposes a non-linear optimization association [START_REF] Sprauel | Contribution of nonlinear optimization to the determination of measurement uncertainties[END_REF]. Nevertheless, it is well known in measurement that the sampling of the real surface and hence best-fitting, imply that the estimation of the surface is provided with uncertainties.

Uncertainties

Research on uncertainty calculation for CMM can be organized into three classes:

• Experimental measurement of gauge blocks: Different artefacts, user-defined parts or geometrical gauges can be used to calibrate and determine uncertainty in CMM [START_REF] Antunes | Estimation of the precision and uncertainty of a calibration[END_REF][START_REF] Weckenmann | Monitoring coordinates measuring machine by user-defined calibrated parts[END_REF]. • Simulations of the measurement process: The Monte Carlo Method is commonly used for simulating complex measurement processes. It has thus been used as an alternative to treat the complex measurement of dimensions and to determine the probability distribution density of measured quantities numerically [START_REF] Schwenke | Assessment of uncertainties in dimensional metrology by Monte Carlo simulation. Proposal of a modular and visual software[END_REF][START_REF] Balsamo | Evaluation of CMM uncertainty though Monte Carlo simulations[END_REF][START_REF] Sommer | A systematic approach to the modelling and uncertainty analysis of measurement[END_REF]. • The statistical approach: Least squares and extreme fits are used in CMM and a variety of researchers have developed several numerical algorithms for three-dimensional data analysis. A statistical relationship method has been developed to determine the validity of the geometric model [START_REF] Kurfess | Statistical verification of conformance to geometric tolerance[END_REF]. The same research group created a method to estimate the uncertainty of measurement on envelope surfaces [START_REF] Choi | Uncertainty of extreme fit evaluation for three-dimensional measurement data analysis[END_REF]. Evaluating uncertainties in the vectorial tolerances using the least squares method has been proposed. The rotation and translation vectors of the surface were combined to form a homogeneous coordinate transformation matrix which, however, was approximated to the first order expansion [START_REF] Yau | Evaluation and uncertainty analysis of vectorial tolerances[END_REF]. A number of projects related to the study of a set of measured points have been carried out. Methods have been suggested to evaluate geometric errors obtained on parts when the manufacturing process is known [START_REF] Summerhays | Optimizing discrete point sample patterns and measurement data analysis on internal cylindrical surfaces with systematic form deviations[END_REF][START_REF] Henke | Methods for evaluation of systematic geometric deviations in machined parts and their relationships to process variables[END_REF]. Finally, research has been done on estimating the measurement uncertainty of small circular features measured by CMM [START_REF] Phillips | The estimation of mea ment uncertainty of small circular features measured by coordinate measu machines[END_REF]. 3.1. The "point, vector, parameter" model

Modelling surfaces using random vectors

E E
All geometric surfaces, which can be described in an equation, can be represented with a set of point/vectors and intrinsic parameters. This assembly of a point and one or more vectors allows locating and orienting the surface. The intrinsic parameters enable us to construct the surface. There are several solutions in modelling. They are based on sets of points and/or vectors representing the same surface (Fig. 3).

However, statistical descriptions of a surface are more efficient if the number of parameters is reduced to the bare minimum (principle of parsimony). Thus, the random vector modelling method is the one proposed in the table of "point/vector/scalar" modelling (Fig. 4).

In the rest of the document, all the scalars in this modelling (point coordinates, coordinates of the vector(s) and parameters intrinsic to the model) will make up a vector shown as P. To illustrate this approach, let's take the case of a cylinder (Fig. 5). The random vector of the cylinder is made up of a point, a vector and a radius.

The random vector P used here is derived from a minimal random vector which was used in best-fitting.

Statistical best-fitting: vector expectation and the covariance matrix

E

The real surface to be estimated is made up of an infinite number of points. Statistically, the infinite number of points of the surface produces a population to be characterized. The set of points evaluated can thus be considered as a statistical sampling of this population. Modelling by random vectors consists in estimating the first and second statistical moment of this vector, in other words, the mathematical expectation and the set of covariances in vector P.

E

The random vector

A random vector P of n dimensions belonging to < n is the vector [P 1 , P 2 , . . . , Pn] T . Its components. P i , are random variables. The random vector components are assumed to be continuous. The Studying the bibliography, we can assert that classic methods of geometrical verification of products do not allow us to claim without a doubt that a product complies. Current research does not directly take into account uncertainties in the verification of geometrical specifications. It is limited, first of all, to determining the magnitude of the measured feature; then, calculations of measurement uncertainties are made separately. To compensate for this problem, we are going to present surface modelling based on the notion of random vectors. In this way we will be able to determine the uncertainty of the complete surface. This method will be presented in this chapter and used in the rest of the document to obtain a statistical index of the compliance of a product to a system of geometrical specifications. probability density function f is defined as:

E E f : < n →]0, 1[ P → f (P) E ˝
The probability that P is found within a space, ˝, is given by the function f defined as follows:

Z Z Z E . . . f (a i )da i = prob(P ∈ ˝)
The standard properties defined in one-dimensional statistics remain valid:

• Probability density f has values in <+ Z Z R
• The value of the repartition function on < n is: . . . f (a i )

da i = < n 1 E RR R i E E i RR R
We assume that the components of random vector a are continuous with a probability density f which admits statistical moments. These moments are defined by the following integral: m k i = . . . a k i f (a i )du i with: k i ∈ N designating the moment of order k of the variable i of vector P.

If the random components P of the vector contain a mathematical expectation, the centered moments can be calculated:

k = . . . [a i -E(a i )] ki f (a i ) da i .
The concept of a covariance matrix can be very useful for representing the second order moments. E

The first two moments of the random vector

The preceding section defined the different moments of the random vector via the hypothesis of existence. Deriving from the basic concepts of multi-dimensional statistics, the expectation vector and the variance covariance matrix of a random vector can be defined.

When P, random vector of n dimensions and components [P 1 , P 2 , . . . , Pn] T , has first order moments, its expectation vector is defined as follows:

E ⎪ ⎪ ⎨ ⎪ ⎩ ⎧ E(P 1 )
E(P 2 ) E(P) .

⎪ .

E(Pn)

E Given that P, a random vector of n dimensions and components [P 1 , P 2 , . . . ., Pn] T which admits second order moments, and thus first order moments, its variance covariance matrix is defined as:

⎛ var(P 1 ) cov(P 1 , P 2 ) . . cov(P 1 , P n ) ⎞ ⎜ cov(P 2 , P 1 ) .

. . .

⎟ cov(p) = . . . . . ⎜ ⎟ ⎝ ⎠ E
The variance covariance matrix is a square matrix of n × n dimensions.

For example, in the case of a plane defined by a point C and a normal vector v (Fig. 6), the expectation vector and the variance covariance matrix are:

E

For the first statistical moment:

P = 1 C x C y C z v x v y v z T
For the second statistical moment: ⎥ ⎥ ⎦

E = ⎢ ⎢ ⎣ cov(P) ⎡ var(C x ) cov(

Best-fitting and estimating uncertainties

The first step in best-fitting is choosing a model to describe the measured surface. The model is an equation of the surface. It must enable us to make a best fit based on a statistical criterion and to express the residues of the best fit. The calculation is carried out by expressing the distance between probed points and the model to best fit. Models used in this paper are based on the description "point, vector(s), parameter(s)" described in Section 3.1 used as a basis for the random vector. These models are expressed in parametric form:

S E E p (u, v)
where u and v represent curvilinear coordinates of the surface. The distance between point M and the best-fit model is marked as d and given by the function d defined as:

d : < 3 → < d(M j ) → d j
The second step is the choice of the best-fitting criterion. This criterion depends on the hypothesis taken on the distribution function of residuals [START_REF] Sprauel | Uncertainties in CMM Measurements, Control of ISO Specifications[END_REF] written as the series f.

The criterion used is maximum likelihood which gives an independent generic approach to function f. The goal is to maximize the likelihood function defined as:

NbP Y ϕ = f [d(M j )] j=1
where NbP represents the number of data points.

ARTICLE IN PRESS

Function ϕ only presents one overall maximum for surfaces produced by our means of production. This property has been shown in previous work [START_REF] Aranda | A study of the likelihood function maximization against the signature shape[END_REF]. The condition to check is to solve the following system of equations:

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ϕ 1 2 = 0 . ⎧ ⎪ ∂P = 0 ⎨ ∂ϕ ∂P . ⎪ ∂ϕ . ∂Pn = 0
where n represents the number of parameters.

In the case of a Gaussian probability density function, the problem can be expressed in the following way:

1 2 X j=1 2 X ∂d j ⎡ ⎤ NbP ⎣ -1 d j ⎦ NbP NbP ϕ = √ 2 e ⇒ j=1 d j ∂P i = 0
If the n parameters, P i , of the geometric element associated to the digitized surface were perfectly defined, the standard deviation could also be estimated in the same way:

∂ϕ u t 1 X 2 v u NbP ∂ = 0 ⇒ = NbP d j j=1
However, such estimator would lead to a biased evaluation of , because a set of n parameters P i has already been derived from the acquired data. Therefore, the standard deviation of the measurement has to be computed with the following expression, also called residue of the least squares optimization:

u t 1 X 2 v u NbP = NbP -n d j j=1
This deviation, , can be propagated to deduce the covariance matrix of the estimated parameters P i , using classical differential expressions of the uncertainties. From the diagonal components of the covariance matrix, the error bars of P i are then easily calculated, since the statistical distribution of these random variables corresponds to a Fisher-Student law. Moreover, the covariance matrix is also useful to propagate the uncertainties of the fitted surfaces to any derived geometric element. It has to be pointed out that the propagation process of the standard deviation requires a precise definition of the derivatives ∂d j /∂P i . For that reason, a nonlinear algorithm, based on a Gauss-Newton type method, has been implemented to solve the optimization conditions. It does not require any coordinate transformation and is not sensitive to the selection of initial intrinsic parameters. It allows proceeding all the classical surfaces, i.e. lines, planes, circles, cylinders, cones and spheres.

Usually, the measurement uncertainties are assumed the same for all the digitized coordinates. Nevertheless, this hypothesis in not required to apply the maximum likelihood criterion. In fact, if the standard deviation, j , of each distance d j is well known, a weighted least squares criterion can be minimized to estimate the parameters of the fitted surface.

Field of probability of the presence of matter

Objective and description of the method

--→ --→
The next sections of the paper will define the probability that the real surface is, totally or partially included in the specified tolerance zone. For that purpose, the probability, written prob OM , that any point M of the 3D space is within the matter side of the real surface will be expressed. This model is called Field of probability of presence of the matter (FPPM). At the current step of the method, the random vector of the surface has already been estimated. It represents the mean value and the set of variance/covariances of the geometrical parameters of the vector/point/dimension obtained during the best-fitting (described in Section 3). The expression of the FPPM is done in three steps. First, the distance between any point of the 3D space and the mean surface has to be calculated. Its standard deviation has then to be derived from the variance/covariances of the random vector. Hence a parametric formulation of the mean surface must be written for these two operations. Finally the value of the FPPM is obtained by a simple integration of the Gaussian law. In this part, a formulation enabling us to define the function prob OM will be presented.

Parametric formulation of mean estimated surface

--→ E E P E
The first step in calculating probability prob OM is to find a parametric form for the mean estimated surface. This formulation will enable us to develop the estimated mean surface from the first statistical moment of random vector P. In some cases, a Cartesian formulation may appear more appropriate as in the case of a plane. However, such logic does not allow for a generic approach to the problem of verification.

The idea put forward here is to develop a parametric equation S E (r, s) representing the complete associated surface. It will be defined with the random vector P, where r and s are the curvilinear coordinates of the surface. Fig. 7 shows an example of a cylindrical associated surface.

The next section will show how to propagate the uncertainties of a random vector to a derived integral element.

Propagating uncertainties 4.3.1. The principle of propagation

The previous sections explained how to evaluate the uncertainties of a random vector by linearly propagating the variance of a one-dimensional random variable. The uncertainty propagation of a random vector of any dimensions toward another vector of any dimensions will be now presented.

E

Given any function g of < n in <, the variance of the value of g for a random vector P of n dimensions is expressed as:

E n n E E ∂P ∂P var[g(P)] =
XX ∂g(P) ∂g(P) cov(P i , P j ) j=1 i=1

i j

E

In the case of a function g of < n in < m , the covariance of the random vector g(P) is expressed as:

E E n m E E
∂P ∂P cov[g k (P), g l (P)] = XX ∂g k (P) ∂g l (P) cov(P i , P j ) i=1 j=1 i j

Using the Jacobian matrix of function g, the propagation of uncertainties can be written in a matrix form as follows:

E E E E E T E E cov[g(P)] = J g/P • cov(g) • J g/P With J E E = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ g m E E E ∂P ∂P P 1 1 1 E ∂P g m E E
-→ -→ P -→ -→ E P E E E T E E
to jth coordinate of P. In Section 4.2, the estimated mean surface of a real feature was represented by the parametric equations S E (r, s). The uncertainties of the random vector P will be propagated, through this equation, to any point of the surface. For a fixed set r and s, shown as u 0 and v 0 , the coordinates of point H are thus obtained within the mean surface: OH = S E (r 0 , s 0 ). The aim is to determine the threedimensional variance of this same point H. It will be calculated from the covariance matrix of the random vector OH. Using the method defined before. This propagation is achieved in the following way:

cov OH = cov[S E (r 0 , s 0 )] = J S/P • cov(P) • J S/P E E P E E E ∂S E E E E ∂S E P 2 E E P x E E P z ⎢ ⎥ With: ⎡ ⎤ ∂S Px Py ∂S Pz ⎢ ∂P 1 ∂P 1 ∂P 1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ J S = ⎢ ∂P x . . . . . . ⎥ ⎢ . . . . . . ... ⎥ ∂S ∂S ∂P k ... ∂P k ⎣ ⎦ OH E P E E E T T E E
In some cases, the parametric formulation of the mean surface requires using a local coordinate system. The make-up of the coordinate system must be formulated and the uncertainties of the initial random vector must be propagated to the reference point. The propagation formula is found to be:

cov -→ = cov[S E (r 0 , s 0 )] = J S/Rf • J Rf/P • cov(P) • J Rf/P • J S/Rf E E E
where Rf represents the local coordinate system, and J S/Rf represents the Jacobian matrix of the parametric equation in relation to the local coordinate system. The J Rf/P is the Jacobian matrix of the construction of the local coordinate system in relation to the random vector P. The next sections present an application of uncertainty propagation in the case of a plane and a cylinder. The Jacobian matrix J Rf/P of the construction of the coordinate system can be derived from the expressions of the coordinates of these vectors. Parametrizing the mean surface of the plane is achieved in Fig. 8. The parametric equation of the plane is:

E E 1 E E E e 1 x z v 0 E E v E ∧ e 1 E E E v E 4
E E E E S P (r, s) = r • e 1 + s • e 2
A fixed point H in the plane is expressed as:

-→ E E -→ E E OH = S P (r 0 , s 0 ) = OC + r 0 • e 1 + s 0 • e 2

OH /Rf

The Jacobian matrix used to propagate uncertainties toward point H is: " 1 0 0 r 0 0 0 s 0 0 0 0 0 0 # J-→ = 0 1 0 0 r 0 0 0 s 0 0 0 0 0 0 0 1 0 0 r 0 0 0 s 0 0 0 0

-→ E E E T T E E
In conclusion, the three-dimensional variance of point H within the estimated mean plane is:

cov OH = J S/Rf • J Rf/P • cov(P) • J Rf/P • J S/Rf E E 1 E E E

Propagating uncertainty in the case of a cylinder

The random vector P of a cylinder is P =

OCx OCy OCz vx vy vz R . To clarify the parametric equation, the local coordinate system Rf (C, e 1 , e 2 , e 3 , 1) can be constructed. The unit, shown as "1" in reference frame, is used as a basis for the parameter representing the radius of the cylinder. The components of the local coordinate system vectors are defined as follows:

e 1 x z 0 E v E ∧ e 1 v E E E v v v z -→ E 1 = p v 2 + v 2 . -v x , e 2 = ||E ∧ e 1 || , e 3 = || E || E
The Jacobian matrix of the construction of coordinate system Rf, named J Rf/P , can be derived from these expressions. Parametrizing the mean surface of the cylinder is achieved in Fig. 9.The parametric equation of the surface is obtained:

-→ E P E E OH = S E (, ) = v + Rn() With E E E n() = cos() • u + sin() • w
" 1 J-→ /Rf = 0 0 0 R • cos() 1 0 0 0 1 0 0 R • cos() 0 0 0 R • cos() R • sin() 0 0 0 R • sin() 0 0 0 R • sin() 0 0 u x • cos() + v x • sin() # 0 0 uy • cos() + vy • sin() 0 0 uz • cos() + vz • sin()
Finally, the three-dimensional variance of point H within the estimated mean cylinder is:

-

→ E E E T T E E cov OH = J S/Rf • J Rf/P • cov(P) • J Rf/P • J S/Rf

Formulating the field of probability of the presence of matter (FPPM)

--→ The goal of our statistical modelling is to clarify a function prob OM of < 3 in < which, at any point M in space, associates its probability of being within matter. This calculation is carried out in four steps (Fig. 10).

In the above section, the result of propagating uncertainties of the random vector toward point H within the estimated mean surface was the three-dimensional variance of this point. In other words, propagating uncertainties allows us to characterize the 3D distribution of probability representing the position, or the coordinates, of point H. As indicated in Fig. 11, point H = q(M) is the orthogonal projection on the mean surface of any point M in a 3D space. To define the FPPM or the probability that M is part of matter, only the dispersion of point H in the normal direction to the surface is required (Fig. 11).

E 2 OH n S T -→ -→ -→ n S OH -→ 2 2 2 --→ E
The dispersion of point H in the normal direction of the mean surface is obtained by projecting onto a unitary normal vector (shown as n S in Fig. 11). The projection is demonstrated by the relation: -→

• -→ = n S • cov OH • n S . To get the total local variance, adding the variance of the residues obtained in best-fitting (the statistical representation of the deviations due to the form defects and the random perturbations of the CMM) will suffice:

Total/ -→ = -→

.n S + Residues .

The standard deviation in the normal direction of the mean surface and the position of point H enables us to give a complete definition of the probability density function to integrate in order to calculate prob OM . To carry out this operation, the integration line has to be clarified. This line is defined by point H and direction n S (Fig. 12).

This representation makes it possible to visualize the bounds used in the integral: [-∞;d(MH)]. Our application will be restricted to surfaces of high quality. In this case, the probability density of the position of a point within an estimated mean surface has a Gaussian behaviour [START_REF] Sprauel | Uncertainties in CMM Measurements, Control of ISO Specifications[END_REF]. The following expression of FPPM is thus obtained:

OM MH E •n Z 2 n S 2 d x ( ) n S --→ 2 prob --→ = S q 1 e -1
Total/ -→ dx where d(x) represents the algebraic distance to the center of Gaussian distribution (point H) (Fig. 13).

Expressing iso-probable surfaces: statistical matter boundaries (SMB)

The FPPM permits to estimate the probability of any point of the 3D space to be inside the matter. In other word, this is a statistical representation of the limit of the matter. The FPPM do not express a finite boundary of the matter but a fuzzy one. However, it is possible to express a finite boundary for a fixed risk. This is done by the concept of statistical matter boundaries (SMB). With a mathematical point of view, this is an iso-probable surface of the FPPM. In a physical point of view, the SMB represents the estimation, for a fixed risk, of the envelope containing the matter. This section will show how to use the FPPM in order to calculate statistical matter boundaries accurately. First, the definition of SMB will be presented. Next, a practical method of calculation will be proposed, followed by examples of a cylinder and a plane.

Definition

--→

The statistical matter boundaries are defined as the isoprobable surfaces of the presence of matter. From a mathematical point of view, an SMB is delimited by the following space: M/p OM = p˛, p˛ ∈ ]0, 1[ . In other words, the SMB is the solution set of the following equation:

OM E •n Z 1 2 E E --→ h i 2 -1 d(x) p˛ = q e 2 Total/n dx -∞
2 Total/n --→ where OM is the unknown.

Practical calculation of an SMB

It is difficult to solve the equation of an iso-probable surface. However, it is possible to provide a solution in parametric form taking the uncertainty of a point on the estimated mean surface. Since the position of a point within an estimated mean surface has a Gaussian behaviour, the limit in variations of this surface can be calculated for a fixed risk. Maximum variations are calculated by a constant k times standard deviation. --→ The parametric equation of an SMB is written: SMB(r, t) = E

P q 2 n S -→ -→ S E (r, s) + k • Total/ -→ • ns (r, s)
, where ns (r, s) represents the vector normal to the mean surface with curvilinear coordinates r and s, and k the constant fixing risk.

Case studies

v E --→ -→ E E --→ -→ v E --→ E E

SMB of a cylinder

The least complex case study is the SMB of a cylinder. A local coordinate system does not have to be constructed. Consider a point C, a vector E, a scalar R and vector n() representing, respectively, a point within the axis of the cylinder, the direction of the cylinder axis, the radius, and a normal vector rotating around the cylinder axis. The cylinder can thus be defined in the following parametric equation (Fig. 11

): OM(, ) = OC + • v + R • n().
After propagation and projection along the normal to the cylinder, at given point, the SMB equation is as follows: SMB(, ) • n() which makes it possible to plot the SMB of the cylinder (Fig. 14).

= OC + • E + R • n() + k • 1 n() • cov(OC) • n() T --→ --→ E E E E
The result is a hyperboloid with an elliptic base. Physically, the hyperboloid shape can be explained by the uncertainties on orientation parameters of the estimated mean cylinder. This distortion will be found on every best-fit surface containing orientation parameters.

The SMB of a plane

The case of a plane is interesting because it is much easier to evaluate its SMB from a Cartesian equation. The random vector of a plane is in the following form: The surface obtained is a paraboloid (Fig. 15) where sections parallel to the mean plane are elliptical. As the cylinder case, the distortion obtained is depending on the uncertainty of the orientation parameters of the mean estimated plane.

Conclusion

In this paper, a statistical approach of CMM measurements has been proposed to characterize any real surface.

In this approach, a geometric element is no longer defined by a perfect feature but is described by a FPPM. This field, based on the concept of random vectors, defines the probability of any point to stay inside the matter. It integrates the random perturbations of the CMM, the uncontrolled factors of the sampling strategy and the form defects of the measured surface.

This concept can be applied to any standard surface (plane, sphere, cylinder, cone, . . .) and extended to models accounting for typical form errors thus increasing the significance of the best-fit method. This method has been tested with experiments on classical geometrical specifications [START_REF] Sprauel | Uncertainties in CMM Measurements, Control of ISO Specifications[END_REF]. The increase of the number of best-fit parameters to describe and fit the measured surfaces leads to a better characterization of the form defects and thus reduces considerably the uncertainties [START_REF] Eslter | Measurement as inference: fundamental ideas[END_REF].

The second part of our work will focus to the implementation of a virtual gauge method integrating the field of probability of the presence of matter. It will permit checking the geometrical conformance of a real part defining the risk of acceptance or refusal of the product. It will automatically include all errors or random perturbations and avoid any intermediate geometrical construction which amplifies the initial measurement uncertainties.

Fig. 1 .

 1 Fig. 1. Verification section.

Fig. 2 .

 2 Fig. 2. Fields of geometrical verification with a CMM.

Fig. 3 .

 3 Fig. 3. Point/vector descriptions of a plane.

Fig. 4 .

 4 Fig. 4. Point/vector/scalar modelling.

Fig. 5 .

 5 Fig. 5. Random vector of a cylinder.

Fig. 6 .

 6 Fig. 6. Random vector of a plane (point and orientation vector).

Fig. 7 .

 7 Fig. 7. Parametric formulation of a cylinder.

Fig. 8 .

 8 Fig. 8. Parametrizing a plane.

. 3 . 2 .

 32 Uncertainty propagation in the case of a planeThe vx vy vz . To clarify the parametric equation, a local coordinate system (C, e 1 , e 2 , e 3 ) can be constructed:

Fig. 9 .

 9 Fig. 9. Parametrizing a cylinder.

Fig. 11 .

 11 Fig. 11. Variance of point H.

Fig. 10 .Fig. 12 .

 1012 Fig. 10. The steps in calculating prob(M). Fig. 12. Integration principle.

Fig. 14 .

 14 Fig. 14. The SMB of a cylinder.

R 2 E+ 2 • 1 + 2

 2212 n() • cov(OC) • n() T + 2 • n() • cov(OC) • n() T i + res 2

E 1 PFig. 13 .

 113 Fig. 13. Integrating the probability density function.

Fig. 15 .

 15 Fig. 15. SMB of a plane.


F Gand

J Verrière

S Deck

Effect of upstream turbulence on single and dualstream jets. Assessment of Zonal Detached Eddy Simulation (ZDES)

This paper aims at investigating the effect of upstream turbulence (for instance coming from the engine's internal geometry) on the jet flow development. Simulations are performed using ZDES [1] combined with synthetic turbulence generation methods. Two cases are studied: an incompressible, singlestream jet and a compressible dual-stream one. It is shown that the upstream turbulence reduces the RANS-to-LES transition and improves the prediction of the location of the shock-cells and global jet flow development, which advocates a systematic consideration of realistic nozzle exit conditions in eddy-resolving simulations.

Introduction

Hybrid simulations, and more generally eddy-resolving ones, are now widely used in the aeronautic industry [START_REF] Deck | High-fidelity simulations of unsteady civil aircraft aerodynamics: stakes and perspectives. Application of zonal detached eddy simulation[END_REF]. The commonly used RANS/LES approaches (switching from RANS in attached areas to LES in detached ones) rely on the natural flow instabilities to trigger the development of the resolved turbulence in shear layers. These instabilities are rather weak in the case of jet flows, which leads to large delays in the RANS-to-LES transition and poor aerodynamic and aeroacoustic predictions. Some improvements of DES formulation have been proposed to mitigate this delay [START_REF] Deck | Recent Improvements of the Zonal Detached Eddy Simulation (ZDES) Formulation[END_REF] [START_REF] Shur | An Enhanced Version of DES with Rapid Transition from RANS to LES in Separated Flows[END_REF].

Another issue related to the use of eddy-resolving simulations for jet flows is the consideration of freestream turbulence in the jets core, which is inevitable in experiments and comes from the upstream parts of the engine in flight configurations. Few published studies have tackled this issue with different methods [START_REF] Brunet | Random Flow Generation Technique for Civil Aircraft Jet Simulations with the ZDES Approach[END_REF] [5] [START_REF] Tyacke | Large-Scale Multifidelity, Multiphysics, Hybrid Reynolds-Averaged Navier-Stokes/Large-Eddy Simulation of an Installed Aeroengine[END_REF]. It is noteworthy that, in all these cases, the introduction of upstream turbulence also improved the RANS-to-LES transition.

Therefore, the objective of this paper is to investigate the effect of adding freestream turbulence targeted at reproducing realistic inflow conditions on the turbulence development in mixing layers within the framework of ZDES simulations of jets. This is assessed on two configurations: a single stream, low speed jet and a more realistic dual-stream nozzle at high Reynolds number.

Zonal Detached Eddy Simulation (ZDES)

In this paper, hybrid RANS/LES simulations are performed using the Zonal Detached Eddy Simulation (ZDES) [START_REF] Deck | Recent Improvements of the Zonal Detached Eddy Simulation (ZDES) Formulation[END_REF]. Modes 1 and 2 of the ZDES are used. The ZDES mode 1, introduced to treat separations triggered by the geometry, relies on a user-defined zonal decomposition of the computational domain in RANS and DES areas. The hybrid length scale

I ZDES d ~ used for the Spalart-Allmaras (SA) model is equal to:      ∆ = areas DES in ) , min( areas RANS in ~IZDES DES wall wall I ZDES C d d d (1)
where d wall is the wall distance; C DES is the original DES97 constant [START_REF] Spalart | Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[END_REF] and

ω ω S I ZDES vol I ZDES = ∆ = ∆ ∆ ∆ ∆ = ∆ = ∆ õr z) y x ( ~1/3
, 𝑆𝑆 ̅ 𝜔𝜔 is the average cross section of the cell normal to the vorticity vector ω. In the present study, the subgrid length scale based on the local vorticity vector direction Δ ω was chosen since there is a preferred mean vorticity direction. Furthermore, the near-wall functions of the SA model are removed in LES areas for mode 1:

f v1 =1, f v2 =0, f w =0.
The ZDES mode 2 was developed to deal with separations over smooth surfaces and works as a global approach where the model itself defines the RANS and DES areas, which makes the treatment of complex geometries straightforward for the user. To ensure that attached boundary layers are treated in RANS whatever the grid density, a shielding function f d is used -similar to the one employed in DDES [START_REF] Spalart | A New Version of Detached-eddy Simulation, Resistant to Ambiguous Grid Densities[END_REF]-to shield the attached boundary layers. The hybrid length scale reads: . As explained in Ref [START_REF] Deck | Recent Improvements of the Zonal Detached Eddy Simulation (ZDES) Formulation[END_REF], the nearwall functions of the SA model are not modified in mode 2 of ZDES.

Single stream jet

Description of the test case and simulations

The round jet investigated is similar to one investigated experimentally in Ref. [START_REF] Davoust | Dynamics of m = 0 and m = 1 Modes and of Streamwise Vortices in a Turbulent Axisymmetric Mixing Layer[END_REF]. The Reynolds number based on the nozzle exit velocity U 0 = 20 m/s and diameter D = 0.15 m is equal to Re D = 2.1 10 5 . During the experiments, the boundary layer inside the nozzle was tripped using sand paper so that the exit boundary layer is fully turbulent, of thickness δ 0 = 3 mm = 0.02 D and momentum thickness θ 0 = 0.36 mm = 0.0024 D. The simulations were performed with the elsA software developed at ONERA [START_REF] Cambier | The Onera elsA CFD software: input from research and feedback from industry[END_REF] which solves the compressible Navier-Stokes equations on structured multiblock meshes. A sketch of the computational domain and boundary conditions is presented in Fig. 1. 40 grid points are used to discretize the nozzle boundary layer, 360 points are used in the azimuthal direction and the mesh contains 40.10 6 cells in total.

The time integration is performed using an implicit LU-SSOR algorithm and a second-order accurate backward Gear scheme with 5 Newton sub-iterations to ensure a drop of one order of magnitude of the inner-iteration residuals.

The diffusive fluxes are discretized using a second-order-accurate centered scheme. The convective terms are treated with the hybrid centered/upwind second-order-accurate AUSM+P scheme [START_REF] Mary | Large Eddy Simulation of Flow around an Airfoil Near Stall[END_REF] using MUSCL extrapolation of the third order.

To generate a realistic turbulent inflow content, the Synthetic Eddy Method (SEM) [START_REF] Jarrin | Reconstruction of Turbulent Fluctuations for Hybrid RANS/LES Simulations using a Synthetic-Eddy Method[END_REF] is used. The SEM is based on the generation of velocity fluctuations carried by synthetic eddies which are superimposed to a target mean flow to obtain a synthetic unsteady turbulent field used to feed an inlet boundary condition.

The formulation relies on the prescription by the user of a target Reynolds stress tensor, which is related here to the jet core turbulence rate Tu: 𝑢𝑢 𝚤𝚤 𝑢𝑢 𝚥𝚥 ����� = 𝑇𝑇𝑢𝑢² × 𝑈𝑈 0 2 × 𝛿𝛿 𝑖𝑖𝑖𝑖 (where U 0 is the inlet mean velocity and δ ij is the Kronecker delta). The values tested for the turbulence length scale σ are given in Table 1. For all simulations, the time step is 2 10

Results

In this paper, the focus is put on the major modifications in the jet development induced by the introduction of upstream turbulence. A more complete analysis and validation of these results as well as the comparison between mode 1 and 2 of ZDES and DDES is presented in Ref [START_REF] Gand | Investigation of turbulence development in incompressible jets with zonal detached eddy simulation (ZDES) and synthetic turbulent inflow[END_REF].

The shape and size of the injected eddies in the nozzle inlet plane of simulation ZDES2-SEM can be educed from Fig. 2 (b), to be compared with the uniform inlet in simulation ZDES2. The insets in Fig. 2 depict the turbulence development in the early shear layer. The spatial wavelength and transition length seem to be modified by the upstream turbulence. The effect of upstream turbulence on the jet flow development is presented in Fig. 5. The turbulence rate at the nozzle exit in simulation ZDES2-SEM is low but comparable to the one in the experiments, and this seems to improve the potential core length prediction compared to the simulation ZDES2 with zero turbulence rate. The increase of the turbulence rate in the simulations results in a shorter jet (ZDES2-SEM2). When larger synthetic eddies are injected (ZDES2-SEM3 compared to ZDES2-SEM with the same inlet turbulence rate), they are less dissipated by the grid in the nozzle, thus a higher turbulence rate is achieved at the nozzle exit and the jet is slightly shorter. Velocity spectra along the jet centerline in Fig. 6 show the progressive attenuation of the jet column mode at St D = 0.6 when the turbulence rate is increased. The dual-stream jet considered in this paper is the MARTEL configuration which was investigated experimentally at the Pprime institute in Poitiers, France. This configuration is illustrated in Fig. 7 and consists in an isolated heated axisymmetric coaxial jet, without external flow. The details of the operating point computed are given in Ref [START_REF] Verrière | Zonal Detached-Eddy Simulations of a Dual-Stream Jet[END_REF].

The simulations were carried out with the elsA software [START_REF] Cambier | The Onera elsA CFD software: input from research and feedback from industry[END_REF] on a structured multi-block grid of around 205 million grid points (illustrated in Fig. 8), with 416 points in the azimuthal direction. In the mixing layers, the radial spacing is about ∆r = δω/40, and the axial and azimuthal spacings verify ∆x,z < δω/2, where δω is the shear layers vorticity thickness extracted from a preliminary RANS computation.

Fig. 8 ZDES grid and inlet boundary conditions for the Martel test case

The time integration is identical to the one used in the previous section (see section 3.1). The convective fluxes are treated with the centered second order scheme of Jameson with second and fourth order dissipation coefficients equal to χ2 = 0.5 and χ4 = 0.008.

In order to introduce upstream turbulence, the Random Flow Generation (RFG) technique developed by Smirnov et al. [START_REF] Smirnov | Random Flow Generation Technique for Large Eddy Simulations and Particle-Dynamics Modeling[END_REF] is used at the nozzle internal ducts inflow boundary conditions (see Fig. 8). This method is based on the superposition of harmonic functions and has already been used in elsA simulations in previous work [START_REF] Brunet | Random Flow Generation Technique for Civil Aircraft Jet Simulations with the ZDES Approach[END_REF] [START_REF] Marty | LES and hybrid RANS/LES simulations of turbomachinery flows using high order methods[END_REF]. The parameters investigated in this paper are given in Table 2, the same parameters are used for both primary and secondary duct inlets. The target Reynolds stress tensor is related to the turbulence rate by 𝑢𝑢 𝚤𝚤 𝑢𝑢 𝚥𝚥 ����� = 𝑇𝑇𝑢𝑢² × 𝑈𝑈 0 2 × 𝛿𝛿 𝑖𝑖𝑖𝑖 , where U 0 is the mean velocity at the inlet of the duct considered.

The physical time step used for unsteady simulations is ∆t = 1.5.10 -7 s. The spectral analysis has been carried out using a time signal of T = 36 ms = 116 D/U, where D is the bypass exhaust diameter and U is the bypass exhaust velocity.

Results

Fig. 9 illustrates the salient features of the baseline Jω simulation. The reader is referred to Ref [START_REF] Verrière | Zonal Detached-Eddy Simulations of a Dual-Stream Jet[END_REF] for a detailed validation and analysis of this simulation. The secondary mixing layer is quickly destabilized, whereas the primary one seems to suffer from a delay in the RANS-to-LES transition which may be due to the low velocity difference between the primary and secondary streams [START_REF] Verrière | Zonal Detached-Eddy Simulations of a Dual-Stream Jet[END_REF]. In the secondary stream, a shock-cell patterns is observed. The evolution of the injected turbulence in the nozzle ducts is depicted in Fig. 10. The intensity and size of the velocity field distortions are correlated with the RFG parameters and no major damping of the fluctuations is observed as they are transported within the ducts. The RANS-to-LES transition in the secondary mixing layer can be assessed in Fig. 11. As observed in the results for the single-stream jet, secondary peaks in the profiles of turbulent kinetic energy are somewhat smoothed when upstream turbulence is introduced -even if the main peak is stronger, which suggests a more intense transition. The quicker RANS-to-LES transition leads to a better agreement between the simulations and the experiment regarding the mixing layer growth rate. The azimuthal distributions of vorticity in the jet cross plane at x/D = 1 presented in Fig. 12 indicate the level of resolved turbulence outside the mixing layers in both streams. This results in increased mixing in the streams core, which in turns modifies the shock-cell positioning as shown in Fig. 13 and reduces the phase-lag between simulations and experiments. This is a major improvement in the simulations and shows the necessity of taking into account upstream turbulence for a quantitative jets aerodynamic performances prediction.

Summary and conclusion

In both cases investigated in this paper, the introduction of upstream turbulence at the nozzle inlet in the jets core has led to improved agreement with the experimental data with regards to:

(1) The jet flow development, in terms of potential core length and shock-cell location, thanks to a better reproduction in the simulations of the experimental aerodynamic conditions; (2) The RANS-to-LES transition in the shear layers, which is partly attributed to the penetration of the injected fluctuations inside the nozzle boundary layers which promotes the development of LES content in the shear layer. Of interest, it was not observed that this behavior impaired the shielding of the attached boundary layers treated in RANS mode with the ZDES mode 2 (which provides a RANS/LES switch without any user intervention). Therefore these results promote a systematic introduction of upstream turbulence in eddy-resolving simulations of jets, and emphasize the need for more detailed experimental data on this matter. However, only aerodynamic data has been analyzed so far. Future work will be devoted to the assessment of the jet noise modification due to the upstream turbulence rate.

  to DDES, the mode 2 of ZDES provides a specific definition of the subgrid length scale according to the flow resolution and a threshold value f d = 0.8 determined in Ref [1]. When f d <f d0 , max ~∆ = ∆ II ZDES (Δ max being the characteristic mesh length necessary to ensure the correct behavior of f d ). When f d >f d0 , the subgrid length scale revolves to

Fig. 1

 1 Fig. 1 Single stream jet case, computational domain

Fig. 2

 2 Illustration of the simulations and close-up on the mixing layer early development More quantitatively, the RMS profiles along the lip line in Fig.3show that the dual-peaked shape characteristic of vortex pairing is attenuated. This is attributed to a reduction of the RANS-to-LES transition length when upstream turbulence is introduced in the jet core, due to the penetration of resolved fluctuations into the nozzle boundary layer treated in RANS which promotes the development of LES content in the early stage of the mixing layer development. The attenuation of the vortex pairing is further evidenced by the velocity spectra in Fig.4near the nozzle exit, since the hump near St θ = 0.013 associated with the vortex pairing is smoothed in cases with upstream turbulence.

Fig. 3 Fig. 4

 34 Fig. 3 Evolution of shear stresses along the lip line

Fig. 5 Fig. 6 4 . Dual-stream jet 4 . 1 Fig. 7

 564417 Fig. 5 Evolution of the mean velocity and RMS of velocity along the jet centerline

Fig. 9

 9 Fig. 9 Main features of the Jw simulation

  Jw-T10-H/10 (c) Jw-T5-H/10 (d) Jw-T5-H/5 Fig. 10 Radial velocity fields

Fig. 11

 11 Fig. 11 Fan mixing layer streamwise evolution

Fig. 13

 13 Fig.[START_REF] Verrière | Zonal Detached-Eddy Simulations of a Dual-Stream Jet[END_REF] Time-averaged streamwise velocity in the fan and core streams (red lines in Fig.7)

  -6 s. The unsteady results are time averaged over 200 ms (~25 D/U) after a transient stage estimated of around 120 ms (~16 D/U).

	Case	RANS/LES approach	Inflow condition
	ZDES2	ZDES (all domains : mode 2 Δ ω )	Steady uniform Pt, Tt
	ZDES2-SEM	ZDES (all domains : mode 2 Δ ω )	SEM Jarrin σ=0,02 Tu=7,5%
	ZDES2-SEM2	ZDES (all domains : mode 2 Δ ω )	SEM Jarrin σ=0,02 Tu=35%
	ZDES2-SEM3	ZDES (all domains : mode 2 Δ ω )	SEM Jarrin σ=0,05 Tu=7,5%

Table 1

 1 Simulations performed for the single-stream jet case

Table 2

 2 Simulations performed for the Martel test case. H is the height of the duct considered

Acknowledgments

The simulations presented in section 3 were performed using HPC resources from GENCI-TGCC (Grant 2014-t20142a7215).

Fan

Core


ONERA, The French Aerospace Lab

Wall-modelled Large Eddy Simulation (WMLES) may resolve turbulent fluctuations in the outer layer of attached boundary layers when such a description level is needed (mild flow separation with upstream history effects, dynamic loading and aeroacoustics predictions…). The inner layer dynamics, whose resolution cost is prohibitive at high Reynolds numbers, is modelled ( [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF]). One possible WMLES strategy is a near-wall RANS / outer LES treatment in line with hybrid RANS/LES methods suited for complex geometries and flows. The zones of interest for LES and the zones were a RANS description is sufficient may be defined by the user, as for instance in Zonal Detached Eddy Simulation (ZDES, [START_REF] Deck | Recent improvements in the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF]), with its mode III suited for WMLES ( [START_REF] Deck | Zonal Detached Eddy Simulation (ZDES) of a spatially developing flat plate turbulent boundary layer over the Reynolds number range 3150 ≤ Re θ ≤ 14000[END_REF]). Mean skin friction predictions within 5% of the Coles-Fernholz correlation are achieved in a flat-plate zeropressure-gradient boundary layer with Δx + =200 and Δz + =100 resolutions by properly choosing the RANS/LES interface wall distance, e.g. y = 0.1 δ or y + = 3.9 Re τ 1/2 (fig. 1). The latter interface position provides in addition a higher resolved fraction of C f ([4]). The physical relevance of the resolved turbulent content (see fig. 2) may be assessed from streamwise velocity spectra (fig. 3). Comparing with experimental data ( [START_REF] Deck | Zonal Detached Eddy Simulation (ZDES) of a spatially developing flat plate turbulent boundary layer over the Reynolds number range 3150 ≤ Re θ ≤ 14000[END_REF]) suggests an excess of streamwise velocity fluctuations right below the RANS/LES interface, especially with the 3.9Re τ 1/2 interface. To better understand whether the fluctuations are physically relevant or whether they are spurious artefacts caused by the WMLES strategy, their frequency-dependent convection velocity is determined here (see fig. 4) by a spectral evaluation method introduced in [START_REF] Renard | On the scale-dependent turbulent convection velocity in a spatially developing flat plate turbulent boundary layer at Reynolds number Re θ = 13000[END_REF] in a Wall-Resolved LES case and based on the following equation:

)) ( Im( ) ( 2 ) ( f u x u S f uu S f = f c U    (1)
Two different interface locations will be compared, with special attention paid to the location of the U c = <u> line (analogous to the critical layer in resolvent analysis [START_REF] Sharma | On coherent structure in wall turbulence[END_REF]) relative to the spectral energy sites (fig. 3) in order to discuss the physical nature of the fluctuations.

Fig. 1 Fig. 2

 12 Fig. 1 Mean skin friction coefficient evolution Fig. 2 Streamwise velocity isosurface coloured by wall distance near Re θ = 13 000 (ZDES, 3.9 Re τ 1/2 interface)

Fig. 3 Fig. 4

 34 Fig. 3 Streamwise velocity spectrum, Re θ = 13 000 (ZDES, 3.9 Re τ 1/2 interface, solid line U c = <u>)