professeur, LISPEN M Richard Bearee

Président M Nabil Anwer

M Hicham Nouria

M Adel Olabi

Mohammad Vakilinejad

Olivier Gibaru

Adel Olabi

Nationale Supérieure D'arts

Métiers

Keywords: Industrial Robots, Robotic Manufacturing, Non Geometrical Errors, Compliance Behaviour, Ultrasonic cutting, Honeycomb Material, Workcell Optimization, singularity avoidance

Amélioration de la précision des robots industriels pour la découpe de matériaux composites

Improving the accuracy of Industrial robots for machining process in composite materials composite

List of Abbreviations

Introduction

Progress in various fields of engineering such as electronics and mechanics has extended the horizons of engineering applications. Robots are one of the perfect examples of devices made feasible by simultaneous achievements in once non-relative subjects. Along with the technological and scientific advances, the movement in fictional movies and literature in the late 70s took rule as a second pillar for picturing the future of human lives impacted by robots.

In brief, although the depicted future in which distinguishing humans from their identically structured and emotional robots are yet to come, however, robots are gradually occupying the duties and positions of expert human in many applications.

Industrial robot (IR) arms are a derivative of robots designed and heavily involved in different industrial applications. However, the applications in which industrial robots were involved, have been subjected to changes. Initially, IRs were only employed for applications such as pick and place and painting. These applications are known for their requirement of repeatability in place of accuracy. Previously, the demands in the industrial applications were comparatively unchangeable in time. The geometrical and fundamental aspects of the workpieces were rarely subjected to updates. Therefore, the changeability and adaptability of production lines were not considered as a vital property.

Nowadays on the other hand, based on significant advancement in areas such as numerical simulation software, computational capacity of computers and material science combined with the rapidly evolving market in sectors such as aerospace and transportation, the manufacturing industry is undergoing a revolution to come up with designs and configurations where changeability is the number one feature.

In light of this industrial revolution, special attention has been paid to the IRs because of notable characters inline with the objectives of recent manufacturing needs. They come with a comparatively large workspace, easy to program with a lighter weight compared to conventional machining devices and less costly.

But however, replacing machining devices with their IR competitors does come with certain drawbacks. Notable disadvantages are also involved which introduce further challenges to this fundamental change. Due to the lower accuracy values of IRs in comparison to the machining devices, their involvement in certain machining processes, especially those coming with high and varying machining interactions requires further improvements to satisfy the dictated industrial tolerances.

Chapter 1. Introduction

This PhD thesis is dedicated to this goal in partnership with Le Creneau Industriel company located in Annecy, France. Along with other industrial centers and companies, Le Creneau Indstriel is also willing to step in the terrain of robotinvolved machining production lines. Le Creneau is a producer of CNC machines mainly focused on providing machining systems and solutions for high accuracy aeronautical industry.

Their focused area mainly involves ultrasonic cutting of honeycomb materials, carbon fiber trimming, and hole drilling. The base of the conventional CNC machine structures in the Le Creneau industry is a gantry structure containing three transnational actuators armed with a shoulder having two rotary degrees of freedom to execute 5 axes machining tasks. Le creneau industriel is aiming to vast its applications and reduce the cost of their proposed CNC solutions and shipments by replacing the heavy and large gantry structures with Stäubli TX200 robot arm. Figure 1.1 depicts both bases for a CNC machining device. To investigate the feasibility of replacing the robot arm with the conventional gantry structures, the outline of this thesis is positioned in the following manner:

Chapter 2 surveys different potential sources of inaccuracy in robot-involved machining processes. A literature review of the subject is prepared. The main contributors of robot positioning errors are mentioned. The portion of contribution for three non geometrical error sources, namely compliance behavior, thermally induced inaccuracy and nonlinear transmission errors of robot gearboxes for our case study is evaluated. The evaluation is carried out based on experimentally driven data for loaded and unloaded Stäubli TX200 robotic arm.

Due to the importance of the error sources observed experimentally, two non geometrical error sources, compliance behavior of robot structure and nonlinear errors of robot transmissions systems were chosen for further investigations in chapter 3. Compliance behavior was firstly focused. To compensate for the errors resulted from this error source, an elastostatic model of the case study robot is developed involving different aspects including the self-gravity effect of robot structure, the impact of external load and the self-balancing mechanism attached to the robot. As a vital element of the elastostatic model, the stiffness values of robot joints are experimentally evaluated. The model is then used to compensate for the robot path for a given loading condition. Experimental values are presented to show the efficiency of the elastostatic model performance.

The remainder of the chapter deals with the efforts dedicated to identify and reduce the nonlinear transmission errors of the robot joints. In this part, the periodic errors of a laboratory case study robot is identified using a proposed novel path in Cartesian space massively reducing the error identification process time compared to the suggested methods in the literature. The identified error is used for an offline compensation procedure. The performance and efficiency of the identification and compensation procedures are validated based on experimental results.

As one of the main applications of Le Creneau industriel is providing solutions for ultrasonic cutting of Aramid honeycomb materials, chapter 4 provides a closer look at specific characters of this machining process executed by Stäubli TX200 robotic arm. To thoroughly examine the robot behavior in the machining condition, the machining interaction forces and geometrical error values were quantified by several test sessions. A force model is generated to estimate the machining forces for given features in the honeycomb chamfering process. This model is then used to simulate the deviations and displacement of robot Tool Center Point (TCP) in a chamfering process of a honeycomb material.

Due to the high value of the compliance behavior of the machining tool in this machining process, a novel method is proposed to counter the compliance effect for a robot-tool compliant system. The compliance model in the triangular ultrasonic cutting knife, which performs as the machining tool in this application, is developed using the structural stiffness matrix (SSM) method. The model is examined on-site with experimental values. The most compliant direction of the tool is modeled as a virtual joint of the robot structure. The effect of an external load on this new 7 axis robot is simulated using the equivalent torsional stiffness of the tool, and experimentally achieved stiffness values of the robot joints.

Due to the great influence of the machining tool structure, the loading condition applied to the tool is modeled using a non uniformly distributed load along its length. An offline procedure for decreasing the geometrical errors observed in the honeycomb pieces is proposed based on an iterative algorithm. The procedure uses the generated machining force and tool compliance model. The compensated and uncompensated error values are reported.

Chapter 5 presents a strategy to optimize the placement of a carbon fiber workpiece subjected to a trimming process. The criteria of this optimization problem are to reduce the effect of compliance behavior and nonlinear transmission errors in the final workpiece. A novel cost function based on the variance of the induced errors Chapter 1. Introduction is proposed. To avoid singular configurations of the robot, a new approach is introduced based on penalization of cost function depending on how far a proposed configuration is, from three different singular conditions that a wrist partitioned robot can encounter.

The optimization method employed in this chapter is a modified free pattern search algorithm with generalized grid construction. This grid aides the optimization algorithm for better decision making on choosing the best direction for moving forward. Simulation results are presented to show the performance of both the proposed cost function and optimization method. Figure 1.2 presents the outline of this thesis. Error Sources In Robotic-Based Machining

Improving Accuracy in

Introduction

Robot inaccuracies can rise from different sources. According to handful of research works such as [START_REF] Kim | Robotic Machining: A Review of Recent Progress[END_REF][START_REF] Schneider | Experimental investigation of sources of error in robot machining[END_REF][START_REF] Gong | Nongeometric error identification and compensation for robotic system by inverse calibration[END_REF] these sources can be classified into two different categories which are namely Geometrical and Non geometrical error sources:

• Geometrical Errors

These errors arise from the tolerances of the robot elements manufacturing, assembly, and installation process. They are present in both links and joints thus they always exist regardless of the working condition of the robot. Rather than manufacturing tolerances resulting in deviation of fabricated robot parts from their nominal values, joint, and axes misalignment are also categorized in this area.

• Non geometrical Errors

Other types of errors which can not be compensated by permanently modifying the geometrical model of the robot are termed non geometrical errors. The most important feature of these errors is their dependency on further factors such as applied load and the smoothness of trajectory in provoking compliance behavior of the robot, environmental factors such as time-variant temperature changes causing deformation of robot links or configuration dependent errors such as irregularities in joint variables.

Reducing geometrical errors of a robot is usually based on updating the mathematical model serving as forward kinematics by using experimentally obtained offsets instead of nominal values for robot manipulator Nubiola and Bonev, 2013. Geometric calibration methods can differ in various manners. Mathematical model adopted for base-to-end frame transformation (Direct geometrical model) can vary from updating conventional Denavit-Hartenberg parameters (also called DH parameters) to exponential transformation in Lie Theory used in [START_REF] Fu | Analysis of unified error model and simulated parameters calibration for robotic machining based on Lie theory[END_REF]. However, the idea behind all the proposed approaches is to compare the actual placement of the robot end-effector with the desired placement calculated from the nominal geometrical model. The difference between the desired and actual values are used to permanently update the geometrical model of the robot to overcome geometrical errors. The literature review presented by [START_REF] Wu | Geometric calibration of industrial robots using enhanced partial pose measurements and design of experiments[END_REF] can serve as a guideline for Chapter 2. Error Sources In Robotic-Based Machining inclined readers in geometrical calibration topic.

Even though the impact of Geometrical errors of the robot is undeniable in robotic based applications, a glance on the literature reveals that this error source is mainly regarded as the main contributor for an unloaded situation of the robot. [START_REF] Schneider | Experimental investigation of sources of error in robot machining[END_REF] report a contrasting result out of their investigations on a robotic arm executing a machining task. Moreover, the continuous progress in the fabrication methods offers more accurate elements with lower geometrical tolerances promising decreased geometrical error values.

A generic classification for non geometrical errors is seemingly impossible due to the vast number of sources fulfilling its definition. Furthermore, these sources are not limited to the mechanical aspects of an Industrial Robot (IR) system. Other sources originating from electrical and programming aspects also take part to a certain level in contribution to non geometrical errors. The effect of limited bandwidth mentioned in [START_REF] Kim | Robotic Machining: A Review of Recent Progress[END_REF] is an example of a non geometrical error possessing a nature other than mechanical aspects. However, from a mechanical perspective to the non geometrical errors affecting robot accuracy, the noteworthy error sources pointed out by the majority of researchers can be classified as illustrated in Figure 2.1.

FIGURE 2.1: Classification of Mechanical Robot Pose Errors

In light of the above-mentioned points, in this chapter of the thesis, the effect of non geometrical errors is investigated to reveal their effect and contribution in the misplacement of robot tool center point (TCP). The remainder of the chapter is organized as follows:

In section 2.2, we took a closer look at the compliance behaviour of the robot. A literature review of the origin and effect of this error source is carried out. The compliance behaviour is divided into two main subcategories which are self gravity effect and external couple effect. Using experimentally driven data of a machining 2.2. Compliance Behaviour of Robot 7 robot case study, the portion of each source in the total misplacement of TCP due to compliance behaviour is determined.

Section 2.3 concerns the effect of irregularities observed in the robot placement due to cyclic errors arising from the transmission system of robot joints. State of the art in this subject is presented to ensure the validity of error origin according to the literature. The amplitude of this error type in the condition of the designed experiment is evaluated.

Section 2.4 is dedicated to the time-variant thermally induced errors. A literature survey on this topic has been reported. An experiment is executed to verify the effect of thermal changes in the accuracy of the robot. The evolution of temperature on the robot shell has also been reported. The chapter ends with noteworthy comments and conclusions and an illustration of portions of the contribution of studied effects throughout the chapter.

Compliance Behaviour of Robot

Replacement of conventional machining tools with IRs while offering advantages such as extended workspace and more importantly lower prices, introduces new obstacles. In the wide range of publications and research articles contributing to the robot accuracy improvement, it is hard to name an item in which authors have not mentioned the compliance behaviour of robot as one of the main contributors to the robot pose inaccuracy if not the most important among them. The effect of compliance behaviour can be responsible for more than 70% of robot inaccuracies [START_REF] Bu | Stiffness analysis and optimization in robotic drilling application[END_REF].

The compliance behaviour of robot originates from the deformation of robot structure due to link weights and external couples applied on the robot. These loads are then translated into axial and torsional tensions along the robot body and thanks to the serial configuration of robot the resulted strains magnify by marching towards the end effector of robot.

Among different elements in robot structure, the compliance effect of joints is believed to be the main contributor to the total compliance behaviour. In Rezaei and Akbarzadeh, 2018 authors have pointed that according to their simulations on compliance behaviour of a parallel robot, the effect of joint deformation in the total misplacement of the robot is about 7 times greater than the one of links. Results presented in [START_REF] Klimchik | Robotic manipulators with double encoders: accuracy improvement based on advanced stiffness modeling and intelligent control[END_REF] show that by only eliminating the joint deflection errors authors were able to compensate at least more than half of the compliance error in the robot workspace.

Each revolute joint of an IR is composed of two main elements. An actuator which is a rotary motor armed with an encoder to send feed back to its corresponding control unit and a transmission system which reduces the angular velocity and displacement of motor shaft to the robot link by a specific constant ratio. Due to this arrangement the angular deflections from the link side and compliance of transmission system stay hidden from the encoders responsible for joint variable accuracy [START_REF] Kim | Robotic Machining: A Review of Recent Progress[END_REF] Thus in some applications to achieve the required accuracy, robot In light of above mentioned facts, one can conclude that although elastostatic behaviour is among most important error sources in IR applications, but it is also case sensitive. Different robot arms are designed and manufactured with different elements, geometries and materials thus they come with different weights. The different gear boxes implemented in transmission systems differ in stiffness and reduction ratio. These facts will result in drastic differences in elastostatic behaviour of robots from one another.

To evaluate the effect of elastostatic behaviour of our case study robot Stäubli TX200, we have conducted an experiment which is explained in detail in the following section.

Compliance Behaviour of Stäubli TX200

As described previously, the compliance behaviour rises from two different sources. One being the effect of the weights of robot structure on pose error of Tool Center Point (TCP) which hereafter is termed as Self-Gravity effect. The other source is the effect of external forces applied on the robot end effector which referred to as External Couple effect. It is trivial that based on the mass of robot and the magnitude of exerted external load, the amount of contribution of each source vary.

One of the main sources of external couple is the weight of additional instrument attached to the robot flange. Even though this part can be regarded as self-gravity effect by assuming that the instrument is a part of robot body, but considering the fact that for each machining process this instrument is subjected to change, in this work this part is regarded as a part of external couple.

To examine the effect of each source on Stäubli TX200 robot arm, robot was programmed to follow a 1000(mm) × 1200(mm) rectangle lied on a plane parallel to XY plane of robot base frame (Z = 448.5(mm)). robot was programmed to follow the path with constant orientation and Cartesian speed of 1 (cm/sec) to avoid any excitation of dynamic behaviour of robot structure during the test. The path was followed in two different conditions of the robot. Once with no external load to verify the self gravity effect and once with a constant load of 500(N) in gravity direction. Configuration of robot and force imposition along with the TCP placement captured by a Laser Tracker (LT) are presented in Figure 2.3. According to Figure 2.3, the compliance effect of robot is highly dependent on the TCP location. For both loading conditions the displacement of TCP is less than half on the closer side of the rectangle to robot than the further side. This is evidently due to the lever effect of robot for comparatively more stretched configurations. But since one of the advantages of using IRs, is their extended workspace, eliminating such zones of the robot workspace for improving its accuracy in machining tasks, is not an effective solution especially for processes dealing with large scaled workpieces.

Figure 2.3 also illustrates that in this loading condition, the self-gravity effect on robot is higher than external couple effect. In extreme conditions in the test, the external couple effect is responsible for 1 (mm) added to the self gravity effect which accounts for self gravity effect. This is due to the massive elements and links utilized in this robot which on one hand boasts the total stiffness of the robot structure but however imposes a large self-gravity effect. But it is worth mentioning that this robot is equipped with gravity balancing system attached to joint 2. This system massively reduces the torque resulted from structural gravity over joint 2. Further more, an additional weight compensation option able to counter the gravitational effect of additional mass on robot flange is implemented in the robot controller. More details on the self gravity compensation system and additional weight compensation option are explained in Chapter 3.

Cyclic Joint Transmission System Errors

During the experiments on the compliance behaviour of the robot presented in Figure 2.3, some sort of fluctuational displacement of TCP distributed along the path were noticed. These errors were not originated from the dynamic behaviour of robot Chapter 2. Error Sources In Robotic-Based Machining since the movement of robot was slow. The errors were not affected by the load applied on the robot end effector. Figure 2.4 depicts a closer view of the actual path of robot for both loading conditions. Different researchers have pointed out that these errors are coming from the deviation between nominal and actual values of the joint variables. This deviation comes from the fact that the transmission system laying on the way of the input and output signal of the joint variable is not ideal. Structural compliance, backlash, hysteresis and friction are among the most important causes of the nonlinear behavior of transmission systems [START_REF] Kircanski | An experimental study of nonlinear stiffness, hysteresis, and friction effects in robot joints with harmonic drives and torque sensors[END_REF][START_REF] Oh | Influence of the joint angular characteristics on the accuracy of industrial robots[END_REF] .

Industrial robots can wield variety of transmission systems depending on their application such as harmonic-drives, Cycloidal drive and gear trains using several stages. Harmonic drives were developed in 1995 for aerospace applications but their high-ratio and comparatively small size caused them to spread in many applications in technology. Harmonic drives employ a non-rigid gear called flexspline for speed reduction whose stiffness is lower than that of conventional transmissions containing nonlinear relation between the input and the output torques which makes it more challenging to control the system [START_REF] Taghirad | Modeling and parameter identification of harmonic drive systems[END_REF]. Due to the vast usage of harmonic drives their behavior have been studied and modeled by many researchers addressing different effects such as hysteresis and friction [START_REF] Tuttle | Understanding and modeling the behavior of a harmonic drive gear transmission[END_REF][START_REF] Tuttle | A nonlinear model of a harmonic drive gear transmission[END_REF][START_REF] Kircanski | An experimental study of nonlinear stiffness, hysteresis, and friction effects in robot joints with harmonic drives and torque sensors[END_REF][START_REF] Taghirad | Modeling and parameter identification of harmonic drive systems[END_REF][START_REF] Dhaouadi | A new dynamic model of hysteresis in harmonic drives[END_REF][START_REF] Kennedy | Modeling and control of the Mitsubishi PA-10 robot arm harmonic drive system[END_REF][START_REF] Yamamoto | Modeling and compensation for angular transmission error in harmonic drive gearings[END_REF][START_REF] Preissner | A high-fidelity harmonic drive model[END_REF][START_REF] Tjahjowidodo | Theoretical modelling and experimental identification of nonlinear torsional behaviour in harmonic drives[END_REF][START_REF] Li | New approach for modeling and testing of harmonic drive in robotic systems[END_REF][START_REF] Liu | Nonlinear attributes modeling and analysis of harmonic drive manipulator joint[END_REF].

Another broadly used reducer type in robotic applications is Cycloidal drive (or cycloidal gear reducers). They are commonly used in equipment where precise output and large drive payloads are needed. Recently, with the increasing demand of high efficiency and high speed reduction and torque ratio transmission devices in industry, employment of cycloidal drives have become popular in the automation field as robotics, machine tools and automatic machinery [START_REF] Kao | Computer-aided kinematic error analysis of a two-stage cycloidal drive[END_REF]. Cycloid drives rely on an eccentric motion to convert the rotation of the input shaft into a wobbly cycloidal motion of the planet wheel; this motion is then converted back into a rotation of the output shaft concentric to the input shaft by means of some cylindrical pins fixed on the output shaft and in this way the speed reduction is accomplished [START_REF] Gorla | Theoretical and experimental analysis of a cycloidal speed reducer[END_REF].

Research over Cycloid drives started at 1956 by Bossier and Kingston in [START_REF] Botsiber | Design and performance of the cycloid speed reducer[END_REF] As Cycloid drives use epitrochoid [START_REF] Yang | Design and application guidelines for cycloid drives with machining tolerances[END_REF] curves other than involute gears, manufacturing tools are not yet specialized as for manufacturing involute gears [START_REF] Kao | Computer-aided kinematic error analysis of a two-stage cycloidal drive[END_REF]. Thus many researcher have dedicated their efforts to develop the profile generation of these drives using different geometrical approaches (refer toYang and [START_REF] Yang | Design and application guidelines for cycloid drives with machining tolerances[END_REF][START_REF] Pollitt | Some applications of the cycloid machine design[END_REF][START_REF] Yan | Geometry design of an elementary planetary gear train with cylindrical tooth-profiles[END_REF][START_REF] Lai | Geometric design of roller drives with cylindrical meshing elements[END_REF][START_REF] Hsieh | Study on Geometry Design of Rotors Using Trochoidal Curve[END_REF] or other aspects such as undercutting avoidance (refer to [START_REF] Sensinger | Unified approach to cycloid drive profile, stress, and efficiency optimization[END_REF]. Latest modification for tooth profile of cycloid drive has been done in [START_REF] Ren | Tooth modification and dynamic performance of the cycloidal drive[END_REF] where they modified the gear geometry in accordance to the contact stress. Output angular and speed error of cycloid drivers have been studied mostly analytically by several researchers (refer to [START_REF] Blanche | Cycloid drives with machining tolerances[END_REF][START_REF] Huang | On the Surface Design, Tooth Contact Analysis, and Optimum Design of Cycloidal Drives with Modified Tooth Profiles[END_REF][START_REF] Kao | Computer-aided kinematic error analysis of a two-stage cycloidal drive[END_REF]. In a recent study in Pham and Ahn, 2017 have analyzed hysteresis characteristics of a cycloid reducer using a nonlinear spring with dead zone with aide of finite element method. [START_REF] Li | Analysis of a Cycloid Speed Reducer Considering Tooth Profile Modification and Clearance-Fit Output Mechanism[END_REF] have proposed an analytical model to simulate the performance of cycloid speed reducer taking into account the presence of clearances and eccentricity errors among with elastic deformations.

Thermally induced errors

During the task execution of IRs specially with relatively high rotary speeds, the internal heat sources such as motors, bearings, and ambient heat sources will cause the robot arm structure to slightly deform, which will result in not only link expansion but also structure distortions [START_REF] Yin | Real-time thermal error compensation method for robotic visual inspection system[END_REF] The conventional gantry machining devices usually come with heavy elements wrapped up in box-shaped thick armors which increase their thermal resistance against the ambient. In other words, once they reach thermal stability with the environment in working conditions, the temperature changes in the structure will be negligible. But on the other hand, IRs composed of stretched out links covered with comparatively thin armors are thermally prone to environmental temperature changes and reaching stability in working condition [START_REF] Heisel | Thermal behaviour of industrial robots and possibilities for error compensation[END_REF]. The expansion and distortion of robot links may not be solely in great of significance but the chain structure of articulated robotic arms accumulate and magnify the geometrical changes from the base to the TCP pose.

In a study on thermal behavior of robots done in Heisel, Richter, and Wurst, 1997, has been shown that for an IR with an endurance of 120(kg) of nominal loading and a working envelope of approximately 5 meter of diameter, the repeatability of robot which was specified to be 0.03 (mm) by its manufacturer increased up to 1 (mm) due to thermal effects in constant ambient temperature. They have also shown the significant influence of the joint velocities on the thermal behaviour of the robot. But on the other hand their results revealed no consistency between the payload of the robot to the inaccuracies arising from thermal effects. This drastic effect has caused Chapter 2. Error Sources In Robotic-Based Machining researchers to come up with compensation techniques vastly different from one another. A handful of approaches are reported below. [START_REF] Yin | Real-time thermal error compensation method for robotic visual inspection system[END_REF] have proposed a real-time compensation method for thermally induced errors. The method does not require any thermal sensor and is based on measuring the center point of a standard sphere using a laser sensor attached to the end-effector of the robot. Several geometrical parameters of the robot are then updated based on the measured error. Although they have improved the max/mean performance error of their IR case study from 0.441(mm) to 0.136 (mm), but the method requires a constant real time functioning of a laser sensor attached to the robot end effector.

In [START_REF] Gong | Nongeometric error identification and compensation for robotic system by inverse calibration[END_REF], authors have armed the robot arm with thirteen temperature sensors. By executing some experiments, they used a principle component analysis (PCA) to generate an empirical model able to correlate robot parameter thermal errors with the corresponding temperature field. But one note that this method uses cabling-required sensors, real-time measurements and computations beside offline model generation.

Research upon thermal errors is not limited only to Industrial Robots. In [START_REF] Eastwood | Compensation of thermal deformation of a hybrid parallel kinematic machine[END_REF] have improved performance of a typical HPKM (Hybrid Parallel Kinematic Machine) by combining a generic kinematic modelling 'toolbox' with a real-time thermal monitoring system. Assuming each strut to be rod with a heat source at the end, they have proposed a 'CTE factor' to the theoretical expansion equation to overcome the problems such as the effect of shrouding of the struts. Having an assumption of length extension of each rod by real-time measurement of temperature of every heat source and using an iterative inverse kinematic model for compensations, progressive thermal error was limited to 0.03 (mm).

Although we are seeking accuracy of industrial robots, however due to the undeniable importance of absolute accuracy of the measurement machines, the thermal behaviour treatment in this research area specially those of similar structural configuration such as Articulated Arm Coordinate Measuring Machines (AACMM) might be inspiring for adoption in IR applications. [START_REF] Santolaria | Calibration-based thermal error model for articulated arm coordinate measuring machines[END_REF] have proposed a method to calibrate thermal errors for an AACMM. An empirical model characterized the error behavior with respect to the temperature by analyzing the measurement arm error when measuring a standard artifact placed in several positions. Their method uses two temperature sensors, one attached to the body of the arm and the other measures the ambient temperature. Different geometrical parameters were modeled using polynomial functions of the temperature, covering the nominal working temperature range of the arm. The major difference of these application compared to IRs is the missing of internal heat source such as actuators making it feasible to use only one temperature sensor attached to the measurement arm to fully define the temperature along their structure.

Using Finite Element Method (FEM), Li and Zhao, 2016 created a deformation error model of robot axes in different temperature condition. They have concluded that structural deformation of robot changes linearly as function of ambient temperature. To model the effect of internal heat source effect, they measured the steady state temperature of the robot shell on the zones shrouding drives. Having the values of steady state temperatures and neglecting the heat distribution to the environment, they generated a model defining temperature distribution on the surface of body. Despite of all the efforts put to generate CAD models, for the dynamic compensation, they used a method mathematically similar to the one of Yin et al., 2014. Poonyapak and[START_REF] Poonyapak | Towards a predictive model for temperature-induced deformation of an industrial robot[END_REF] used an innovative way to eliminate the warmup cycle temperature-induced errors when an industrial robot is not operated in a thermally stable condition. Applying their method on a KUKA KR-15/2 Industrial robot, They used two charge-coupled device (CCD) cameras pointed by a laser diode attached to the end-effector of the robot. The coordinates of the points of the collision of laser ray on two CCDs were translated into the geometrical dimensions of the robot using robot geometrical model. On each pose of the warm-up cycle, they pictured the robot with a thermal camera yielding the temperature distribution on robots body. The method is based on empirically estimating a linear coefficient of thermal expansion for two main links. Several important observations of their experiment are as follows: They noticed that the heating curves for different locations on the robot surface are all first-order exponential heating curves. Allowing robot arm to cool-down for 10 hours was not sufficient for the robot to fully thermally stabilize in ambient conditions.

In measurements conducted in [START_REF] Kluz | Effect of temperature variation on repeatability positioning of a robot when assembling parts with cylindrical surfaces[END_REF] the researchers have concluded that for small variations in the ambient temperature changes, the pose error behaviour of the robot acts randomly and thus difficult to compensate. However, the have pointed out that by changing the task location in the workspace of the robot to appropriate positions, thermally induced errors will decrease.

Another way to reduce the effect of temperature on the robot's performance is to modify or optimize the trajectory of robot in order to minimize the generated heat (loss) or time of a robot cycle with the existing trajectory generation lows (refer to [START_REF] Hollerbach | Dynamic scaling of manipulator trajectories[END_REF][START_REF] Betts | Survey of numerical methods for trajectory optimization[END_REF][START_REF] Lin | Formulation and optimization of cubic polynomial joint trajectories for industrial robots[END_REF][START_REF] Bestaoui | On-line reference trajectory definition with joint torque and velocity constraints[END_REF][START_REF] Olabi | Enhanced trajectory planning for machining with industrial six-axis robots[END_REF][START_REF] Luca | A sensitivity approach to optimal spline robot trajectories[END_REF][START_REF] Bobrow | Time-optimal control of robotic manipulators along specified paths[END_REF][START_REF] Swevers | Contributions to Contact Modeling and Identification and Optimal Robot Motion Planning[END_REF] which subsequently results in a lower temperature in steady state situation. In [START_REF] Guilbert | Optimal trajectory generation for manipulator robots under thermal constraints[END_REF][START_REF] Guilbert | Optimization of complex robot applications under real physical limitations[END_REF] researchers proposed a method to generate the optimized velocity profile with a minimum time criterion subject to thermal constraints.

Experimental Investigation

In light of above surveyed methods in the existing literature, one can conclude that compensation of thermal errors usually require real time implementation of external sensors and their corresponding mathematical models of compensation strategies. At this point we have decided to conduct experiments on our case study robot (Stäubli TX200) to observe the order of its thermally induced errors to verify the necessity of compensation o this error source at the mentioned expenses.

To observe how likely our case study IR is to retain its accuracy in a working shift due to internal heating effect, we programmed the robot to constantly move its joints for a duration of about 6 hours at the highest possible angular speed back and forth. 10 points were spotted on the exterior of the robot armor for temperature Chapter 2. Error Sources In Robotic-Based Machining data capturing. These points were chosen very close to the motors and transmission system of the joints which are considered the main heat sources of in the robot structure. Figure 2.5 shows the location of these ten points marked on the robot. During the test session, robot was commanded to take five positions in the workspace. Each time the position of the end-effector and the temperature of the aforementioned ten points were recorded.

Measurement Devices

To measure the Cartesian position of end-effector of robot, we used a Laser Tracker 1 . The LT used in these measurements, was Radian model from API company . This device provides uncertainty of measurements about 20 µm in the range of our experiments. Temperature was measured with an IR-360 infrared thermometer manufactured by VOLTCRAFT providing a precision of 98% of the measured value ±2( o C) (Corresponding to a maximum deviation of ±3.2( o C) regarding our recorded data).

Temperature variation on robot Armor

During the continuous junk movement, at certain time spots, robot was stopped to measure the temperature of the marked points. Initial temperature of both workshop ambient and robot structure was 22.5( o C) stating that robot was initially in complete thermal equilibrium with environment. Throughout the observation session, fluctuations in ambient temperature did not exceed 1.5 ( o C).

Figure 2.7 reveals several noteworthy points. As one can see the temperature changes in all the points follow the general exponential of thermal stability graph of a mass having temperature difference with the environment. The thermal stability can be defined by the time constant of the curve which is the moment when the By a closer look to Figure 2.7 it can be concluded that measurement points attached to the comparatively smaller actuators (such as point 10 attached near actuator of joint 6) possess lower time constants compered to the ones attached to the bigger joint actuators(points 2 and 3).

By comparing the difference between curves attributed to points 2 and 3 in Figure 2.7 attached to link 1 and those of points 4 and 5 attached to link 2, we can conclude that unlike the box-shaped links such as link 1, for lengthy links such as link 2, the temperature distribution is affected only by the neighbouring actuators in each end. Thus a non uniform distortion of the link geometry is expected in the working condition of links of this type.

Another conclusion drawn from from Figure 2.7 is than even more than 6 hours Chapter 2. Error Sources In Robotic-Based Machining of constant movement of robot joints are not enough to reach a complete thermal stability especially specially for the points near principal joints (joints 1,2 and 3).

Repeatability Measurement Approach

To evaluate the repeatability of our robot, five spots were chosen representing a rectangle and its midpoint in a promising zone of the robot workspace (Figure 2.8). Each time that robot was temporarily stopped to measure the temperature of the selected points on its structure, robot was also commanded to reach the Cartesian positions defined in table 2.1 (values defined in robot base frame). To report the inaccuracy induced from the thermal behaviour of robot structure, the achieved position from the measurement of each point is compared to the position of the corresponding point in the first measurement. This allows to eliminate the errors rising from other sources such as geometrical imperfections and self gravity of the robot structure. Thus the repeatability at time t (R(t)) is defined as the distance between the measured position at instant T and the one measured prior to the robot movement (t = 0):

R i (t) = (x i (t) -x i (0)) 2 + (y i (t) -y i (0)) 2 + (z i (t) -z i (0)) 2 (i = 1 : 5) (2.1)
The evolution of repeatability parameter defined above is depicted for the chosen five stops is depicted in Figure 2.9. FIGURE 2.9: Evolution of position Repeatability of robot as function of Warming Time Figure 2.9 shows that for a zone located in the prime of the workspace of robot, the increase of inaccuracy attributed to the thermal behaviour of robot is less than 0.3 (mm). Moreover it reveals that for the first hour of test execution, the induced inaccuracy is about robot nominal repeatability (0.06(mm)) thus we can conclude that for this period, thermal behaviour does not practically impose any inaccuracy.

We should emphasize that this test was executed under extreme conditions regarding joint actuators velocities and hence we can expect slighter increase of robot inaccuracy for machining tasks which usually demand low values of angular velocities in the joint space.

Conclusion

Throughout this chapter, we studied the effect of non geometrical errors on the IR pose aimed to perform machining processes. Compliance behaviour of robot structure was studied in Section 2.2. It was mentioned that the compliance behaviour for static and quasi-static robot tasks can be divided into two different categories which are namely self-gravity effect which is due to the imposed torques on the robot joints to withstand the weight of robot links and the external load effect which is the contribution of the applied load on TCP to the joint deflections.

During the experiments carried out on the compliance behaviour of the robot, we noticed cyclic irregularities that were distributed all along the desired path. These cyclic errors were independent of the external path. By attributing these errors to the misalignment and gear wear effect in the joint transmission systems, we concentrated our attention on a closer look at this error type in Section 2.3. This chapter was then followed by a literature review about the sources of this error and more generally the topic of imperfections in the transmission system behaviour.

Chapter 2. Error Sources In Robotic-Based Machining

To observe the effect of temperature variation along the robot structure and its accuracy due to task execution, in Section 2.4 robot apparatus was programmed to continuously move its joints idly at maximum speed and the effect of generated heat by the robot actuators was reported in terms of both temperature changes in several points on the body and augmented inaccuracy due to the thermal expansion of the links.

According to the analysis of the gathered data in the simulations and experiments reported in this chapter, for a Stäubli TX200 industrial robot arm executing a quasi-static task with an external vertical load of 500(N) (which is not an unrealistic assumption while an ordinary spindle weights around 40(kg)) and a configuration illustrated in Figure 2.3, the portion of contribution of each surveyed non geometrical error source is stated in the pie chart shown in Figure 2.10. The conclusion drawn from Figure 2.10 is the vivid influence of the compliance behavior on robot pose displacement which accounts for about 87% (Self-gravity and external load effect). According to the figure, The effect of compliance behaviour is 7 times greater than the remaining sources combined.

Chapter 3

Identification And Compensation Of Non Geometrical Robot Errors

Introduction

Adopting robot arms for machining processes comes with both advantages and challenges. One of the main challenges of IRs utilized for machining tasks is their inaccuracy compared to their conventional alternatives. As concluded in chapter 2, the main sources of non geometrical errors in order of their amount of contribution, are compliance behaviour (resulting from both self gravity and external load effects), cyclic errors attributed to joint transmission systems and thermally induced errors.

Due to the importance of these error sources, this chapter is focused on identification and compensation of the first two mentioned non geometrical error sources, namely Compliance behaviour and Cyclic joint transmission system errors. These sources were responsible for about 96% of positioning error of robot TCP in the experimental conditions stated in chapter 2.

To overcome the compliance effect of robot structure, we propose an elastostatic model capable of estimating and consequently compensating the static and quasistatic positioning errors of robot TCP in section 3.2. The model involves both self gravity and external load effects. The section describes the mathematical modeling along with experimentally determined properties such as joint stiffness values. The generated model has then been put into practice to practically observe its efficiency in a real loading condition. The section ends with the introduction of a novel way of joint stiffness identification based on robot performance in quasi-static conditions.

Cyclic joint transmission errors are studied in 3.3. State of the art in identification and compensation of this error source is presented firstly and throughout the section, we have introduced a novel method based on an optimal trajectory for robot TCP to identify and mathematically model this error. Later on, an offline strategy is proposed to reduce the effect of cyclic transmission joint errors. The identification model and compensation strategy are applied to a case study robot to validate their effectiveness.

Modeling and Compensation of Compliance Behaviour of Robot

Series structure of IRs despite of offering a magnificent versatility for different applications induces a major drawback for the accuracy of robot. Inaccuracy resulted 20 Chapter 3. Identification And Compensation Of Non Geometrical Robot Errors from the flexibility of the joints and links cause a large pose inaccuracy due to selfgravity and machining forces and remain the major obstacle to widespread use of robots for precise machining applications.

This section consists of first a literature review over different approaches taken by researchers for modeling and compensating compliance effect. Further comes explanations of our generated elastostatic model , at last we present practical results for model verification.

State of The Art

According to [START_REF] Marie | Elasto-geometrical modeling and calibration of robot manipulators: Application to machining and forming applications[END_REF], there are two main approaches to overcome the compliant errors of IRs. First approach is to change the physical design of the robots to use shortened links, increased sections and parallelogram loops as Lijin and Longfei, 2017 have proposed. Trivially this approach is a costly one especially when the robot end users lack the possibility or authorization to modify any physical aspect of the robot. Second approach is to calibrate the robot to handle the errors of this origin. Calibrating for joint and link compliance is based on generating a model which takes into account the elasticity behavior of the robot.

The effect of the inherent compliance of the robot causes two main imperfections to its performance. Deflection of the TCP in both position and orientation in static and quasi-static working mode and vibrations and fluctuations in dynamic working conditions. Enormous efforts have been put to model industrial robots with their compliant behavior. Depending on the complexity and approaches of such methods, [START_REF] Klimchik | Stiffness modeling for perfect and non-perfect parallel manipulators under internal and external loadings[END_REF] have categorized the existing methods into the following three main groups: Finite Elements Analysis (FEA), Matrix Structural Analysis (MSA) and Virtual Joint Modeling method (VJM). We make use of their state of art survey to briefly explain each group:

• Finite Element Analysis (FEA)

In this approach, robot structure is studied with its actual geometry and material. Using a CAD software each robot element is divided into numerous elements and nodes. By applying equilibrium equation on each element, the reaction forces and torques in each node is computed and using the material properties, the displacement and deformation of each element is computed. Thus the behaviour of entire robot structure can be simulated for any given loading condition.

The advantage of this approach is its accuracy and the thorough information about the most critical zones of robot in terms of strains and deformations. But the computational cost required by this method has prevented this model to be widely used in robot compliance studies.

To overcome the computation cost problems of this approach while benefiting from its accuracy, [START_REF] Klimchik | CAD-based approach for identification of elasto-static parameters of robotic manipulators[END_REF] have proposed a way to apply the FEA simulations on different elements of robot in order to achieve the stiffness matrix of the corresponding element. Thus this matrix can be used by drastically lowering computational costs with an acceptable accuracy.

In either cases of using FEA model however, a detailed information of robot elements (geometry and constructing materials) must be in hand which are rarely shared with the robot's end user.

FEA method has been used more to model the behavior of parallel robot mechanisms since they are simpler in terms of composing elements and comparatively more accurate applications [START_REF] Bouzgarrou | Rigidity analysis of T3R1 parallel robot with uncoupled kinematics[END_REF][START_REF] Briot | Global topology optimization of industrial robots with the linearization method[END_REF][START_REF] Corradini | Evaluation of a 4-degree of freedom parallel manipulator stiffness[END_REF]Krut, 2003 Klimchik, Pashkevich, and[START_REF] Klimchik | CAD-based approach for identification of elasto-static parameters of robotic manipulators[END_REF][START_REF] Rezaei | An investigation on stiffness of a 3-PSP spatial parallel mechanism with flexible moving platform using invariant form[END_REF][START_REF] Wang | Dynamics Model of 4-SPS/CU Parallel Mechanism with Spherical Clearance Joint and Flexible Moving Platform[END_REF][START_REF] Cao | A method for stiffness analysis of overconstrained parallel robotic mechanisms with Scara motion[END_REF][START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF].

• Matrix Structural Analysis (MSA)

In this method, instead of using the true shape of robot elements used in FEA method, elements are presented by simple compliant elements such as beams and their stiffness matrix relating the displacement of the element to the applied load. The size of global stiffness matrix and consequently the computation complexity is dependent of the degrees of freedom attributed to the composing element [START_REF] Rezaei | Compliance error modeling for manipulators considering the effects of the component weights and the body and joint flexibilities[END_REF]. Marie, Courteille, and Maurine, 2013 have presented the efficiency of this approach for an industrial robot arm for a practical task. Klimchik, Pashkevich, and Chablat, 2019 have presented the fundamentals of this approach for robotic application.

• Virtual Joint Method (VJM)

Mathematically simpler among all methods, VJM method assumes robot links to be completely rigid. Thus according to this method, the robot compliance behaviour can be modeled by only joint deformations. Robot joint are modeled as linear torsional springs which deform only about its corresponding joint axis. Due to this assumption, the stiffness matrix used in this approach has the same size as the number of joints. However in some research works, not all the robot joints are considered to be rigid which result in reduced order of stiffness matrix [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF].

Although this approach requires fewer computations compared to previous methods, but its effectiveness in reducing robot compliance pose errors, offers the optimum cost-efficiency for wide applications in robotics. [START_REF] Salisbury | Active stiffness control of a manipulator in cartesian coordinates[END_REF] was the first to introduce VJM method for industrial robots and transformation between joint and Cartesian stiffness. But this job was neither focused on the error compensation coming from the flexibility nor stiffness identification of the robot. He introduced a method for online force control in which the desired stiffness values of the joints will be updated during a robot task. This approach requires direct communication to the motor drives to command torque values which the possibility does not exist in our application. The effect of joint flexibility on the pose error was first addressed, experimented and calibrated in [START_REF] Whitney | Industrial robot forward calibration method and results[END_REF]. They determined the stiffness of joints 2 and 3 of a PUMA 560 serial robot. They applied several forces and recorded the displacement of desired points by a theodolite to estimate for joint compliance.

Even though it has been more than three decades of the introduction of this method, yet many recent research items have adopted this model for modeling robot compliance behaviour [START_REF] Du | Compliance error calibration for robot based on statistical properties of single joint[END_REF][START_REF] Theissen | Closed-forceloop elastostatic calibration of serial articulated robots[END_REF].

Elastostatic Model Generation

This section is dedicated to generate a mathematical model that is able to calculate the torque, angular deviation of joints and displacement of TCP due to self gravity and external forces based on VJM methodology. The model is subdivided into three parts: part one is the explanation of the effect of gravity on each joint for a base downward (floor mounted) robot version. Part two states the modeling of the balancing system on joint number two of the robot and the last part computes the induced torque regarding external forces.

Self Gravity Effect of Robot Structure

As shown earlier in chapter 2, the effect of self gravity on our case study robot is more than the effect of the rest of non geometrical error sources combined. Thus in the elastostatic model generation, the first part is dedicated to this very error source.

As the self gravity effect is mainly due to the imposed torques resulted from the link weights, it is necessary to use the exact values of link masses and their center of gravity (CoG). These information have been provided by Stäubli company and due to confidential issues are not presented in this thesis. However the provided coordinates of CoGs were defined in their own defined coordinate systems which are attached to the actual placement of actuators (Figure 3.1) rather than following DH frame construction strategy.

First step to generate this part is then to define the mass center point (Center of Gravity) of each link in its corresponding DH coordinate system. In the Stäubli technical data Sheet, Position and Orientation of each frame relative to the previous (stäubli) frame ( Āi-1 A i ) is provided. As shown in Figure 3.1, A i represents the origin of frame i attached to link i of the robot and X i and Z i depict the direction of the frame and finally G i illustrates the position of Center of Gravity of the link. Red frames shown in the figure represent the corresponding DH frames attached to the robot links.

To define the Coordinates in DH system, first we convert all the data in the base frame (Frame A 0 in figure 3.1). By having the transformation matrix between each consecutive system (T i-1 i ), transformation matrix defining frame i in base frame is: To find the transformation matrix between the base frame and our conventional DH frames, Transformation matrices defining each coordinate system in the base frame (A 0 ) are easily computed using the DH values of robot. These transformation matrices can also be computed using the Direct Kinematics of robot having joint variables all equal to zero but spacial care should be taken about the translation between the DH base frame and the one of Technical data sheet. We recall that the base frame in DH representation is coincident with frame 1 in Figure 3.1 since the value of joint 1 is equal to zero.

T 0 i = T 0 1 × T 1 2 × ... × T i-1 i (3.1)
Let D 0 i be the transformation matrix defining the DH coordinate system attached to link i in base frame. The coordinate of mass center point of each link (CG i ) described in its corresponding DH frame is determined using the following expression:

CG i = (D 0 i ) -1 × G i 0 (3.2)
By having the coordinate of CoGs defined in DH frames, one can calculate the torque applied on the each joint axis caused by the gravitational forces while the robot base has been mounted on the floor by using a recursive algorithm starting from joint 6 and by applying the equilibrium equation over each link. The value of torque imposed on join i is determined by the following expression:

   T i = T i+1 + -→ d i × m i -→ g + -→ r i × ∑ N j=i+1 m j -→ g , (i < 6) T i = -→ d i × m i -→ g , (i = 6) (3.3)
Where T i is the torque vector applied on joint i, -→ d i and -→ r i are vectors connecting the origin of frame i to the Center of Gravity of link i (CG i ) and the origin of the frame i + 1 respectively. And m i stands for the mass of link i. Finally, the resulted torque vector is projected along the axis of the joint (z i ) since we neglect the lateral bending of the transmission system:

τ i = T i . -→ z i ) -→ z i (3.4)
Figure 3.2 depicts the described parameters for computing the imposed torque on joint 2 for an arbitrary robot configuration due to the gravitational force. The output of this part of model is thus the 6 × 1 array containing the torque components projected to their corresponding joint vector:

τ sg = [τ 1 , τ 2 , τ 3 τ 4 , τ 5 , τ 6 ] T (3.5)
But however, these loads may not be the values directly applied on the transmission system of the joints. Thus for passing the self gravity originated torques directly applied of joints in VJM method, there remains another process. In some industrial robots (Stäubli TX200 included), the manufacturers implement a balancing system whose main duty is to reduce the gravity effect on a principle joint (usually joint 2). This balancing system eliminates a part of the gravitational load on its way to the gear box of the targeted joint.

In Stäubli TX200 robot arm, joint 2 possesses a balancing system for which the counter torque applied on the joint from the balancing system must be computed. Nest section is dedicated to formulate this balancing system.

Balancing System

To reduce the effect of gravity on pose errors in heavy IRs, some robot manufacturers have implemented gravity compensator systems in robot structures. These systems which use either spring suspension mechanisms (such as Stäubli RX170 surveyed by [START_REF] Olabi | Improving the accuracy of industrial robots by offline compensation of joints errors[END_REF] or pneumatic ones, as presented in [START_REF] Klimchik | Calibration of industrial robots with pneumatic gravity compensators[END_REF], are usually attached to the second joint of the targeted robots.

Stäubli TX200 robot uses a balancing system armed on joint 2 to reduce the effect of gravitational torque applied over this joint. This system consists of a spring attached to a fixed point on link 1 at one end and close to the extremity on link 2 on the other end (Figure 3.3). The spring is pre-stretched to apply a pre-load P, thus by considering that a rotation of θ causes the length of the spring expand from initial length L to final length L , the total value of the imposed fore F from the balancing system to link 2 is expressed by the following equation:

F = P c + k(L -L) (3.6)
Where k is the stiffness of the spring. As shown in Figure 3.4 the system is designed to be completely vertical while the joint 2 is at zero position.

The final length of spring L is computed using the following expression:

L = L 2 + 2r(r + L)(1 -cos θ) (3.7)
Rotating joint 2 by θ(rad) generates an angle between the direction of force F and vertical direction, shown as α in figure 3.4. The value of this angle is expressed as:

α = tan -1 (r + L) sin θ (r + L) cos θ -r (3.8)
The torque applied on joint 2 resulting from the balancing system is given by: FIGURE 3.4: Force imposition of balancing system on joint 2

T b = F × r × sin α (3.9)
Where every parameter is either a constant coming from the specification sheet of the manufacturer or a function of the joint variable θ 2 . The portion of the gravitational torque imposed on joint 2 is then computed as follows:

τ sg 2 = τ 2 -T b (3.10)
To observe the effect of this balancing system over the performance of the robot, we simulated a movement of the robot in which joint 2 rotates from -90(Deg) to +90(Deg) while having other joints at their zero position constantly with no external load applied on the end-effector. Fig. 3.5 depicts the difference of applied torque on joint 2 for two cases of unarmed and armed joint with balancing system in blue and red lines respectively. As one can see balancing system massively compensates for the self-gravity effect of the robot on the applied torque on joint 2 which is comparable with the one generated in Olabi et al., 2012. Figure 3.5 shows that the applied torque on joint 2 without any external force can rise up to about 3000(N.m) without using the balancing system (blue curve) while the additional system (red dash) has decreased this value to about 400(N.m). This reduction of about 87% proves the efficiency of the balancing system. But however, the self-gravity error source is the main contributor to robot TCP pose error and thus requires further improvements by the elastostatic model.

In this section we implemented the effect of balancing system on reduction of self gravity toque induced on joint 2. Thus at the end of this part, the actual torque values which are directly imposed and tolerated by the robot joints resulted from the link and actuator weights are explicitly computed. In the next part we add the external force effect to the model to achieve the net load values applied on joints for an arbitrary given condition for a robot task.

Effect of External Couple

To take into account the effect of external couple resulted from machining forces, one should convert the exerted couple into the resulted torques applied on robot joints. Although these coupes result in both force and torques on the joints and links of the robot, but due to the assumption of VJM method, only the projection of torques in the joint axes directions are required for the model (Figure 3.6).

Figure 3.6 illustrates the definition of τ e i where subscript e stands for external (which states the origin of the imposed torque) and i refers to the corresponding joint number. According to [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF], the following expression is used to compute τ e as a function of external force:

τ e = J T × F (3.11)
Where J is the Jacobian matrix of robot arm and τ e and F are the 6 × 1 arrays containing the imposed torque of robot joints and external force respectively with the following arrangement:

τ e =         τ e 1 τ e 2 τ e 3 τ e 4 τ e 5 τ e 6         , F =         F e x F e y F e z T e x T e y T e z         (3.12)
By having the torques resulted from the external couple computed, the total torque imposed on the joints (τ t ) is trivially the summation of both self gravity and external couple effect:

τ t = τ e + τ g (3.13)
The computation of the torques tolerated by joints were necessary to simulate the compliance effect. next section describes how one can achieve the deflection of robot structure and consequently TCP displacement due to compliance effect by using VJM method.

TCP pose error computation

As we assume linear torsional springs to model robot joints, the deviation of each joint is only a function of its corresponding projected torque described in the previous part. Thus the relationship between the joint torques and the angular deviations can be expressed by the following equation:

        τ t 1 τ t 2 τ t 3 τ t 4 τ t 5 τ t 6         τ t =         k 1 0 0 0 0 0 0 k 2 0 0 0 0 0 0 k 3 0 0 0 0 0 0 k 4 0 0 0 0 0 0 k 5 0 0 0 0 0 0 k 6         [K] ×         δθ 1 δθ 2 δθ 3 δθ 4 δθ 5 δθ 6         ∆Θ (3.14)

Modeling and Compensation of Compliance Behaviour of Robot

29

Where in angular stiffness matrix [K], elements k i represent the angular stiffness of ith robot joint and array ∆Θ, as presented in the equation, contains the value of deviations of all robot joints (δθ i ).

Thanks to the diagonal configuration of angular stiffness matrix [K], the inverse of this matrix is the inverse of each element of the diagonal. Thus by having the torque values applied on the joints, the angular deviations can be computed easily as:

∆Θ =          k -1 1 0 0 0 0 0 0 k -1 2 0 0 0 0 0 0 k -1 3 0 0 0 0 0 0 k -1 4 0 0 0 0 0 0 k -1 5 0 0 0 0 0 0 k -1 6          [K] -1 ×τ t (3.15)
By assuming that joint deviations are small, to translate the angular deviations in the joints to the displacement of the TCP, Jacobian matrix is used According to [START_REF] Khalil | Modeling, identification and control of robots[END_REF], jacobian matrix (J) relates the differential of position and orientation vector to the differential of joint angles by the following expression:

δP 3×1 ω 3×1 = J 6×6 × ∆Θ (3.16)
In which the part δP is the position displacement of the TCP (δP = [∆X, ∆Y, ∆Z] T ) and ω is an array containing rotational disorientation which however can not be used without further treatments to indicate changes of any rotational vector (For more information on this topic, reader can refer to [find a reference]).

By replacing ∆Θ from equation 3.15 into equation 3.16, the relation between TCP pose error and imposed torques on the joints is expressed as:

δP 3×1 ω 3×1 = J 6×6 × ([K] -1 6×6 × τ t ) (3.17)
In equation 3.17, the Jacobian matrix (J) is a function of robot configuration (J = J(θ i | i=1:6 )) and torque values stored in τ t array are also a function of robot configuration and external forces which are the inputs of the model. Thus for having a complete VJM model for compliance behaviour of robot, the only missing elements are the stiffness values of robot joints placed in compliance matrix ([k] -1 ). Next section thoroughly explains our methodology for determining the stiffness values of robot joints.

Joint By Joint Stiffness Identification Joint 1 Stiffness Identification

Machining processes are often executed having work piece laying over a table with horizontal surface in front of robot. Thus compliance of joint 1 substantially affects the robot performance in machining process since it is most of the time the only joint undertaking couple resulting from machining workpieces lying on horizontal planes (parallel to XY of robot base frame). To compensate for the effect of Joint 1 in our Elastostatic model, determining the stiffness value of this joint is undeniably important.

Measurement Device

To achieve the stiffness value of joint 1, we made use of Faro Arm Measurement device. In this method, robot arm was completely bent to lie over horizontal plane so that by applying a force parallel to Y axis of robot, maximum lever will cause maximum torque on axis 1. Faro Arm was mounted on a table next to the extremity of the second link in a way that three points on the extremity of link 2 of robot were conveniently achievable in the workspace of Faro. Faro arm was fixed to the considerably heavy table by a pneumatic system preventing the base of the apparatus to move by imposing a strong suction resulting from air vacuuming on the contact zone of Faro base and table surface. Aforementioned 3 points are chosen to be the head of three screws on the robot shell. The end-effector of Faro is an accurate sphere with diameter of 5(mm) enabling its TCP to place in same point by meeting the head of the screws with different orientations. Figure 3.7 depicts the details of this measurement procedure:

Force Imposition Method

To apply torque on Joint 1, a force parallel to axis Y of the base frame of robot was applied on the end-effector. End-effector used for this measurement session was a Ushaped iron plate to ensure that the connecting rope can move along a circle arc hole to self-align with the horizontal direction under tension. Force was imposed using a rope pulley system connecting an emptied oil container in one end to the endeffector of the robot on the other end. An analog force sensor was placed serial along the rope to have an accurate value of the rope tension. Though this placement will cause the rope to buckle due to the weight of sensor, but after a certain tension value this effect fades away. Force was increased and decreased by adding and removing weights to the container and the value was recorded using the analog force sensor. Figure 3.8 shows the force imposition procedure for this experiment:

Identification Procedure

Applying force at the end-effector causes torque over joint 1 and subsequently induces angular deviation of this joint. This method is based on observation of this deviation corresponding to the applied torque. Each three point are placed on a circle whose center lies on the line of axis 1 of robot. Thus by assuming complete rigidity of link 2, due to angular deviation of joint 1, these points are displacing on the circumference of their corresponding circle.The By assuming small deviations on joint 1, the following relationship converts the displacement of measurement i (∆X i ) to the angular deviation of joint 1 (∆θ 1 ):

∆θ 1 = ∆X i R i (3.18)
Where R i is the radius of the circle attributed to point i. Thus to measure the angular deviation, it is necessary to determine the radius of the circles for each three measuring spots. Joint 1 was rotated about 20 [Deg] with increments of 1[Deg] in unloaded situation to achieve a set of points on the circumference of circles generated by different positions of this three points. Figure 3.9 illustrates the methodology of radius determination for the three measurement points. These data are used to fit the best circles using Least Square method. Table 3.1 illustrates the resulted radius.

Due to the fact that only the displacement of each point from its initial position is of importance, The gathered points are recorded in the working frame of Faro Arm which was constructed arbitrarily and with no actual relation to the one of robot. Different values of force were applied on the robot while recording the corresponding displacement of each 3 measuring points. Using radius determined with the fitted circles and a lever of 2200(mm) for torque imposition obtaind from robot's direct geometry model for the given configuration plus an additional 6(mm) of the designed end-effector width, figure 3.10 shows the resulted pairs of angle deviation and applied torque for each set of loading and each three points of measurement. It is worth mentioning that the mean/max value of the displacements (Data related to all three points combined) are 0.39/0.80(mm). Compared to the uncertainty of the Faro Arm which is claimed to be 0.05(mm) by the manufacturer, measured data can be fully trusted.

By assuming linear stiffness model for joint 1, the stiffness of joint 1 is the slope of the fitted first order polynomial on the (τ i , ∆θ 1 i ) depicted in Figure 3.10. Following expression is used for least square fitting:

τ = k 1 × ∆θ 1 ⇒ k 1 = Σ i τ i * ∆θ 1 i Σ i ∆θ 2 1 i (3.19)
The resulted stiffness for joint 1 is then:

K 1 = 3.0929 × 10 6 ( N.m rad ) (3.20)

Stiffness Identification of Joint 2

Stiffness of joint 2 is also estimated using an approach close to the one of Joint 1. But measurement device is a Laser Tracker to accurately measure the displacement of reflector. Exact value of the force was then available from the gauge of the force sensor. Figure 3.11 illustrates the overall configuration of the robot arm and force imposition system for this test. The two important radius R 1 and R 2 ; namely Actuation and Observation radius are shown in Figure 3.11. Robot arm was stretched to maximize the torque applied on axis 2 but due the visibility of the reflector which was attached on the bottom side of link 2, having enough of distance to hang the bucket and other attachments and also to keep the last link horizontal to ensure a vertical force applied on the end-effector, the value of joint variables were decided to take the following values (values are in Degree ):

θ 1 = 0 θ 2 = 85 θ 3 = -10 θ 4 = 0 θ 5 = -15 θ 6 = 180
We use the nominal values of the robot to compute for R 1 . To identify the radius of R 2 we used the same approach as the one of joint 1. Following comes the values of these properties: Translating the row data into angle deflection and applied torque, we achieve the graph presented in Figure 3.12 for stiffness of joint 2. Applying least square method on the measurement data, the stiffness value of joint 2 becomes:

R 1 (mm) R 2 (mm)
k 2 = 4.1988 × 10 6 (N.m) (3.21)
Compliance Effect of Balancing system on Joint 2

In this part we survey the effect of the balancing system on the compliance behavior of joint 2. Though we have modeled all the joints of the robot with a linear rotational spring, observations over the non-linear effect of the balancing system is of great importance. Several questions may arise at this point: How large can be this nonlinear effect in terms of portion? Is it dominant over the linear part of the stiffness? consequently, within the precision of our application, should we model this nonlinear effect in our elastostatic model? Considering the fact that we have measured the stiffness of joint 2 in a certain point in its working interval (θ 2 = 85[deg]), how has the experimented stiffness value been effected by the balancing system? In this part we are aiming to answer these questions.

As described previously in section 3.2.3, the balancing system attached to joint 2 is composed of a linear spring attached in one end to link 1 and to the extremity of link 2 on the other end. Thus due to this configuration, one can conclude that the stiffness effect of this system can be modeled as a parallel torsional spring added to the stiffness of joint 2 k 2 . The torsional spring modeling the stiffness effect of balancing system is termed as k b . The modeled configuration of stiffness of joint 2 is illustrated in Figure 3.13. according to equations 3.7, 3.8 and 3.9, the torque applied on joint 2 form the balancing system (T b ) is a function of joint variable θ 2 and so is the rotational spring attributed to this system. Thus the total stiffness value of this joint is the summation of the two springs:

k T2 (θ 2 ) = k 2 + k b (θ 2 ) (3.22)
Where K T2 (θ 2 ) is the total stiffness value of joint 2, K 2 is the linear portion of the stiffness presenting the equivalent stiffness of transmission system and link rigidity and K b (θ 2 ) represents the nonlinear stiffness of balancing system.

Nonlinear stiffness of the balancing system can be calculated by the derivative of the torque expression by the joint value:

k b (θ 2 ) = ∂T b ∂θ 2 (3.23)
Using equations 3.7, 3.8 and 3.9, the expression of the stiffness of balancing system as a function of joint variable θ 2 turns to be:

k b (θ 2 ) = r×cos(θ 2 )×(P c -k(L- √ L 2 -2r(L+r)(cos(θ 2 )-1)))×(L+r) √ L 2 -2r(L+r)(cos(θ 2 )-1) + r 2 sin 2 (θ 2 )×(P c -k(L- √ L 2 -2r(L+r)(cos(θ 2 )-1)))×(L+r) 2 (L 2 -2r(L+r)(cos(θ 2 )-1)) 3/2 + k×r 2 ×sin 2 (θ 2 )×(L+r) 2 L 2 -2r(L+r)(cos(θ 2 )-1) (3.24)
Where all the parameters apart form θ 2 are constants and defined in equations 3.7 and 3.8. Figure 3.14, shows the value of balancing system's stiffness in an interval of -90

[deg] to 90[deg] of joint 2.
According to the stiffness evolution of balancing system shown in Figure 3.14, the equivalent torsional stiffness of this system (k b ) is highly dependent on the joint variable θ 2 . But however, the maximum value achieved by this parameter is less than 2200(N.m) which compared to the total joint stiffness (k T2 ) determined by experiments in section 3.2.5 is about 0.05%. This portion is even less for configuration of robot used to experimentally evaluate stiffness of joint 2 (θ 2 = 85[deg]).

Due to the negligible effect of torsional stiffness of balancing system on total stiffness of joint 2, its equivalent torsional spring k b is eliminated from the stiffness FIGURE 3.14: Evolution of the non linear stiffness of robot balancing system attached to joint 2 of joint 2. I.e:

k T2 ≈ k 2 (3.25)

Stiffness Identification of Joint 3

To identify the stiffness of joint 3, we mounted two reflectors on link 2 and 3. The reason is to eliminate the effect of compliance of joint 2 from the results so that the deflection of joint 3 will be more observable. Force imposition system for this step is the same as the one used for previous joints. To use the maximum available length to apply torque on joint 3 and simultaneously reduce the lever over joint 2, we chose the following configuration for the robot arm (values are defined in Degree):

θ 1 = 0 θ 2 = -30 θ 3 = 120 θ 4 = 0 θ 5 = 0 θ 6 = 180
Figure 3.15 shows the placement of the stands of reflectors (blue circles) and the overall configuration of the robot and force imposition system for this step.

As shown in Figure 3.15, the coordinates of the reflectors namely points P 1 and P 2 are gathered for each loading force. In this part, according to the assumptions of VJM model, since there is no torque applied on joints 1 and 4 (deviations of joints 5 and 6 do not affect the placement of P 1 and P 2 ), the displacement of the both points lie on planes parallel to XZ plane of robot frame. Hence for generating equation in further steps, these points are assumed to be on the XZ plane of robot frame. Figure 3.15 shows the normal view to this plane and the modeled schematic for robot structure deformation.

To identify the stiffness of joint 3, we need to extract the value of ∆θ 3 shown in Figure 3.15 for each loading condition from the gathered data. Before explaining the computational steps for achieving the foregoing parameter, we define other parameters shown in the figure: L 2 is the length of link 2. we use the nominal value of this parameter which is 950(mm). L 3 is the distance between P 1 and origin of joint 3 (shown as J 3 in Figure). This parameter is determined using the circle construction method by measuring the position of reflector 1 while having axis 3 rotating. The value of this parameter is L 3 = 803.7(mm). L 1 is the distance between P 2 and origin FIGURE 3.15: Configuration of the robot and force imposition system for joint 3 stiffness identification process of joint 2 and is equal to L 1 = 937.4(mm). ∆θ 2 and ∆θ 3 are the deflections caused by compliance of axes 2 and 3 respectively. D 1 , D 2 and d 3 are the distance between the initial and final position of points P 1 , P 2 and J 3 respectively.

To extract the value of ∆θ 3 , we first make use the value of D 2 captured from the LT, to compute ∆θ 2 :

∆θ 2 = D 2 L 1 (3.26)
d 3 and subsequently the final position of J 3 (J 3 ) are computed having the value of ∆θ 2 . D 1 is computed using the recorded coordinates from Laser Tracker. The final position of point P 1 (P 1 ) is the intersection of the following two circles:

C 1 : (x -x P 1 ) 2 + (y -y p 1 ) 2 -(D 1 ) 2 = 0 C 2 : (x -x J 3 ) 2 + (y -y J 3 ) 2 -L 2 3 = 0 (3.27)
The mentioned set of equations are solved using fsolve function in MATLAB software. Finally, the expression defining ∆θ 3 is:

∆θ 3 = atan( y J 3 -y P 1 x J 3 -x P 1 ) -∆θ 2 (3.28)
As in this step of measurement we have a set of data for the deflection and applied torque on joint 2, we can again estimate the stiffness of this joint and compare the result to the one of previous section. Applying Least Square method releases the following values for the stiffness of two joints (measured data and the fitted model for this step are shown in figure 3.16), :

k 3 = 1.0159 × 10 6 ( N.m rad ), k 2 n = 2.9642 × 10 6 ( N.m rad ) (3.29)
We can see a difference of about 30% between the stiffness of joint 2 resulted from the current and the previous test. We must mention that the displacement of reflector 2 was at maximum 0.1(mm) (average of 0.03(mm)). Considering the uncertainty of the LT in a distance of about 5(m), the resulted value for the stiffness of joint 2 from this measurement is NOT reliable. But the exact fact that the placement of reflector 2 is very close to "stationary" increases the reliability of the resulted value of stiffness of joint 3.

Stiffness Verification of joint 5

To identify the stiffness of this joint we used the Faro Arm instrument but with a new end-effector designed especially for this joint(Figure 3.17).

FIGURE 3.17: Robot Configuration and introduced frame attached to link 4 for joint 5 stiffness identification

As shown in figure 3.17, we used a container-weight system to apply torque on joint 5. This system was connected to an analog force sensor enabling us to have the exact value of the force applied.

The force applied on the end-effector will cause the end effector to displace and this effect is not only due to the deformation of joint 5. Thus recording the absolute displacement of end-effector is not an accurate solution to verify the stiffness of this joint.

To tackle this problem, we used a new strategy to eliminate the effect of compliance of joints prior to joint 5. For each loading condition, a new frame attached to link 4 of robot is reconstructed and the coordinates of predefined points on the end-effector is recorded in the updated frame. The deviation of the joint 5 is equal to the angle between the initial and final normal of the plane ∆θ:

∆θ 5 = arccos | -→ n . -→ n | | -→ n | * |n | (3.30)
In which:

-

→ n = -----→ (P 1 -P 2 ) × -----→ (P 1 -P 3 ) (3.31)
Three different configuration of robot were chosen to execute the measurements. Using Least Square method to calculate the stiffness of the joint 5 reveals the following results for each configuration: Figure 3.20 shows the experimental data and fitted model of these measurements. Data show that stiffness of this joint is very sensitive to the configuration of the robot and also its behavior is highly nonlinear for low values of the torques. But one important point should be noticed that the measurement system (Faro Arm) has 0.1 (mm) of uncertainty which is not considered accurate specially for low values of torque that can not cause an effect of more than 0.1 (mm) in displacement of points. In our elastostatic model, the effect of axes 4 and 6 is neglected due to the short lever of TCP on these axes. This is of evidently a function of end effector geometry and robot configuration, but however the normal distance of TCP from axis 6 for example, is usually within a range of 0 to 50(cm). This lever will cause a torque of 100(N.m) for a perpendicular force of 200(N). This torque would cause a deviation of maximum 0.03[deg] on joint 5 which is by far weakest joint among the four examined joints.

Through out this section, we determined the stiffness of robot joints which are main contributors to the compliance effect of robot. According to equation 3.17, developing the elastostatic model based on VJM method is accomplished by determining the elements of stiffness matrix. To verify the efficiency of this model and determined joint stiffness values, the strategy is to compare its performance with gravity compensation option implemented in robot controller. To go further on this strategy, firstly a brief description of feature is given and later comes a detailed explanation on the verification procedure.

Stäubli TX200 Gravity Compensation Option

Staäubli TX200 robot arm is capable of compensating the gravitational effects on TCP displacements. This feature can counter the effect of link masses plus the weight of additional device attached to the robot flange (for the latter the operator has to define the mass and coordinates of center of gravity in the user interface of robot controller). The details of the of the method which this feature uses for to compensate gravity effects (VJM, MSA or FEA) is not in hand but however through the user interface of this option, we concluded that only the stiffness of joints 2 and 3 are utilized for this feature. This seems reasonable since joints 1, 4 and 6 do not undergo high values of gravitational loads and contribution of joint 5 is considerably lower compared to joints 2 and 3 due to shorter lever. The compliance values (inverse of stiffness values) of joints 2 and 3 defined by Stäubli used in the feature were also achievable through the interface and are defined as:

c 2 = 2.88 × 10 -5 ( [deg] N.m ) c 3 = 7.23 × 10 -5 ( [deg] N.m ) (3.32)
To comply with our notation, the compliance values must be converted into stiffness values having N.m/rad dimension.

k( N.m rad ) = 1 c( • D N.m ) × 180 π ⇒ k 2 cont = 1.989 × 10 6 ( N.m rad ) k 3 cont = 7.925 × 10 5 ( N.m rad ) (3.33)
A glance of the stiffness values show a considerable difference between the experimentally achieved values with the ones of Stäubli controller. Thus two sets of stiffness values are in hand with different sources. Hereafter, the stiffness values achieved by experiments are termed as Experimental set and the values defined by Stäubli company are referred to as Stäubli set. Table 3.3 summarizes the two different sets of stiffness values.

The infinity sign in the value sets in table 3.3 show that the corresponding joint will be considered an ideal joint. An ideal joint is a joint which shows no deflection to any external torque imposition or in other words it has a stiffness value equal to infinity.

Value Set

k 1 k 2 k 3 k 4 k 5 k 6 Exp. 3.0929 × 10 6 4.1988 × 10 6 1.0159 × 10 6 ∞ 2.43 × 10 5 ∞ Stäubli ∞ 1.989 × 10 6 7.925 × 10 5 ∞ ∞ ∞ TABLE 3.3: Joint Stiffness Values used in Elastostatic Model (Defined in N.m/rad)

Elastostatic model verification

The elastostatic model explained earlier in this chapter, mainly contains two different parts. One part is the mathematical expressions such as equations 3.3, 3.7 and 3.17 which are implemented in Matlab software and the second part is the stiffness values which can be either the experimental set of Stäubli set of values.

To verify the performance of the elastostatic model with determined stiffness values, from the stated different sources, an experiment has been designed and conducted. This verification strategy is based on comparing the performance of the following three different approaches for compensating a path resembling a practical mission of the robot: Comparing the efficiency of the compensated path using approach 1 with the one of 2 would reveal which set of stiffness values are more accurate. Comparing the results achieved by approach 2 with approach 3 would proof whether the generated mathematical model for the elastostatic model is well implemented. Finally, we will conclude that by the materials in hand, what would be the percentage of the compliance errors of robot behaviour we can compensate for.

Compensation Algorithm

The explanation of the compensation method is brought through the strategy of the elastostatic model examination. But as a brief statement, the compensation strategy is based on computing the compliance resulted error (∆P in equation 3.17) and consequently program the robot to add the error value with an opposite direction to its desired path.

To go further with the verification procedure, it is noteworthy to add that to be able to compare the performance of the elastostatic model with the load compensation option, we make use of ability of load compensation in compensating for external mass attached to robot flange. In other words, comparing the performance of elastostatic model with load compensation option in terms of self gravity is straight forward. But however, the performance of generated elastostatic model in compensating for external load can only be verified and compared with the load compensation option when the external load is in direction of gravity. Only in his case, this load can be defined as an external mass attached to the robot flange and thus the load compensation option of robot controller can counter the effect. This rectangle is shown in blue as the desired path (X d ) in Figure 3.21. TCP of robot is simulated to track the circumference of this rectangle with an external force of 500(N) in gravitational direction. Due to the external force and self gravity, elastostatic model estimates the black path for this movement. This path is termed the Actual path (X a ) as it will be the practical path carried out by the robot according to the model. Each point on the desired path is related to its corresponding point on actual path by an error parameter Er which is calculated by the following expression:

Er = X d -X a (3.34)
To compensate this error coming from both self gravity and external load, we will add up the error to the coordinates of desired path to have a trajectory for which if the robot is programmed to track, the actual path will lie very close to the desired path. This path is called the Compensated Path (X c ) shown in red in figure 3.21 and is obtained by the following expression:

X c = X c + Er (3.35)
In next step we describe the tests prepared for elastostatic model verification. Several path are generated and examined based on two different stiffness values and different compensation models.

Test Protocols

To verify performance of different approaches stated in section 3.2.7, ten different tests were designed. For tests which utilize the generated elastostatic model to compensate the path, Matlab software has been used to generate the compensated path and tests which adopt the load compensation feature of robot controller, only a mass equivalent to the corresponding load is defined in the robot controller. Tables 3.4 and3.5 show the details of the designed tests including their corresponding compensation method. Two different load parameters are used in tables 3.4 and 3.5. Actual load is the one physically applied on the end-effector of robot and compensation load, is the one that is used in either models to compute for the compensated path.

The predefined rectangle is used as the desired path. Robot was programmed to follow the path in a quasi-static situation (2% of the nominal speed value) to avoid any kind of dynamic effects and position of the TCP was recorded using API TIII Laser Tracker during robot motion. Figure 3.22 depicts the test configuration for this experiment.

As shown in the Figure 3.22, the actual load was applied on the robot end effector using the method described for stiffness identification of joint 3 to ensure that the external force remains always in the gravitational direction.

The resulted path from the ten different tests are brought in Figure 3.23.

To compare different path results, we assign two values to each path namely D 1 and D 2 shown in figure 3.21 to quantify the efficiency of each. These values are the errors corresponding to the further and closer sides of the rectangle to the robot frame respectively. Values of the defined parameters for each resulted path for the ten tests are presented in the table 3.6.

Discussion and Conclusions

• According to the descriptions of tests 7 and 9, they correspond to the compensation approaches 2 and 3 defined in section 3.2.7. As the two paths resulted from the mentioned tests are practically identical, we can conclude that the mathematical implementations of elastostatic model, self gravity and balancing system attached to second joint are functioning properly.

• The comparison of tests 7 and 9 also reveals that with a high level of certainty the load compensation feature of Stäubli uses VJM model for compliance error.

If the feature uses other methods, it is clear that VJM model is as efficient as the model.

• The path resulted from tests 7 and 9 show almost no error in the front side and 0.3(mm) in the back side of the rectangle while the path resulted from test 6 holds errors of 0.9(mm) and 0.5(mm) for the both mentioned sides. By adding the fact that test 6 represents approach 1, we can draw the conclusion that the stiffness values determined by Stäubli is more accurate compered to the ones resulted from experiments. The Following points may explain this outcome:

-The stiffness identification method used by Stäubli R&D group consists of dismounting the transmission system from the robot, mounting it on special fixture and applying the stiffness procedure. This enables one to bypass the link deformation effect and its impact on the stiffness identification.

-For a static stiffness identification procedures, factors such as joint frictions and hysteresis effects are neglected. these factors may affect the results as reported by [START_REF] Ruderman | Modeling and identification of elastic robot joints with hysteresis and backlash[END_REF] .

• In the middle point of the back side, a sudden fall is vivid in all the test paths. This sudden drop of about 0.2(mm) in z direction is the effect of backlash. Figure 3.24 shows a closer view of this point with the corresponding point in the Joint value trajectory: FIGURE 3.24: Effect of backlash in the experimental path Figure 3.24 shows that the created step is actually a result of accumulation of backlash of the three following joints: Joints 2,3 and 5 as they are undergoing a change in their direction of rotation.

The final conclusion of this verification test is basically the update of the stiffness values of robot joints. We are certain that the elastostatic model is functioning properly but regarding the joint stiffness values we remain with the experimental values for joints 1 and 5 and replace the experimental values for joint 2 and 3 with the Säubli values. One should note that during the verification test, the effect of joint 1 on the compliance behaviour was eliminated due to the loading configuration.

In light of above mentioned points, hereafter, any mention of elastostatic model through out this thesis report, addresses the set of stiffness values expressed in table 3.7 for simulations and compensations.

k 1 k 2 k 3 k 4 k 5 k 6 3.0929 × 10 6 1.989 × 10 6 7.925 × 10 5 ∞ 2.43 × 10 5 ∞ TABLE 3.7: Joint Stiffness Values used in Elastostatic Model (Defined in N.m/rad)
In next section we will introduce a novel way for determining the stiffness values of robot joint on real working conditions. We should note in advance that results presented in next section serve only as an introduction to this model and won't be used for the elastostatic model due to lack of experimental validations for the elastostatic model based on these results.

The procedure used for robot joint stiffness identification in section 3.2.5 was based on static behaviour for only one robot configuration for each joint (apart from joint 5 for which the identification took place for 3 different configurations). According to [START_REF] Jang | Calibration of geometric and non-geometric errors of an industrial robot[END_REF] joint stiffness values of joints might differ with the robot configuration.

Active Stiffness Identification

In this part we make use of the gathered data for rectangle test presented in previous section to build up a procedure for which a wide zone of robot workspace imposes its impact on stiffness values identification. Other noteworthy difference of this method compered to the previous one is the active condition of robot. I.e. for static condition, the effect of friction and joint breaks may affect the results but however in moving condition of robot the effect of these phenomenon will decrease. The name Active stiffness identification adopted for this method expresses the moving condition of robot rather than its completely static condition.

Identification Procedure

In the active stiffness identification procedure which is explained in this section, we make use of the data of the paths resulted from tests number 1 and 2 presented in Figure 3.23. We recall that these two paths are resulted from unloaded and loaded situation of robot arm respectively with no compensation method applied. Thus according to the linear VJM based model adopted in this work, the different between the two path is only resulted from the external 500(N) in gravitational direction on the end effector of robot. Thus by subtracting the two resulted paths, the rest of error sources such as self gravity, geometrical errors, frame construction procedure and nonlinear behavior of reducers as well as backlash will be eliminated.

To subtract the two resulted paths, each side of the rectangle were discretized by each millimeter dividing the whole rectangle into 4400 points, 1000 points on each lateral sides (S 1 & S 3 in Figure 3.24) and 1200 points on each longitudinal sides (S 2 & S 4 ).

To eliminate the effect of transmission joint errors, each side Longitudinal sides were fitted with order 4 polynomials using LS method and order 3 was chosen for the lateral sides as shown in figure 3.25. The difference between z component of the corresponding points on each path is termed as E standing for error resulted from the applied force. Next step of this identification method is to find the stiffness values of joints which best simulate the deviation in Z axis resulted from a negative vertical force of 500(N) to be as close to the one resulted from the measured and fitted paths. In other words the problem changes to an optimization problem to minimize the following function:

Er = N ∑ i=1 E(i) -Es(i) 2 N (3.36)
Where E(i) is aforementioned error function generated from the measurement results for point i, Es(i) is the error estimated from elastostatic model for the same point and N stands for the number of points which is in this case 4400.

Since the effect of self-gravity and balancing system of Joint 2 are eliminated due to the relative difference of the two paths, the elastostatic model will be simplified only to the external force effect. To compute the estimated error E using elastostatic model, the parameter τ t in equation 3.17 must change to τ e defined in equation 3.11. Thus the following expression is used for estimating the displacement error resulted from external load effect:

δP t ω t test 2 - δP g ω g test 1 = J × ([K] -1 × (τ g + τ e )) -J × ([K] -1 × τ g ) = J × ([K] -1 × τ e ) τ e =J T ×F -----→ δP t -δP g ω t -ω g = (J × [K] -1 × J T ) × F (3.37)
Thus for computing the error parameter E which is the third component of array δP t -δP g , one only requires the following information: Jacobian matrix J for the given configuration, compliance matrix [K] -1 and applied force F. Out of which, only compliance matrix is assumed unknown.

To simplify the problem, we will only take into account the effect of joints which are considered to have the largest portion of the displacement of TCP. To proceed further, joint torques are the best parameters to preliminary detect the joints with the most effect. Figure 3.26 depicts the torque resulted from the force on each joint while passing the rectangular path on the right and error value E for the entire rectangle circumference on the left. Due to the static measurements on joint 5 maximum deviation of q 5 corresponding to this torque value is around 8 × 10 -4 (rad); having a lever of 200(mm), maximum 50 Chapter 3. Identification And Compensation Of Non Geometrical Robot Errors deviation resulting from this joint is about 0.15(mm) which is in the order of nonlinear errors of reducers and can be neglected. Maximum torque applied on Joint 4 is also about 50(N.m). This joint will also be canceled out for the further steps for the sake of simplicity.

Further step for active stiffness identification is to find the best stiffness values of joints 2 & 3. To observe the effect of the two variables (K 2 &K 3 ) on the estimation of the deviation resulting from the elastostatic model, for a reasonable interval, error parameter (Equation 3.36) has been computed. Figure 3.27 shows the mentioned parameter as function of K 2 and K 3 . In comparison between Stäubli values and optimal set of stiffness values resulted from the active stiffness identification identification, we have noted that the estimated error values depicted in Figure 3.28 differ up to 0.2(mm) at maximum for a condition where TCP experiences up to 1.1(mm). Even though the optimal values predicted more accurate error values (since they were found using exact same error values), we stay with Stäubli since this difference is in the order of cyclic errors observed along the path.

In next section we have dedicated our effort into mathematically model the cyclic errors which as previously stated in chapter 2, are due to the imperfections of transmission systems implemented in robot joints.

Cyclic Joint Transmission System Errors

There are different approaches to overcome joint errors. Robot performance can be improved by using either direct-drive motors (as in some SCARA robots) or high precision gear trains characterized by very small backlash, high efficiency, low inertia, low friction, high stiffness and low weight; or by placing high-resolution encoders at the output of the gear trains. However, these solutions cannot be applied simultaneously and would raise the manufacturing cost of an industrial robot considerably [START_REF] Slamani | Characterization and experimental evaluation of gear transmission errors in an industrial robot[END_REF] Conventional way to measure the transmission error of a transmission system is by evaluating single flank gear meshing. This method has been developed by Ming in [START_REF] Ming | Measurement of transmission error including backlash in angle transmission mechanisms for mechatronic systems[END_REF] to be able to measure for backlash value. But their method requires a direct access to the motor drivers and attaching external encoders. In spite of having all these possibilities they have shown that the value of backlash is a function of meshing point in the gear reducer.

An experimental study done by Taek [START_REF] Oh | Influence of the joint angular characteristics on the accuracy of industrial robots[END_REF] shows that in both broadly used reducer types in robotic applications (harmonic and Cycloidal drives) two different errors appear. The so called Major one is the result of misalignment between the center of the gear assembly and input shaft in the reduction gear transmission. The second error type, smaller in amplitude, is a wave propagation superimposed on the major error and it is due to the gear pitch error, also termed as Teeth-to-teeth error. In his study, he attached encoders to the joints. For vertical joints, an electrical lever and for the horizontal ones he used an autocollimator system to generate reference which shows this method requires a complete authority and ability to dismantle the links of robot which is not the case for many applications.

Being pioneer in using a Laser Tracker for transducer error studies, in [START_REF] Slamani | Characterization and experimental evaluation of gear transmission errors in an industrial robot[END_REF] have shown that bidirectional repeatability is strongly affected by backlash. Their results show that backlash is highly dependent on the robot configuration in view of wear effect As they observed the maximal value of backlash at home position of the robot. By programming the robot to pass a 1 (m) straight line in which in the second half of the line joints 2 and 3 following the exact mirror behavior of the first half, they showed that the overall error shape remains the same thus these errors are only due to the gear errors on the mentioned joints.

In LISPEN laboratory, the experience of working on identification and compensation of nonlinear gear transmission behavior dates back to 2012 starting from the work of [START_REF] Olabi | Improving the accuracy of industrial robots by offline compensation of joints errors[END_REF] Although the authors did not present any practical results in the mentioned paper, but the final result of the compensation of the nonlinear defeats of transmission system was presented in Olabi, 2011. There the author identified the gear defeats of joints 2 and 3 of a Stäulbi RX 170BH robot arm and compensated the articular trajectory related to a straight line in Cartesian coordinate system. The results are depicted in Figure 3.29. Their proposed method was based on constantly measuring the position of a certain point attached to the corresponding link of the desired joint. Then the robot was programmed to sweep the working interval of the desired joint with increments of 0.05[deg], stop for 3 seconds on each pose so that the laser tracker can record the positioning data of the reflector (Figure 3.30). In spite of the great performance of this method in terms of limiting the fluctuations of the TCP induced by gear errors within the repeatability of the robot, the approach was time consuming. Each joint requires about 2000 points causing a downtime of about 8 hours (almost a working shift) of the robot to identify the first three joints. In this part we propose a new method to avoid the time consumption problem. In this method we have proposed a new approach to define the trajectory of the robot for transmission error identification.The error of all the desired joints can be identified at the same time, which can substantially reduce the downtime period of the robot. We target the acceptable range of the gear-caused fluctuations to be within the repeatability of the robot manipulator.

Cyclic transmission Identification Methodology

The proposed method consists of the following steps: All the steps are described in detail along with a practical case study. In this experiment we have identified the gear errors of a robot manipulator for the three first joints using an optimal trajectory and subsequently compensated these errors for a straight line trajectory in a demanding zone of the robot workspace.

Case Study

The proposed method was applied on a Stäubli TX90 robot in LISPEN laboratory with a floor-mounted configuration. It is mounted on a considerably rigid stand to avoid any base deformation due to the gravitational reactions for different configuration of the robot. To measure the TCP coordinate values we used an API laser tracker providing 20(µm) of accuracy for the test configuration. Figure 3.32 shows the configuration of the experimental setup. To convert the acquired data from the laser tracker into the robot base frame, the base frame was identified by rotating joints 1 and 2 as explained by Besset, Olabi, and Gibaru, 2016.

Optimal Trajectory

In general case of 6 axis robots, due to the lever effect and higher gear ratios, the effect of first three joints are assumed to be dominant in the positioning error of the robot. Since searching for the best trajectory in Cartesian space leads eventually to a varying increments of discretized points in articulation space, the proposed trajectory is resulted from discretion of a certain interval of joint values in joint space. The trajectory is composed of a number of poses in which the robot has to stop and the coordinate of the end-effector is measured. Equation 3.38 expresses the general format of the joint variables of the trajectory:

θ i = [θ 1 (0) + i * ∆θ 1 , θ 2 (0) + i * ∆θ 2 , θ 3 (0) + i * ∆θ 3 , 0, 0, 0] (3.38)
Where:

θ 1 (0) = -50(deg) ∆θ 1 = 0.05 * 2 9 (deg) θ 2 (0) = 0(deg) ∆θ 2 = 0.05(deg) θ 3 (0) = 90(deg) ∆θ 3 = -0.05(deg)

Error Modeling

Deviation of the actual joint values from the theoretical values for joint j (∆θ j ) originate from different sources such as joint compliance and imperfect DH values used in the inverse kinematic of the robot along with the periodic gear-transmission fluctuations. To classify the nonlinear behavior of the transmission system apart from the other error sources, we use the following notation:

∆θ j = ∆θ j n + ∆θ j or (3.39)
Where subscripts n and or stand for the error part induced from the nonlinear behavior of transmission system and accumulation of the other sources of error respectively.

To distinguish the latter from the other sources, a 3rd order polynomial is fitted to the residual error of each joint and subsequently subtracted. Due to our assumption, this polynomial line represents the effect of self-gravity and geometrical imperfections. Figure 3.34 illustrates the total error value of the three first joints of the robot 56 Chapter 3. Identification And Compensation Of Non Geometrical Robot Errors along with the polynomial curves representing the so-called or part of the error resulted form optimal trajectory. FIGURE 3.34: Total joint error values (∆θ j ) in blue along with ∆θ j or part of error in Red for joint j (j : 1, 2, 3)

To cancel out the effect of other error sources, we subtract the fitted polynomial from the total error values. The residual part is what we claim that is substantially caused by the nonlinear behavior of transmission gears. Due to the periodic nature of this error, we have chosen to fit summation of sinus functions to express the error for a given joint j:

∆θ nj = N ∑ k=1 a k sin(b k θ j + c k ) (3.40)
The fitting process was carried out using Curve Fitting toolbox in Matlab. To choose The value of N (Number of sinus functions for each joint), we started from 1 and marched forward up to the first number whose amplitude of Nth part is less than 10 -5 (rad). N = 7 was the first number satisfying this condition for all the joints. Figure 3.35 shows the gear errors with the fitted sinus functions resulted from curve fitting.

Error properties (Amplitude, harmonies and their phase differences with the corresponding joint ) of joint 1,2 and 3 are presented in tables 3.8, 3.9 and 3.10 respectively.

Off-line correction procedure

So far we discussed about the identification and quantification of gear induced errors. Hereafter, the achieved mathematical expressions representing gear errors are termed Error Model. For a desired trajectory, primarily, the desired Cartesian properties should be converted into the articular space using inverse kinematics of the robot and only in this space one can estimate and subsequently eliminate the nonlinear gear induced errors of the trajectory by this proposed approach. Elimination of the error is only the subtraction of the estimated error form the desired joint variables:

θ j Compensated = θ Desired -∆θ estimated (3.41)
Figure 3.36 illustrates the compensation procedure where IK and FK stand for inverse and forward kinematics respectively. Using forward kinematics to achieve the Cartesian properties of the compensated trajectory is an optional step depending on the the way that robot is commanded. Mentioned fact reasons presenting of this step with with dashed lines in Figure 3.36.

Results and Verification

Validity and performance of the proposed method was examined by executing a simple test. In this test, joint values corresponding to a 25(cm) straight line were compensated with the achieved error models of the three first joints. To represent a general condition, this test was designed with the following features:

• The trajectory of the line was designed in the Cartesian coordinate

• The position of the line was chosen to have least intervention with the optimal trajectory to verify if the joint gear error signatures remain constant

• Wrist joints (4,5 and 6) were commanded to move constantly (By keeping the orientation constant for the flange frame in Cartesian frame)

To satisfy all the mentioned criterion, we programmed the robot to follow a straight line of 25(cm) with the starting point of X s = [500, 50, 0](mm) and ending point of X e = [750, 50, 0](mm). The path was generated point-by-point with increments of 0.2(mm) leading to 1250 points. Robot undergoes a 3 (sec) pause on each point to allow the tracker to record the position of end-effector and also to avoid any Figure 3.37 states that the fluctuations of the mentioned straight path has been reduced from 0.15(mm) (Corresponding to a maximum deviation of 0.09(mm) from the path origin) to 0.09(mm) (Maximum deviation of 0.05(mm)) while the repeatability of the robot is claimed to be 0.03(mm) by the robot manufacturer. The difference between the targeted path accuracy and achieved one can be described by three main points:

• The effect of gear transmission error of wrist joints were not identified and compensated yet the wrist was involved in the line motion.

• The stated repeatability by manufacturer is conditioned to a constant temperature but since the straight line test takes more than an hour of constant movement of the robot, the repeatability is expected to increase even though the ambient temperature was constant [START_REF] Heisel | Thermal behaviour of industrial robots and possibilities for error compensation[END_REF] .

• One should note that manipulation with the joint variables is limited to the unmodifiable parameter of resolution value of each joints causing a considerable round-off error in the range of estimated gear error orders.

60 Chapter 3. Identification And Compensation Of Non Geometrical Robot Errors

Conclusion

In this chapter, efforts were dedicated to identifying and compensating the compliance behaviour of robot structure and cyclic joint transmission systems as main contributors to non geometrical robot pose errors. Due to the considerably higher influence of compliance behaviour on robot inaccuracy, this source was surveyed priorly in section 3.2.

To model the compliance behaviour, we made use of Virtual Joint Method (VJM) approach to generate an elastostatic model according to its convincing performance presented in the existing literature. The generated elastostatic model involves both self gravity and external force effect. Section 3.2.3 was attributed to self gravity effect in which the procedure of implementing the mass model of the robot (Links weights and coordinates of centers of gravity), gravitational torque imposed on robot joints for a given configuration and effect of balancing system attached to robot were thoroughly explained.

Due to the necessity of joint stiffness values in the elastostatic model, experiments were conducted to determine these values. Section 3.2.5 presented the procedures of force imposition and wielded measurement devices for each joint. The nonlinear effect of the self-balancing system on the stiffness of joint 2 was also modeled in this section using a varying stiffness torsional spring. However, according to investigations, its effect was minor compared to the stiffness of joint 2 and thus neglected for further considerations.

To verify the generated elastostatic model, experiments were conducted by using the load compensation option offered by the robot arm supplier Stäubli company. Since this compensation module is designed to counter care the self-gravity and additional device attached to robot flange (Spindle or ultrasonic wave generator e.g.), it offered an opportunity to validate our elastostatic model. The verification process revealed the validity of the elastostatic model. However, the stiffness values obtained in the experimental procedure were replaced with the pair provided by Stäubli due to the more accurate simulations using the latter.

In Section 3.3 cyclic transmission system errors were focused. A novel method was proposed to mathematically model this error based on an optimal path to reduce the time of the identification compared to the existing method in the literature. An offline compensation method was proposed base on estimating and eliminating the cyclic errors for a required path in the workspace of the robot. Experiments on a lab-scale robotic arm (Stäubli TX90) showed that the proposed identification and compensation process was capable of reducing the amplitude of these fluctuations by twice the repeatability of the robot. The method can be used to increase the path accuracy of the robot manipulators especially those wielded for machining applications where path deviations cause undesirable quality imperfections.

Chapter 4

Robot-Tool Compliance Treatment In Machining Process

Introduction

The intervention of Industrial Robot (IR) arms in machining applications is rapidly widening its scopes. IRs which were used only for a handful of applications mainly based on their repeatability (e.x. pick and place and painting) are now replacing dexterous operators and precise machining devices.

A machining robot cell contains at least three elements [START_REF] Olabi | Amélioration de la précision des robots industriels pour des applications d'usinage à grande vitesse[END_REF]: Robot arm, machining tool, and workpiece. The design of the robotic cell is highly dependent on the targeted machining process. Some processes such as milling metal blocks come with high machining forces and thus extra care should be exercised to improve the machining accuracy by means of compliance behaviour treatment of robot. However, in applications of such type, the compliance behaviour of machining tools, e.g. drill bits is negligible compared to the portion introduced by the robot arm.

Imperfections observed in the resulted workpieces in a robot-involved machining process come from different sources. On the robot side, these errors, as mentioned in previous chapters, arise from geometric and non geometric sources. But as reported earlier, the compliance behaviour of the robot is the dominant cause of tool misplacement for a loaded robot arm for a heavy IR. Moreover, in some machining applications, the effect of machining tool compliance can not be neglected in manufacturing errors.

In this chapter, we aim to address a procedure for compensating errors induced by the compliance behaviour of a robotic-based machining process in which the compliance effect of the machining tool can not be neglected compared to the one of robot structure. While literature is thorough with the machining processes where tool compliance is not an issue (such as drilling and milling), the goal of this chapter is to introduce an outline to counter the compliance effect of the machining tool as well as robot structure in machining processes combining compliant tool.

As a case-study machining process, the chamfering process on honeycomb material was adopted since this machining process comes with a compliant tool and also was in line with the objectives of the industrial partner of this work Le Creneau Industriel. The remainder of the chapter is organised as follows:

Firstly, a state of the art concerning robot machining topic and approaches proposed by recent research works for increasing the accuracy of the robotic involved Chapter 4. Robot-Tool Compliance Treatment In Machining Process machining process is presented in section 4.2. In section 4.3 we take a closer look at the case study machining process, namely the ultrasonic cutting process of honeycomb materials. In this section, the advantages of these materials along with the challenges introduced by their complex structure in manufacturing and numerical simulations of behaviour of these materials are noted.

In section 4.4, a detailed overview of this machining process in Le Creneau Industrial is presented with the description of different elements involved in this process. Section 4.6 presents a quantification of the machining errors involved in the ultrasonic cutting process of Aramid honeycomb materials. The measurement procedure of these so-called machining errors is explained in section 4.8. Data capturing strategy and measurement devices used to construct a machining force model are presented in section 4.9. To observe the effect of robot structure compliance effect in ultrasonic machining errors, section 4.10 is dedicated to the elastostatic simulation of robot structure in a certain machining condition with experimentally-driven machining conditions.

To introduce the compliance effect of the machining tool, section 4.11 details the procedure of the compliance behaviour modeling of the cutting knife using structural stiffness matrix approach. The compliance behaviour of the entire robot-tool structure is surveyed in section 4.12 to evaluate the potion of contribution of each part in the observed machining errors. To march through a more accurate simulation of knife behaviour, in section 4.13 we describe a procedure for simulating loading condition on the cutting knife surface. Section 4.14 presents a compensation method to decrease the machining errors in the ultrasonic honeycomb cutting process. Experimental results validating the approach are presented. The chapter ends in section 4.15 with noteworthy conclusions achieved throughout the chapter.

state of Art in Robotic Machining

Nowadays, robotic-involved machining operations cover processes such as trimming (Slamani and Chatelain, 2019), grinding [START_REF] Chaoui | Improving Surface Roughness in Robotic Grinding Process[END_REF], milling [START_REF] Zhang | Investigation on the milling performance of amputating clamping supports for machining with industrial robot[END_REF] etc. These processes usually involve high level interaction couples provoking robot compliance behaviour in both dynamic and static mode.

Garnier, Subrin, and Waiyagan, 2017 have generated a model for robotic-involved drilling application. They have included different parameters of the drilling process such as tool characteristics, feed rate and bit rotational speed to develop a model estimating machining couple (force and torque) imposed on tool and consequently robot arm. They have also proposed an algorithm to increase the milling accuracy by redundancies involved in the process based on elastostatic behaviour of robot structure.

In [START_REF] Wang | Base detection research of drilling robot system by using visual inspection[END_REF], researchers have armed the operating milling robot with a vision system and a laser displacement sensor. They have developed a method based on mentioned accessories to update the transformation matrix between milling workpiece and designed milling device mounted on a Kuka robotic arm to increase the milling accuracy by decreasing the errors involved in robot base-workpiece relative position definition.

Slamani and Chatelain, 2019 have attributed the final manufacturing imperfections of Carbon Fiber Reinforce polymers (CFRPs) to three different error natures. These sources which can be generalized to any machining process performed by IRs, are namely robot-related, process-related and material-related errors.

Robot induced errors which have been deeply studied in previous chapters, rise from geometric and non geometric sources. Process-related errors cover the sources which are process dependent such as machining forces, dominant force imposition frequencies and machining tool compliance effect. Material-related errors involve the misbehaviour of workpiece in the specific machining conditions. As an example, some materials such as CFRPs, possess low thermal conductivity and thus the increase of local temperature in tool-material contact zone provokes errors that are described in (Slamani and Chatelain, 2019).

The errors induced from compliance behaviour of robot and/or machining tool is trivially proportional to the tool-workpiece interaction couple. This is also evident that the interaction couples are highly dependent to the nature and characteristics of machining process in both magnitude and direction. In light of mentioned points, the graph presented in Figure, explains the parameters affecting the compliant-induced errors for a robotic involved machining process.

Machining process of Honeycomb Ultrasonic Cutting

Honeycomb composite structures offer critical characteristics for designers and developers of recent technologies. High strength along with light weight makes them a promising element to increase the fuel efficiency in transportation sectors and defining their wide application in aerospace sector. Moreover, energy absorption for both sound and impact, self-extinguishing of fire, thermal and electrical insulation and dielectric strength are other appealing properties of these materials made them to be extensively used in electric vehicles, surfboards and high-performance boats [START_REF] Gill | Studying the Mechanisms of High Rates of Tool Wear in the Machining of Aramid Honeycomb Composites[END_REF].

Honeycomb composite materials are generally composed of three different core materials: Ceramic ones for high-temperature applications such as catalyst carriers and heat exchangers, Polymer cores such as Aramid materials which are mainly used in components contacting fluid flow and Metal ones such as aluminum are used from cheap doors to advanced aerospace components such as energy absorbing parts [START_REF] Gibson | Cellular solids: structure and properties[END_REF].

In most uses of honeycomb materials, designers should limit themselves to uniform thickness of honeycomb parts. This is due to the fact that using sculpted surfaces will dramatically increase the machining costs. But on the other hand, since a wide zone of applications involving honeycombs especially Aramid ones deal with complex surfaces designed to represent flow stream (such as radomes, engine cowlings, fairings, helicopter blades [START_REF] Gill | Studying the Mechanisms of High Rates of Tool Wear in the Machining of Aramid Honeycomb Composites[END_REF]), furthering the use of these materials is dependent on decreasing their machining costs and defeats.

Difficulties in machining of Aramid honeycombs mainly originate from the their structural complexities. Aramid honeycombs are reinforced by high performance fibers such as Kevlar and Nomex [START_REF] Seemann | Numerical modelling of nomex honeycomb cores for detailed analyses of sandwich panel joints[END_REF]. [START_REF] Gill | Studying the Mechanisms of High Rates of Tool Wear in the Machining of Aramid Honeycomb Composites[END_REF] have pointed out several machining challenges of Aramid fiber reinforced composites because of delamination of fibers from the matrix, incomplete separation of machined material (flagging) leading to re-machined materials, deflection of the honeycomb wall structure and brittle mode machining of the matrix and the embedded fibers. They have focused on the tool wear effect in machining of Aramid honeycomb materials. They examined the hypothesis that the high rate of tool wear in Aramid honeycomb machining is due to the relatively low stiffness of the material causing excessive rubbing and heat generation on the tool. Using high speed video technology, they observed that a shredder device was completely unusable after 96 minutes of cutting due to the wear effect.

To achieve an efficient model able to describe the non-homogeneous behaviour of honeycomb materials, researchers have chosen different approaches such as analytical, experimental and numerical methods. Research and observations over the behavior of honeycomb materials started in 1947 by Norrise and Charles in [START_REF] Norris | An analysis of the compressive strength of honeycomb cores for sandwich construction[END_REF]. However, Gibson and Ashby have provided in [START_REF] Gibson | Cellular solids: structure and properties[END_REF] a thorough source of analytically and experimentally derived equations for a vast number of phenomenon in the behavior of honeycombs such as elastic deformations, buckling, creep and creep buckling, different failure causes and viscoelastic deformations. Abrate, Castanié, and Rajapakse, 2012 have provided a comprehensive overview of the topic categorized by their industrial applications and since it is a comparatively recent work, consists many numerical simulations. But they do not 4.4. Honeycomb Ultrasonic Cutting Process in Le Creneau 65 talk about machining forces and effects since these effects are more raising from local strain results of the machining forces but the theoretical equations are driven for the global loading of a material structure. An observation of local deviations was done by [START_REF] Mora | Measurement of the Cosserat constant of circular-cell polycarbonate honeycomb[END_REF] explaining the behavior of the honeycomb using Cosserat theory. But they used polycarbonate circular cells to construct the honeycomb material base and the loadings were also uniaxial and not local.

Main obstacle for numerical simulations of honeycomb materials, is the buckling phenomenon in the cell walls under local loading which initiates the failure of the structure. To overcome this problem, wall cells are often modeled with two or three dimension elements to represent the thin wall cells of the structure. Such modeling is also referred to as meso-scale modeling. A recent and comprehensive review and contribution over the different simulation approaches is published by Seemann and Krause in [START_REF] Seemann | Numerical modelling of Nomex honeycomb sandwich cores at meso-scale level[END_REF] along with experimental validation using standard tests. They also mention that apart from the complexity of the element constructions for numerical simulations, defining the mechanical properties of aramid papers constructing the wall cells is an issue. This issue mainly raises from the manufacturing processes notably high-temperature calendaring causing Young's modulus to vary in different in machine and cross machine directions. [START_REF] Jaafar | A 3D FE Modeling of Machining Process of Nomex R Honeycomb Core: Influence of the Cell Structure Behaviour and Specific Tool Geometry[END_REF] have simulated and experimented the milling process of Nomex honeycomb material to observe the machining forces and surface quality. Regarding the surface quality, they have mentioned two main defeats which are mainly uncut aramid fibers along the machined surface and tearing of the walls. Cutting process of honeycomb introduces even more shortcomings. Cutting Aramid honeycomb materials with sharp blades are time consuming and environmental polluting (Xiang et al., 2019a). Introducing ultrasonic cutting technology has overcome many of the encountered challenges. Researchers have reported better surface quality in terms of burr generation and fiber tearing defeats as well as lower machining forces while using replacing ordinary cutting with ultrasonic-assisted cutting process (Xiang et al., 2019a;[START_REF] Xiang | Ultrasonic longitudinal-torsional vibration-assisted cutting of Nomex R honeycomb-core composites[END_REF].

Honeycomb Ultrasonic Cutting Process in Le Creneau

Le Creneau Industriel is one of the few companies specialized in ultrasonic cutting technology. They offer solutions for disk and knife ultrasonic cutting tools used for wide range of materials from soft Aramid Honeycomb materials to metallic honeycombs.

Ultrasonic cutting solution offered by Le Creneau contains different essential elements and facilities. A complete solutions is a combination of different assembled elements. Figure 4.2 depicts a robot arm, armed with ultrasonic cutting assembly.

As illustrated in Fig. 4.2, the assembly is composed of a spindle device which takes control over the power supply and switching of the ultrasonic head. Ultrasonic head generates reciprocating movements along its axis with a frequency of about 20(kHz) to ensure a soft cutting edge, a cutting knife and a support which attaches all the assembly to the flange of the robot. Honeycomb row materials usually come in sheet-like bulks and the main applications of ultrasonic cutters are to extract work pieces with simple to complex geometries. Different geometries are achieved with aim of two main cutting tools: Ultrasonic Knives and Ultrasonic Disks.

Honeycomb Chamfering Process

Honeycomb sandwich panels are widely used in applications such as intake barrel panels where beveled edges are inevitable [START_REF] Cunningham | The effects of various design parameters on the free vibration of doubly curved composite sandwich panels[END_REF]. Beveled edges usually terminate in a solid panel edge strip around the circumference of the panel to ease of attachment to the neighbouring structures with fasteners [START_REF] Soovere | Random vibration analysis of stiffened honeycomb panels with bevelededges[END_REF]) (Figure 4.4).

To form a raw honeycomb piece into a desired geometry fitting into the ending form of the panel, Chamfering Operation is used. In a Chamfering operation, a portion of the honeycomb material is removed in such way that an initially right angled face is turned into an acute angle. The final angle between the lower edge [START_REF] Soovere | Random vibration analysis of stiffened honeycomb panels with bevelededges[END_REF] and the bevel edge is termed as Chamfer Angle and is trivially the most important feature of the chamfer. Figure 4.5 depicts a chamfer with a chamfering angle of 30 . In chamfering process using ultrasonic technology, both ultrasonic knives and disks can be used. However, one should note that the length of knife is longer than the radius of an equivalent disk and subsequently, performing a chamfering operation with ultrasonic knife is several times less time consuming thus more of interest in spite of notable advantages of disk tool over knife in resulted surface quality [START_REF] Xiang | Ultrasonic longitudinal-torsional vibration-assisted cutting of Nomex R honeycomb-core composites[END_REF].

To operate a chamfering process with ultrasonic knife, similar to any other machining operation, different machining parameters are defined. These parameters are vital for path and tool orientation data generation in any CNC programming software or interface.

All throughout this work, we refer to different parameters of cutting process of honeycomb material. In this part we define all the terms and definitions of these parameters. A general schema of cutting process of honeycomb for a chamfering operation is depicted in Figure 4.6. Figure 4.6 shows a cutting knife operating to form a chamfer with a desired angle of β (chamfering angle) from the horizontal plane. To ensure the best surface quality, one of the important constrains of programming the operation, is to coincide the lateral cutting surface on the knife with machining plane (Shown in Figure 4.6). Since the geometry of the cutting knife consists of a sharpened "V" edge with a certain angle, the deviation of the tool axis, for having the two mentioned planes coincide, differs from the complementary vale of the chamfering angle (α = π 2β). The angle between vertical axis of the work-piece and the tool axis (α) is termed the Lean angle. Due to reasons such as directions of forces and cutting depth, knife edge is tilted toward the cutting direction of the material forming an angle of θ between the tool axis and the vertical work-piece plane which is referred to as Lead angle.

Machining Errors in Honeycomb Chamfering Process

All the existing research works in the literature in the field of ultrasonic machining imperfection of Aramid honeycomb materials, have focused on microscopic errors such as burr generation, fiber tearing and surface quality (Xiang et al., 2019a;[START_REF] Xiang | Ultrasonic longitudinal-torsional vibration-assisted cutting of Nomex R honeycomb-core composites[END_REF][START_REF] Hu | Research on curved surface forming of Nomex honeycomb material based on ultrasonic NC cutting[END_REF][START_REF] Zhang | Characteristics of surface microscopic of Nomex honeycomb after ultrasonic assisted cutting[END_REF]. In this work, we dedicate our effort to model, quantify and compensate the manufacturing defeats in terms of imperfections observed in resulted work-pieces in comparison with desired piece geometry. In light of the points noted regarding a resulted workpiece of honeycomb ultrasonic cutting process depicted in Figure 4.8, we divide the geometry of a resulted pieces into three different zones. These zones are namely Entrance Zone, Steady State Zone and Exit Zone. These zones are illustrated in Figure 4.9 where a schematic model of a resulted piece is presented. As shown in Figure 4.9, the machining errors of a honeycomb workpiece material vary from one zone to another. We attribute the difference of the resulted error Chapter 4. Robot-Tool Compliance Treatment In Machining Process on each zone to the counteracting tool-workpiece machining load. In the entrance and exit zone, the contact area between the knife and material varies as function of the TCP placement and so does the machining force. On the other hand, the steady state zone signifies the portion of the workpiece in which the contact configuration between the workpiece and cutting knife stays unchanged.

Knife entrance zone

S t e a d y S t a t e Z o n e E n t r a n c e Z o n e E x i t Z o n e

The dimensions of the entrance and exit zones are only a function of the machining features such as lead and lean angle, however, the length of the steady state zone is the trivially the subtraction of the entrance and exit zones of the total workpiece length. Thus depending on the required dimensions of the workpiece, the portion of each zone in the final resulted workpiece varies. But having in mind that the length of the entrance and exit zones will vary in the range of several centimeters (depending on the machining features), the main concern raises for the steady state zone errors. We should also mention that by having an efficient procedure to compensate the errors of the steady state zone errors, the designers of the machining process can eliminate the error of the entrance and exit zone errors by simply adding an offset length to the required machining path to ensure that the required length falls within the steady state zone of knife-workpiece contact. To proceed to introduce different geometrical imperfections of a resulted workpiece, a close-up view of the cross section of the resulted workpiece in steady state zone is brought in Figure 4.10. • Angle β : The actual chamfering angle in the steady state zone.

• Displacement d 1 : The deviation of the cut line on the bottom surface of the work-piece in the straight line zone from the programmed (desired) cutting line.

• Displacement d 2 : The deviation of the cut line on the upper surface of the work-piece in the straight line zone from the programmed (desired) cutting line.

The other phenomenon faced in the experiments is that, although all the tests were programmed in a way that a gap of 1.5(mm) remains between the tip of the knife and the bottom surface of material (Distance l shown in Figure 4.10), while machining the honeycomb material after a certain value of chamfering angle, the knife had penetrated not only inside the adhesive layer, but also the protective moquette. Figure 4.11 presents the resulted defeat on the protecting bed after a honeycomb chamfering process. The first hypothesis seemingly capable of explaining the observed errors in side honey comb pieces after chamfering operations is the robot compliance behaviour. In other words, the main contributor to the errors is robot deviations due to the machining forces. However, to completely validate this hypothesis, one should provide answer to the following questions:

• What is the exact error values for a given operation parameters?

• What are the corresponding machining forces applied on the robot?

• Are the machining forces capable to deviate the tool path and orientation to explain the corresponding manufacturing errors?

To answer the above mentioned questions, an experimental test session was executed by variable control method. In this set of tests, some features will remain constant due to the machining policies of Le Creneau, and the important design parameter (Chamfering angle) is varied. These tests must provide the answer to the first two questions. To have the exact values of the machining forces, a force sensor is used to capture the machining forces. The resulted workpieces are then used for error measurements which are explained in further steps. Below comes a detailed description on the test parameters, conditions and involved measurement and performing devices.

Experimental setup for Honeycomb Chamfering process

To achieve a relationship between chamfering design parameters and different elements on the final manufacturing errors observed on the honeycomb material chamfering process, a set of tests are designed to be performed using Stäubli TX200 Industrial robot arm. The captured data of these experiments are interacting machining forces and geometrical errors observed in the resulted workpieces.

The protocols of the test session were as following: Different honeycomb materials were cut using different geometry properties and cutting conditions. Each time, the machining forces were measured using a force sensor and the corresponding machining errors were measured using Faro arm. The force sensor we used in this test sessions was a Kistler Type 9257B Multi Component Dynamometer. The position of the force sensor and working table with respect to the robot base is depicted in Figure 4.12.

Kistler Table FIGURE 4.12: Placement of the force sensor regarding robot base Workpieces were placed on the Kistler table using an adhesive layer to prevent displacement. Sensor surface was also covered by a layer of moquette to increase the safety factor to avoid the collision of the cutting knife tip and the sensor surface. By executing each test, the sensor starts to record the interaction forces between the knife and honeycomb work-piece in a defined coordinate frame. Sensor coordinate system and the mentioned additional elements are shown in Figure 4.13.

As shown in Figures 4.12 and 4.13, the size of the Kistler sensor is limited. Therefore the size of the honeycomb workpieces was about 15(cm) × 12(cm). The surface of the sensor is covered by a layer of moquette protecting the surface. An adhesive layer (shown in red) is attached to the moquette to firmly hold the material. However, in certain tests, the adhesive layer could not resist the high values of the chamfering process.

Examined Honeycomb Materials

During the test session, we used two different honeycomb materials. One was a Nomex honeycomb with density of 38( These two honeycomb materials are shown in Figure 4.14. The workpieces shown in figure 4.14, are actually resulted pieces from two cutting conditions. Each workpiece was subjected to two chamfering process with same machining characters. The goal of executing two processes on each workpiece was to observe the compliance of captured force and error data from the two identical processes.

Cutting Tool

The cutting knife was another variable parameter during the tests. Two different knives were used through the tests with same geometries. The height of the knives were 85(mm) from the root to the tip. The difference between the two sets of knives were their manufacturing company. Knives were either provided by Dukane company termed hereafter as N1 or PRACARTIS company referred to as N2 in this work.

We should note that, the cutting knife in ultrasonic application has a limited lifetime and it is prone to fracture mainly caused by fatigue due to the varying tension followed by ultrasonic waves. For the above described tests, machining forces were recorded by the Kistler dynamometer while executing the tests. Geometrical manufacturing errors were also measured after removing the resulted pieces from the sensor. The procedure of measuring the errors observed in the workpieces and treatment of the force data captured during the test execution to achieve a machining force model are reported thoroughly in the coming two sections (sections 4.8 and 4.9).

Cutting Strategy and Parameters

Error Measurement Procedure

Measurements on honeycomb materials are always challenging due to their thinwall and low-rigidity structures. To avoid undesired contact effects on temporary deformation of the honeycomb body during measurements, researchers have proposed either low-contact force methods or contactless measurements devices [START_REF] Qin | Burr removal from measurement data of honeycomb core surface based on dimensionality reduction and regression analysis[END_REF] to accurately report the geometrical dimensions of such materials. These methods involve expensive facilities [START_REF] Alblalaihid | Fabrication and characterisation of a novel smart suspension for micro-CMM probes[END_REF] and software [START_REF] Neamţu | Training in coordinate measurement using 3D virtual instruments[END_REF] and/or wield sophisticated numerical post processing [START_REF] Qin | Burr removal from measurement data of honeycomb core surface based on dimensionality reduction and regression analysis[END_REF] which are out of the scope of our equipment and core research subject. Thus in this work, we have developed an easy-to-apply on site measurement protocol for measuring different error types. The protocol is described thoroughly in this section.

In this part we explain how we proceed to measure the errors on the resulted work-pieces. According to the previous section, we are dealing with different geometrical errors which do not have the same approach for measurement. The measurement procedure for the three different error types illustrated in Figure 4.10 plus the chamfering angle at the entrance point of the workpiece are reported below.

Entrance Angle

Entrance angle is the angle with which the knife initially penetrates the honeycomb workpiece. This angle is only found on the entrance and therefor does not form a plane. To measure this angle we used a manual high accuracy goniometer as shown in Figure 4.16. The entrance angle was accurately equal to the desired chamfering angle in both materials and knives.

Steady State Chamfering Angle

The steady state chamfering angle is termed to the angle which the ramp generates with bottom surface in straight line cut zone (β in Figure 4.10). Since this parameter is our main concern due to application of Le Creneau, this angle has been measured using a more accurate device and process. We made use of a Faro arm device to measure this angle. Each time, the work-piece was fixed on a marble table whose flatness is within microns. Pieces were mounted on the table using two vises as holders. By touching the marble table with the end-effector of the Faro we constructed the table plane. Due to the geometry of the honeycomb, we used a thin metal sheet to cover the surface of the chamfer to avoid the end-effector sphere fall into local holes created by the cells. By touching several points on the metal sheet, we construct another plane (Chamfer Plane shown in Figure 4.9 in red) which forms the angle β with the table plane. This process was automatized by AutoDesk Software. By measuring one piece of the hard material and one of soft material several times with the mentioned protocol, the repeatability of this measurement protocol was ±0.5(deg) for the hard material and ±1.5(deg) for the soft material. Two main reasons can be mentioned: By fixing the soft piece each time, since the interaction force between the clamp and the piece was different, the visible deformation leading to change all the geometrical aspects was also different. Second reason was the local deformation caused by the applied force from the Faro end-effector to the chamfer surface. The so called local compliance behaviour of soft honeycomb material was witnessed to be responsible for lower accuracy of the angle measurement process for this material.

Displacement Errors

Measuring displacement errors (which are namely d1 and d 2 in Figure 4.10) introduced some other obstacles. For d 1 , the main problem arose from the fact that the very weak tiny cell walls in the entrance zone were mostly damaged, deformed or torn away from the piece while removing the piece from the adhesive layer. This made it breathtaking to decide an exact point which we can refer to as the entrance of the point. The rest of the procedure was only to extend the straight line and measure its normal distance to the mentioned point by a ruler.

Measuring d 2 was even harder: The upper surface does not provide an accessible room for the point of a pen to create a line or leaving a sign for the entrance while placed upside down lying on the table. The procedure was to visually follow a certain row of cells. Figure 4.19 presents these error values for the examined chamfering angles. We should note since the measurement process of displacement errors does not possess a high level of precision, the presented data in Figure 4.19 can be used preferably only for comparison uses.

In this section we presented the errors measured in the machining process of honeycomb ultrasonic cutting. One can easily conclude that the manufacturing errors reported in this section, regardless of the type, increases when the chamfering angle (β) decreases. This clearly shows the augmenting trend of machining forces when the critical design parameter (β) decreases. This indicates that the machining accuracy of this process does not satisfy the tightly-toleranced pieces required by aeronautic field (the required tolerance is usually within ± 1. 5 [deg]). To investigate our effort to find the main contributor of the observed errors, in next section a simulation of the robot compliance behaviour in corresponding loading condition is presented.

Machining Forces in Honeycomb Ultrasonic Cutting Process

Due to the increase of the use of honeycomb core materials especially in sandwich structures, many researchers have dedicated their effort to generate explicit expressions and models able to predict the machining forces for both ultrasonic and ordinary cutting scenarios ( [START_REF] Hu | Research on Cutting Force Model of Triangular Blade for Ultrasonic Assisted Cutting Honeycomb Composites[END_REF]Xiang et al., 2019a;[START_REF] Xiang | Ultrasonic longitudinal-torsional vibration-assisted cutting of Nomex R honeycomb-core composites[END_REF]). However, these models require different property constants of the honeycomb materials which are not conveniently available and moreover, only capable of predicting the evolution of the machining forces by modifying certain machining parameter. In other words, an accurate model, capable of expressing machining forces, within an acceptable error range, in Aramid honeycomb ultrasonic cutting process is lacking in the literature. Thus we chose an empirical approach to build such model based on specific machining conditions of Le Creneau company.

The raw data captured by the force sensor (Kistler dynamometer) in previous section is defined in the sensor's coordinate system and possess some notable characteristics. Figure 4.20 shows the evolution of machining forces in ultrasonic cutting operation for hard honeycomb material in chamfering angle of 50 [deg] defined in the Kistler sensor coordinate system. .20 reveals several notable characteristics of the honeycomb ultrasonic operation forces using knife tool. There are two zones in which all the components of machining forces are null. These zones correspond to the time before the first knife-material contact and after complete exit of the knife. Knife entrance zone depicted in the figure, corresponds to the time window between the first contact of the knife with honeycomb work-piece and the instance when knife has completely penetrated the material. In this zone, the active knife length increases by time and so do the machining forces.

F z max F z min
All three force components shown in Figure 4.20 fluctuate in time. These fluctuations reveal the non-homogeneous nature of the honeycomb material resulted from variation of the number of Aramid polymer sheets by marching the knife along cells. The amplitude of these variations can go up to about the average value of a certain force component, thus for each component two parameters are defined which are namely minimum and maximum of absolute force value (for example Fz min and Fz max for Z component). To prepare for the extreme conditions, we will choose the maximum value of each component for further steps.

One should note that the data capturing process is triggered manually by operator and it's a selective parameter which will not impact the results as far as the process is triggered before knife-material meshing instance. Nevertheless care should be exercised to trigger and stop the data acquisition in proper timing to avoid bigsized data cells.

For applications such as simulating the compliance behaviour of the robot structure under loading condition of honeycomb ultrasonic machining process, one can use either the force vector defined initially in Kistler coordinate system (By assuming to have the transformation between the Kistler to the robot base frame) or more conveniently, machining forces defined in the tool frame which is equivalent to the TCP frame of the robot.

The presented data in this part, were data corresponding to a set of experimental tests carried out on a hand full chamfering angle on both soft and hard honeycomb materials. The required chamfering angle for final shape of the honeycomb material, however might be a different angle other than the ones experimented empirically. Thus we decided to build an estimating model, based on the gathered experimental values to enable us to estimate the machining forces for any arbitrary chamfering angle within the practical range. The proposed generic force expression for any direction of the both workpiece frame and tool frame is assumed to have the following format:

F c cos(β) + F s sin(β) + F 0 (4.1)
To find the appropriate values of F c ,F s and F 0 Least Square method was applied on the achieved experimental data in the described force measurement procedure. Following expressions estimate the machining force for a given chamfering angle of β expressed in the Kistler table coordinate system presented in Figure 4.22 for Hard honeycomb material : Using same approach for the Soft honeycomb material results in the following expression for forces induced in chamfering operation:

   f x (β) =
   f x (β) =
31.9 sin(β) + 2.6 cos(β) -38.9(N) f y (β) = -10.5 sin(β) + 14.4 cos(β) + 11.8(N) f z (β) = As one can easily conclude from the forces evolution shown in Figure 4.21, the machining forces in chamfering operation on honeycomb material using ultrasonic knife, increase in magnitude in each direction regardless of the workpiece material. Figure also reveals that the machining forces in soft material is about 5 times less compered to the ones of hard material. The conclusions comply with the geometrical errors evolution resulted from the error measurement procedures in section 4.8.

Defining the machining forces in the Kistler coordinate system comes with advantages such as avoiding additional transformation matrices when using general jacobian matrix for robot elastostatic behaviour. But however, for simulating the compliance behaviour of knife structure, which is one of our objectives through out this chapter, machining forces defined in knife coordinate system is inevitably necessary.

In applications where we aim to observe the behavior of knife in real cutting conditions we make use of the force data recorded in the Kistler table coordinate system. To proceed, we need to translate the force vector into the knife frame which is composed of a Lateral knife direction, Longitudinal direction which is equivalent to the tool axis, and third will be the common normal of the aforementioned vectors, also termed as Vertical direction of the knife. Knife frame is shown in red in Figure 4.22, where axis X, Y and Z represent the longitudinal, lateral and vertical directions of knife respectively.

To convert the recorded data into knife frame, a new rotation matrix R Kn f Kist is defined to transform the force vectors captured in the coordinate system of Kistler table to the one of cutting knife. Expression below was used to compute the latter:

R Kn f Kist = R Stb Kn f -1 * R Stb Kist (4.4)
Where R Stb Kist is a constant rotation matrix from robot base frame to the Kistler coordinate system and R Stb Kn f is the rotational part of the transformation matrix between To build a model of machining forces expressed in knife coordinate system same approach is used with the one explained for Kistler table. For Hard material, the machining forces defined in Knife coordinate system are as following:

   f Longitudinal(x) (β) =
106.3 sin(β) + 37.8 cos(β) -118.7(N) f Lateral(y) (β) = -301.8 sin(β) -84.2 cos(β) + 331.8(N) f Vertical(z) (β) = -391.9 sin(β) -98.9 cos(β) + 459.1(N) (4.5) Adopting same approach for the machining force of chamfering process in Sotf honeycomb material, gives the following model:

   f Longitudinal(x) (β) = -14.0 sin(β) -4.2 cos(β) + 13.9(N)
f Lateral(y) (β) = -0.1 sin(β) + 17.1 cos(β) + 0.3(N) f Vertical(z) (β) = -75.9 sin(β) -13.7 cos(β) + 84.9(N) (4.6) Figure 4.23 depicts the experimental data and generated model estimating the machining forces in chamfering process of both soft and hard honeycomb material defined in Knife coordinate system presented in Figure 4.22.

According to the graphs exhibited in Figure 4.23, the vertical part of the machining force is the dominant component applied on the cutting knife. Having in mind that this component is applied on the direction perpendicular to the knife body with the lowest thickness, further explains the importance of this component on the compliance behaviour of tool.

Another noteworthy point drawn from Figure 4.23 regards the force component applied on the longitudinal direction of knife structure. Due to the movement of the knife with respect to the workpiece, a negative value is expected for this force component. This expectation is satisfied from the machining forces in hard honeycomb material. However, in some points, a positive value is shown for this force component in soft honeycomb material. We must add that the force values on this points are between 0 to 0.5(N) which is negligible compered to the flactuational behaviour of machining forces observed in honeycomb chamfering process described earlier.

By having the exact values of machining forces and resulted errors in the honeycomb workpieces, necessary elements are in hand to evaluate the elastostatic behaviour of robot structure in the machining condition. Next section describes the simulation analysis of robot compliance behaviour in a chamfering process.

Simulation of Robot Behaviour in honeycomb Machining Process

In this section, we dedicate our effort on simulating the robot impact on errors caused by machining forces applied on robot TCP in the ultrasonic cutting operation on end-effector's orientation and positioning. Simulated mission is the machining operations of final piece shown in Figure 4.7 .

To program a machining mission, a certain procedure must be followed. Firstly, the machining features such as paths and curves are translated into the path and orientation of machining tool using an available Computer aided tools. Required movement of machining tool in then converted to the robot joint movement planning using the definition of robot end-effector, work-piece placement and orientation with respect to the robot base frame and robot inverse kinematics.

The joint values, TCP definition and machining forces applied on TCP, will enable us to simulate the compliance behaviour of the robot and consequently drawing conclusion of the error portion rising from this error source. Each step is brought in detail in the remaining part of this chapter.

Overall Configuration of the test

In this experiment, the row honeycomb block is placed on a vertical table stand supported by four columns screwed on the workshop ground. The table is located in front of the robot having the Robot's base frame x axis is normal to the surface of the table crossing it at x = 1678(mm). Placement of work-piece and the table stand with respect to the robot is shown in figure 4.24. Questions may rise on the accuracy of the actual placement of the supporting table defined in the software and the actual positioning of the table in workshop as a potential source of error. In real machining operation, the transformation matrix between robot base frame and the table is updated each time using high-accuracy LVDT1 senor before executing the machining operation. Another strong evidence for breaking this hypothesis is the accurate value of entrance angle in the previously described error measurements.

Top Right B o t t o m Left

TCP Definition

The machining head used for ultrasonic cutting consists of four assembles parts. These parts are spindle, ultrasonic head and cutting knife attached to the robot flange using a designed support. Head assembly is shown in Figure 4.25.

Definition of the TCP consists of determining 3 components of position of TCP defined in the 6th frame of robot ( the one attached to the surface of the robot flange) and the orientation of TCP frame (also termed as end-effector frame) leading to a TCP frame is constructed with the following procedure: Origin of this frame is the TCP point, z axis is along the spindle axis pointing to the interior of side of tool (This definition is necessary for the Siemens controller used for this operation) and y axis is parallel to the one of 6th frame of the robot. In this case, the final orientation of the TCP frame can be resulted of a rotation of -135(deg) of the 6th frame about its y axis. Figure 4.26 illustrates the details of the TCP frame. It is worth mentioning that this frame does not move with respect to the 6th frame i.e. by having the spindle rotating along its axis, the orientation of the frame stays steady.

T tcp 6 =     cos( -3π 4 ) 0 sin( -3π 4 ) 229.57 0 1 0 0 -sin( -3π 4 ) 0 cos( -3π 4 ) 457.84 0 0 0 1     (4.7)

Path Generation

Path generation is the process of determining Cartesian position and orientation of TCP frame with respect to the workpiece which will consequently lead to calculation of the joint values using the inverse geometrical model of the robot arm integrated with the achieved flange-TCP transformation matrix. Defining Cartesian position and orientation is carried out using a computer aided manufacturing (CAM) software. This software is also able to compute for the joint values but these results are only used to simulate the robot movement during the test to ensure no collision is likely to occur. Final output of this software is the Cartesian position and orientation which will serve as an input to the Siemens controller manipulating the robot to execute the final action.

Joint values resulted from the CAD software and the one resulted from our inverse kinematic model using the aforementioned TCP definition are equivalent. The starting and the ending point of the complete trajectory are identical since the tool completes the circumference of the square. This point is top left vertex and the tool follows the following sequence to execute the cutting process:

Top → Right → Bottom → Left
To calculate the joint variables of the robot arm for the machining process, both position and orientation of the TCP frame must be in hand. As explained previously, the generated data from the CAM software contains the Cartesian position and orientation of the TCP frame. The orientation is defined using Euler Angles with α, β and γ notation standing for rotation about x, y and z respectively. Where due to the Euler rotation presentation, the y axis is the resulted y axis of TCP from the first rotation of frame with α degrees about the x axis and z axis is the resulted z from the consecutive two rotations of α along x and β about y . Therefor the rotation section of the transformation matrix between the TCP and base frame can be expressed as:

R tcp 0 = R(x, α) * R(y , β) * R(z , γ) (4.8)
Expanding the rotation transformation matrices will lead to:

R tcp 0 =     cos (β) cos (γ) -cos (β) sin (γ) sin (β) sin (α) sin (β) cos (γ) + cos (α) sin (γ) -sin (α) sin (β) sin (γ) + cos (α) cos (γ) -sin (α) cos (β) -cos (α) sin (β) cos (γ) + sin (α) sin (γ) cos (α) sin (β) sin (γ) + sin (α) cos (γ) cos (α) cos (β)     (4.9)
The three mentioned angles generated by the CAM software for this experiment are shown in figure 4.27.

The fact which should be noticed here is that since the geometrical properties of the workpiece stays constant along a certain edge, the relative configuration of the cutting knife with respect to the edge should also stay steady. Thus we expect a constant rotation of the frame as the TCP travels along a certain edge of the square. While α and β components of the orientation are constant for each edge (As shown in figure 4.27), having the γ component varying dramatically especially for the edge "Top". This trend of the third component trivially changes the orientation of the TCP frame. But we recall the fact that it is not the knife which is attached to the TCP frame. The ending point of the knife is the TCP frame origin but knife is under the command of the spindle angle. In other words, it is the orientation of knife which matters and not the orientation of the TCP frame.

Knife frame can be assumed as a frame attached to the knife resulting from a single rotation along the z axis of TCP frame. The spindle rotation angle is also one of the outputs of the CAM software which hereafter will be termed as θ s . Since the spindle rotation direction is the reverse of the one of z axis of TCP frame, the final orientation of the Knife Frame can be computed as follows: Computing the joint variables for a given pose (position and orientation) of the TCP frame is achievable using the inverse geometrical model of the robot-end effector structure. Although all the non redundant robots are solvable for their joint variables for a desired pose [START_REF] Lee | Displacement analysis of the general spatial 7-link 7R mechanism[END_REF], but only those holding simple structures possess an analytical expression for their inverse geometrical models [START_REF] Paul | Robot manipulators: mathematics, programming, and control: the computer control of robot manipulators[END_REF]. Among robots with simple geometry, the wrist-partitioned robots stand as a special case. Since Stäubli TX200 robot is a wrist-partitioned robot arm, we made use of the inverse geometrical model of such robot types explained in Khalil and

R Fr 0 = R tcp 0 * R(z , -θ s ) = R(x, α) * R(y , β) * R(z , γ -θ s ) (4.10)

Loading condition

To simulate the robot behaviour in the test condition, using a comparatively accurate loading applied on robot TCP is trivial. In this operation, we face a four edge block, each undergoing chamfering operation for chamfering angle of β = 30 • (deg).

Efforts were dedicated to measure the machining forces corresponding to this angle with the previously mentioned procedure, but due to high level of forces, the adhesive substance, was not capable of handling the process leading to the piece removal from Kistler table.

Chapter 4. Robot-Tool Compliance Treatment In Machining Process

In this step, we make use of the described force predicting model resulted from the curve fitting process in the previous steps. The different force components described in Kistler coordinate system are calculated below:

f x = 132(N), f y = -131(N), f z = 89(N)
One should note that the generated force model used to estimate the machining forces of this chamfering process, was developed using data captured based on the configuration presented in section 4.7. Therefor, the direction of the computed forces must be modified based on the feed rate direction configuration of each edge.

Elastostatic Simulation

After determination of TCP, values of joint variables and machining forces, we are ready to simulate the action by the previously explained elastostatic model to simulate the error order we should expect from the compliance effect of robot structure.

The compliant induced errors can generally rise from three sources: The effect of machining forces, self-gravity effect of robot structure and the additional torque introduced by the weight of the assembly head on the joint torques. It is worth mentioning that the final error resulting from robot compliance is the superimpose of each mentioned source. These sources may counter-act and cancel out the effect of one another, thus the contribution of each is only judicious when surveyed separately.

As explained in chapter 3, robot controller in armed with a gravity compensation module. This module is capable of counter-caring the gravitational effect of robot links as well as tool assembly head. The only potential source of error is an inaccurate definition of the head assembly mass and CoG 2 to the robot controller by the operator.

As we are using the elastostatic model of the robot arm developed in chapter 3, simulations are able to take into account the three following effects: 1. Gravitational effect of robot elements (Self gravity effect), 2. Mass of the additional machining assembly and 3. The machining forces.

Recalling equation 3.15 from chapter 3, the angular deviation in robot joints (∆Θ) imposed by a set of torque values applied on robot joints (τ t ), can be computed using the following expression:

∆Θ = [K] -1 × τ t (4.11) Where [K] = diag(k 1 , k 2 , k 3 , k 4 , k 5 , k 6
) is a diagonal matrix contains the stiffness value of the robot joints. Set of stiffness values used in this step is listed in table 3.7. The calculated joint deviations were superimposed to the nominal joint value for the corresponding robot pose to achieve the actual joint value:

θ actual = θ desired + ∆Θ (4.12)
Eventually, the TCP displacement can be achieved by comparing the resulted TCP placement of the actual joint variables with the one resulted from desired joint variables using the forward geometrical model of the robot arm.

New set of joint values were fed to forward kinematics of the robot to compute the new orientation and position of the TCP by the assumption of the rigidity of the machining assembly. The deviation of orientation between the initial and final of the TCP frame is defined as the deviation of z axis of TCP frame. Consider z 1 and z 2 are the initial (desired) and final (deviated from the compliance effects) z axis of TCP frame. Thus the deviation in orientation can be expressed by the following expression:

D(rad) = cos -1 (|z 1 .z 2 |) (4.13)
Depending of the compliant sources taken into account to compute the imposed torque array (τ t ), the simulations were executed in three different conditions:

(A) First simulation is based on the assumption that the weight compensation option of Stäubli controller of the robot is able to fully compensate the selfgravity effect of robot elements and also the effect of gravitational effect of the mass of the machining head assembly. Thus this simulation practically shows the effect of machining forces solely.

(B) Second simulation takes into account of up to 20% of error in declaration of the mass value of the machining head. Mass of the assembly and its COG point is defined by the CAD model, potential differences between the actual mass properties and the one of CAD model can come from the differences in material properties of different elements in the real product and the modeled one, countless details of actual spindle which can not be implemented in the CAD model such as tiny electronic pieces, effect of wiring and its relative movement with respect to the assembly causing continuous displacement of COG etc. to observe any possible error rising from the mass declaration, second simulation takes into account up to 10(kg) of inaccuracy of defined mass of the assembly among the machining forces effect.

(C) The third simulations assumes no compensation strategy applied on the process. In this case, the self-gravity effect, machining force impact and total effect of the assembly mass on the robot compliance are included. The result of this simulation compared to the previous simulations will reveal the importance of the compensation strategies on the TCP displacement and deviation. A glimpse on the order of the results achieved in the simulation of TCP displacement and orientation error due to machining forces clearly reveals a vivid incompatibility with the experimental error results. As reported in section 4.8, the chamfering error is expected to rise progressively by decreasing of the chamfering angle. Thus, since an error of almost 8 [deg] is observed in the chamfering angle of 40 [deg], the expected chamfering angle is expected to be more than 8 [deg].

One should note that the angular deviation values resulted from the simulations reported in Figure 4.31 stand for the absolute deviation values. The deviations measured on the experimental part is only the projected potential angular deviation of the tool over the path which is only one component of the absolute deviation and trivially possesses a smaller value. But since the absolute simulated values are considerably negligible compered to the experimental ones, the projection part is skipped in this study.

By assuming the rigidity of the machining head, according to the CAM software, for an arbitrary point over the trajectory, joint 3 must undergo a deviation of 2. The non correlated results from the elastostatic simulations with experimental data well defines that by only considering the elastostatic behaviour of the robot arm we can not explain the observed geometrical chamfering errors in honeycomb workpieces. But however, before implementing the tool compliance model in the robot-tool behaviour, two potential error sources were studied:

• The inaccurate work-piece placement was initially invalidated due to two strong evidences: Firstly, the error measured in the entrance angle in all the workpieces was negligible. Another proof to support its invalidity was the difference in the measured error in different honeycomb materials for identical work-piece placement and machining features.

• The deformation of the work table supporting the kistler table and work-piece was invalidated during a measurement session by measuring its deflection induced by a normal force of 200(N) in work-piece zone using a manual gauge (Figure 4.33). The resulted deflections were negligible compared to the error order observed in honeycomb work-pieces.

In this part, we simulated the compliance behaviour of the Stäubli TX200 robot arm in machining condition of a chamfering process executed on an Aramid honeycomb workpiece. The simulation results were inaccurate of describing the expected machining errors. In next section, we aim to generate a model describing the compliance behaviour of the machining tool (Which is the ultrasonic cutting knife) to be integrated in a further step in the total compliance behaviour of robot-tool system.

Knife Compliance Behaviour

To achieve a model containing the compliance behaviour of robot and machining tool, a close-precise estimation of the behaviour of both parts is required. The compliance behaviour of robot structure has been thoroughly investigated in chapter 3.

In this part, we aim to focus our efforts on generating a model capable of estimating the deformation of the tool structure for a given loading condition. The approach proposed in this section is applied on our case study where the machining tool is a triangular cutting knife for honeycomb chamfering. However, this method can be adopted for any robot-involved machining application using compliant tool.

Machining tools are trivially designed for their attributed machining tasks having different characteristics. They widely differ in size, geometry and more importantly the material they are composed of. As described earlier, cutting process of honeycomb material in our application is performed using an ultrasonic mechanism. Due to the limited power of the wave generator implemented in the ultrasonic device, the mass of the knife is of great importance to ensure the required amplitude of the reciprocating movement in the knife structure. On the other hand, this constraint over the knife mass, prevents the designers to come up with larger dimensions and consequently stiffer designs. Thus ultrasonic knives suffer from considerably higher compliance compared to other machining tools such as drilling bits.

The importance of including the compliance model of the cutting knife is easily concluded in sight the above mentioned points as well as negligible errors from the robot structure as simulated in section 4.10 compared to the actual geometrical errors witnessed in the honeycomb workpieces. To generate a compliance model of the knife, we make use of Structural Stiffness Matrix (SSM) approach. This matrix is capable of estimating the deformation of the knife for a given loading condition.

The ultrasonic cutting knife is composed of a screw in one end to attach to the wave generator and a main body resulted from different trimming on a metal block. The final shape is a triangular plate with a thickness of 2(mm) and sharpened edges on the sides and tip. Figure 4.34 depicts the dimensions of the knife.

The structural Stiffness Matrix (SSM) is defined as the stiffness matrix relating the displacement of the knife tip to the applied load on the knife. The load, generally, can contain both force and torque components, but however, due to the nature of the ultrasonic cutting process, the torque components are neglected. Thus the As shown in Figure 4.35, a force with three components of f x , f y and f z is applied on the knife tip. Let us consider that the mentioned force causes a displacement of d = [δx T , δy T , δz T ] T to the knife tip. Stiffness matrix [K k ] relates the displacement d to the applied force:

  f x f y f z   =   k xx 0 0 0 k yy 0 0 0 k zz   ×   δx T δy T δz T   (4.14)
Due to the geometry of the knife shape, the stiffness values of knife structure is expected to vary from one direction to another. Thus for determining the elements of stiffness matrix, each direction is treated separately. One should note that the diagonal form of the stiffness matrix [K k ] is due to the fact that the coordinate system attributed to knife consists of its symmetrical lines.

For a generic notation, the x, y and z axes of the knife are termed respectively as Longitudinal, Lateral and Vertical directions of the knife structure. Further comes the compliance behaviour of the knife in each mentioned directions.

Knife Compliance Behaviour in Vertical Direction

Due to the loading conditions and attachment type of the knife on the ultrasonic head, a Clamped-Free triangular beam with a rectangular cross section is our proposed geometry for modeling the compliance effect in vertical direction. Thickness of the beam is set to be constant allover the knife surface and equal to t = 2(mm). The height of the triangle is chosen to be h = 70(mm) (midpoint between the knife tip and the ending point of flat zone of actual knife). But deciding about the initial base length of the triangle (b 0 ) needs more considerations. Figure 4.36 depicts an exaggerated deformation of the knife structure caused by a concentrated load f x at the right end. The illustrated loading causes a deflection function of z(x) representing the vertical displacement of the beam surface as function of x component defined in the frame shown in Figure 4.36. We recall governing differential equation in beam theory for the elastic curve of a beam subjected to any loading condition:

d 2 z dx 2 = M(x) EI(x) (4.15)
Where M(x) is the Bending-Moment function corresponding the loading condition, coefficient E represents the Elasticity Modulus of the beam material and I(x) is the Second moment of area of the cross section of the beam with respect to a Central axis perpendicular to the plane of the bending moment. This property (I) is presented as a function of length since our model of the beam does not possess a uniform cross section across the beam.

To find the best value for initial base of the triangular beam (b 0 ), we use the cross section of the actual geometry of the knife. The second moment of area in the cross section of both model and actual geometry should be equal. Figure 4.37 shows the cross section of the two mentioned geometries.

One should note that the cross section of the modeled geometry is a rectangle while the one of the actual geometry consists of two additional triangles on the sides of a rectangle. Since the thickness of the rectangle in both geometries are equal, in Considering the difference of the length of the two triangles to be b 0b. using the second moment of area for a rectangle with sides of b 0b and t with respect to the Central axis crossing the latter side, and the one of a isosceles triangle with a height of d and base of t with respect to the axis of symmetry we have:

(b 0 -b)t 3 12 = 2 (d)t 3 48 (4.16)
Solving for b 0 gives:

b 0 = b + d 2 (4.17)
Due to the dimensions of the knife, the base of the modeled triangular beam of the knife will be 12.75(mm).

Let's consider the following loading condition of the beam geometry: In a section of distance x of the base, the length of the cross section rectangle (b(x)) can be defined as:

b(x) = b 0 (1 - x h ) (4.18)
Thus the second moment of area as function of distance from the base I(x) can be expressed as:

I(x) = b 0 (1 - x h ) t 3 12 (4.19)
The resulting torque applied on the section M(x) is expressed by the following expression:

M(x) = F(h -x) (4.20)
Replacing M(x) and I(x) in equation 4.15 with the obtained results in equations 4.18 and 4.20, reveals the governing differential equation of knife deformation as following:

d 2 z dx 2 = F(h -x) E * b 0 (1 -x h ) t 3 12 = 12hF Eb 0 t 3 (4.21)
Boundary conditions of the problem induce the "no-deformation" and "no-slope" constrain to the problem which results in the following equation expressing the deformation of the beam as function of distance from the clamp support z(x):

z(x) = 6hF Eb 0 t 3 x 2 (4.22)
The numerical values of the constant parameters defined in equation 4.22 are defined is SI system in table 4.2. The value of parameter b 0 has been shown with a different color in table 4.2. The reason behind is the fact that this parameter has been calculated based on the specific characteristics of compliance behaviour of knife structure in vertical condition. The parameter b 0 will be used in further steps to find equivalent rectangle of the polygonal shape of the real knife's cross section. However, depending on the approach and direction of force imposition, the equivalent rectangle and consequently the value of b 0 would be subjected to change.

To evaluate the stiffness value of the knife in vertical direction (k zz ), we make use of the achieved knife deformation expression (equation 4.22) to find the relation between the vertical deformation of the knife tip (δz T ) and the applied vertical force f z :

k zz = f z δz T = Eb 0 t 3 6h 3 = 2.73 × 10 4 (N/m) (4.23)

Knife Compliance Behaviour in Lateral Direction

The compliance behaviour of the knife structure in lateral direction shown in Figure 4.39 is surveyed using same approach of the one described in vertical compliance behaviour. The deformation mechanism for both loading conditions is based on the bending moment imposed by the force along the body. Thus the governing differential for beam deformation is the following:

d 2 y dx 2 = M(x) EI l (x) (4.24)
Where I l (x) represents the second moment of area for lateral loading of the knife structure. However, the great difference in the two mentioned directions is due to the difference of the second moment of area about the central axis of cross section. As shown in Figure 4.40, for lateral banding, the tensile stress distributes along the vertical direction and thus the second moment of area about the z axis is involved. As shown in Figure 4.40, the area of the cross section consists of a rectangle and two triangles (shown as Tr in the Figure) which unlike the cross section related to vertical deformation of knife, the location of their center of area is not located on the symmetry line (axis z in this case). To achieve an equivalent rectangle having Chapter 4. Robot-Tool Compliance Treatment In Machining Process identical thickness t and length of b 0 , we make use of the same approach presented in previous section (vertical deformation).

The second moment of area for the actual cross section geometry is the summation of second moment of area of a rectangle having sides of t and b and two triangles with located away form the symmetry line with a distance of b 2 + d 3 . Thus the total value of second moment of area for the cross section is:

I z = 2 × dt 2 × ( b 2 + d 3 ) 2 + td 3 36 + tb 3 12 (4.25)
And the second moment of area for the equivalent rectangle would be:

I z(eq) = tb 3 0 12 (4.26)
By equating the second moment of area I z and I z(eq) , in equations 4.25 and 4.26, the length of the equivalent rectangle (b 0 ) for the lateral bending behaviour as a function of the geometrical parameters of the actual cross section is found as:

b 3 0 = b 3 + 2d 3 + bd × (3b + 4d) (4.27)
Using the numerical values of b and d for the initial cross section which are 11(mm) and 3.5(mm) respectively, gives the value of b 0 = 14.8(mm) for the equivalent rectangle.

By having the rectangular equivalent of the cross section of the knife structure for lateral deformation behaviour, the governing equation of knife behaviour defined in equation 4.24 can be rewritten as:

d 2 y dx 2 = f y (h -x) Eb 3 0 (1 -x h ) 3 t 12 = 12 f y h 3 tEb 3 0 × 1 (h -x) 2 (4.28)
Applying zero slope and displacement of the knife structure at x = 0 ( dy dx | x=0 = 0,y(0) = 0), the solution of the differential equation 4.28, is the following function:

y(x) = 12 f y h 2 tEb 3 0 × h ln h h -x -x (4.29)
The displacement function y(x) suggests an unlimited displacement for the knife tip (x = h). This is due to the fact that the second derivative of the displacement expressed in equation 4.28 is undefined. One should note that the Bernoulli beam theory is valid for the conditions where deformation of neutral axis is small. On the other hand, the actual geometry of the knife presented in Figure 4.34 states that the knife does not have a pointy tip. Thus to evaluate the deformation of the knife tip in lateral direction, we make use of the value of the displacement resulted from equation 4.29 at x = 0.98 × h. In this case, the lateral stiffness of the knife structure turns to be:

k yy = f y δy T = tEb 3 0 12h 3 × 1 ln(50) -0.98 = 2.95 × 10 5 (N/m) (4.30)

Knife Compliance Behaviour in Longitudinal Direction

The compliance behaviour of knife structure in longitudinal direction, unlike the lateral and vertical directions, can not be studied by the bending beam theory of Bernoulli. In the loading condition applied on this direction, the force does not impose any moment on any principal axes of the knife geometry. As shown in Figure 4.41, the load applied on the longitudinal direction f x , causes a displacement of δx T in the knife tip. To evaluate the displacement of the knife structure due to the longitudinal force ( f x ), we make use of Hooke's law by assuming that the axial stress imposed along the knife structure (σ(x)) will retain the knife material in the elastic Stress-Strain zone.

By recalling the Hooke's law for elastic zone of Stress-Strain, the strain along a material obeys the following expression:

= σ E (4.31)
Where is the axial strain along the structure, σ is the axial stress and E as defined earlier, is the modulus of elasticity of the material. The deformation of the material possessing a length of L (∆L), undergoing a constant strain along its axis, can then be described as:

∆L = × L (4.32)
To model the compliance behavior of the knife structure under the longitudinal force based on the Hooke's law, one should note that due to the varying area of the knife's cross section, the axial stress distributed along the longitudinal axis varies due to the definition of the axial stress:

σ(x) = f x A(x) (4.33)
In which A(x) is the area of the cross section as function of the distance from knife base (x).

To evaluate the deformation of knife followed by a longitudinal load ( f x ), let's consider an element along the longitudinal axis x with a side value of dx as shown in Figure 4.42. Due to the negligible variations of stress (σ(x)) along the small element with the length of dx, one can use the equation 4.32 to express the deformation of the element. Thus the derivative of the deformation of the knife (dL) can be expressed as:

dL = × dx (4.34)
By replacing in equation 4.34 with its equivalent in equation 4.31, the derivative of the deformation can be rewritten in the following format:

dL = f x EA(x) dx (4.35)
The displacement of the knife tip along the longitudinal direction, would be trivially the summation of the derivative of deformation of the elements from the knife base to the tip:

δx T = dL = h 0 f x EA(x) dx (4.36)
The actual shape of the cross section of the knife (A(x)) is the polygonal shape presented in Figures 4.40 and 4.37. To ease the procedure of evaluating the integral presented in equation 4.36, the cross section of the real shape will be modeled by a rectangle with equivalent area. Thus the length of equivalent rectangle (b 0 shown in Figure 4.40) is b 0 = b + d. Therefore the following function expresses the area of cross section as a function of x:

A(x) = b 0 t(1 - x h ) (4.37)
In which the numerical value of b 0 equals to 14.5(mm) (b 0 = b + d). By having the function of area of the cross section as stated in equation 4.37, The solution of the integral becomes:

δx T = lim x→h h f x Eb 0 t ln h h -x (4.38)
Since the resulted expression also contains a logarithmic term approaching to infinity such as the one of lateral deformation of knife structure, the compliance of the knife will be evaluated using a point in vicinity of the knife tip (x = 0.98 × h):

k xx = f x δx T = Eb 0 t h × 1 ln(50) = 5.82 × 10 7 (N/m) (4.39)
To evaluate the knife compliance behaviour in a practical condition, we make use of the chamfering forces of hard honeycomb material for a chamfering angle of 30 [deg] as adopted for the robot compliance simulation. The machining forces defined in the knife frame are estimated using the force model (equation 4.5). Recalling the SSM principal equation for the knife structure stated in equation 4.14, the knife tip displacement can be evaluated using the following equation:

  δx T δy T δz T   =   k -1 xx 0 0 0 k -1 yy 0 0 0 k -1 zz   ×   -32.8(N) 108.0(N) 177.5(N)   (4.40)
Using the achieved stiffness values for different directions of the knife, the displacement array for the described loading condition is:

  δx T δy T δz T   =   -5.6 × 10 -4 0.37 6.50   (mm) (4.41)
As one can easily conclude by comparing different elements of the resulted displacement of the knife tip, the compliance behaviour of the knife in z direction (vertical direction of knife structure) is by far greater than the rest. The displacement in the longitudinal direction is less than a micron. The displacement of the lateral direction is about 20 times less than vertical direction. One should note that not only the stiffness values of these two directions are considerably larger than the vertical direction, but also the chamfering forces in these directions are always smaller than the vertical component.

Experimental Verification of Knife compliance Behaviour

According to the conclusion resulted from the analysis executed on modeling the compliance behaviour of knife structure in different directions, the vertical direction of the knife has shown the most compliant behaviour among all. Due to the importance of the knife compliance behaviour for further implementations, we decided to put it into examination to verify the validity of the generated model in comparison to an actual loading condition.

To ensure that the testing condition is close to the actual working condition, a long knife was attached to a ultrasonic head connector. The connector was fixed on a rigid metal jaw in the workshop. Load was applied using a rope and pulley system. different dead-weights were hanged on the free end of the rope while the other end was tied to a Rilsan necklace whose buckle was placed on the narrow zone of knife tip to ensure a concentrated loading condition. Test setup is shown in Figure 4.43. Two displacement gauges were used to measure the deflection of the knife in two different points. Gauge number 1 is the one measuring the deflection of the tip and Gauge number 2 was mounted to measure at 31(mm) above the root of the knife. Three different weights were attached to the rope and the deflections were captured using the gauges. The mass of the three weights were: 3.1, 6.2 and 8.7(kg). For each loading condition, the deformation of the knife structure in vertical direction is estimated using equation 4.22 for two values of x (x = 31(mm), 70(mm)). Figure 4.44 illustrates the model estimated deflection value along with the experimental data. As shown in Figure 4.44, the estimated values achieved from the generated compliance model of the knife structure in vertical direction complies with closely to 4.12. Robot-Knife Compliance Behaviour 103 the data driven from the experiments. This results confirms that the model can be reliably wielded for further analysis of the robot-knife compliance behaviour.

Robot-Knife Compliance Behaviour

According to the presented modeling of knife compliance behaviour in section 4.11, we concluded that the compliance behaviour of the machining tool in ultrasonic honeycomb machining process can not be neglected on the overall compliance behaviour. In this section, we aim to verify the contribution of both machining tool and robot structure to the observed machining errors in honeycomb chamfering process reported in section 4.6.

The approach to verify the effect of each system (knife and robot structure), is to use a virtual joint with equivalent rotational stiffness to represent the compliance effect of the tool attached to the robot. Thus compliance modeling of the resulted 7axes robot structure would be a generalized format of the one applied on the 6-axis robot.

To proceed further, the most compliant direction of knife structure, which is the vertical axis, is modeled by a rotational spring as shown in Figure 4.45. To achieve the equivalent stiffness value k eq , we make use of the displacement value δz T calculated in section 4.11.1. Since we are using rotational springs in our compliance model, the relationship between the applied force f z and tip displacement δz T , must be inverted to a torque-rotational deviation format. The lever of force f z on the virtual joint is the knife length h and assuming small value of tip displacement δz T , the angular deviation of joint, termed as δθ eq , is equal to:

δθ eq = δz T h (4.42)
Thus the equivalent torque-rotational deviation of the joint can be expressed as:

k eq × δθ eq = f z × h (4.43)
By replacing δθ eq in equation 4.43 with its equivalent achieved in equation 4.42, and using the relationship between δz T and vertical force applied on the knife tip f z achieved in equation , equation 4.43 can be rewritten by the following form:

k eq × ( 6h 3 f z Eb 0 t 3 )/h = f z × h (4.44)
Solving for k eq in equation 4.44 leads to:

k eq = Eb 0 t 3 6h (4.45)
By having the equivalent torsional stiffness of the machining tool, which is now considered as the 7th link of robot, the compliance behaviour of the robot-tool structure will be modeled using the same approach of an ordinary 6-axis robot. Since the compliance behaviour of robot structure is configuration dependent, the compliance behaviour of robot-tool structure will be surveyed in a configuration in which maximum torque is imposed on the robot joint (apart from joint 1).

Let's consider a configuration where the robot is completely bent on a horizontal line (axis x of robot base frame) and a vertical force ( f z ) is applied on the TCP (knife tip). to maximize the torque applied on joint 5, this joint is rotated for 45[deg] to ensure a vertical direction of the force applied perpendicular to the knife. Figure 4.46 illustrates the modeled robot-tool structure undergoing described force imposition condition. The displacement of the TCP, as mentioned previously in chapter 3, can be calculated by the following expression:

δP t ω t = (J 6×7 × [K] -1 × J T 7×6 ) × F (4.46)
In which J and [K] are Jacobian and stiffness matrices and F is the external couple array with the following entries:

F =         0 0 -f z 0 0 0         (4.47)
The stiffness matrix stated in equation 4.46, is a 7 × 7 matrix containing 6 joint stiffness values attributed to the robot joints and an additional equivalent rotational stiffness representing the torsional stiffness of the machining tool:

[K] =           k 1 0 0 0 0 0 0 0 k 2 0 0 0 0 0 0 0 k 3 0 0 0 0 0 0 0 k 4 0 0 0 0 0 0 0 k 5 0 0 0 0 0 0 0 k 6 0 0 0 0 0 0 0 k eq           (4.48)
The Jacobian matrix (J) is computed for the configuration presented in Figure 4.46 using the following joint variables (defined in degrees):

θ 1 = 0 θ 2 = 90 θ 3 = 0 θ 4 = 0 θ 5 = -45 θ 6 = 0 θ 7 = 0
The value of displacement δz tcp , the third element of displacement vector δP t , for the given configuration and as function of the joint stiffness values becomes:

δz tcp f z = 561.3 × k -1 2 + 201.4 × k -1 3 + 38.3 × k -1 5 + 0.49 × k -1 eq (4.49)
To evaluate the portion of contribution of each element in robot-tool compliance behaviour, the stiffness parameters appeared in equation 4.49, are replaced with their corresponding numerical values. The robot joint stiffness parameters are replaced with the values presented in table 3.7 and for the equivalent torsional stiffness of knife (k eq ), the numerical value achieved by replacing knife properties in equation 4.45 is adopted:

k -1 2 = 2.3810 × 10 -7 (rad/N.m) k -1 3 = 1.2618 × 10 -6 (rad/N.m) k -1 5 = 4.1152 × 10 -6 (rad/N.m) k -1 eq = 7.5 × 10 -3 (rad/N.m)
Replacing the above presented numerical values in the compliance robot-knife compliance component in z direction achieved in equation 4.49, reveals the following percentage of contribution for robot structure and knife element:

δz tcp f z = 561.3 × k -1 2 + 201.4 × k -1 3 + 38.3 × k -1 5 11% + 0.49 × k -1 eq 89% (4.50)
The results presented by equation 4.50 indicate that the effect of compliance behaviour of knife structure is dominant to the portion raising from robot structure by more than 8 times. These results have been achieved in a configuration in which Chapter 4. Robot-Tool Compliance Treatment In Machining Process the compliance of robot in vertical direction is considerably higher than the configurations we face for an ordinary machining task executed by robot arm. Therefor, the contribution of the knife structure is expected to be even more for the real cutting condition surveyed in section 4.10.

Due to the presented importance of the knife structure in the compliance behaviour of the machining system, in next step we try to relax the assumption on the loading condition on the knife behaviour to achieve a more accurate knife compliance behaviour.

Modeling of the Loading Condition

In the previous steps, we assumed that the vertical component of the machining force is applied as a concentrated force on the point of the knife. Hence, the mechanism of the machining process does not fully correlate to a concentrated single force applied on the tip. In this part we propose a more general regime of loading to observe if the simulated behavior of the knife would describe more the practical results. The proposed general regime follows the consists of the following assumptions:

• Instead a single force vector, or combination of several concentrated forces, this model is based on the distributed loading condition which seems more realistic due to the physics of the problem.

• The distributed load subjects the penetrating part of the knife inside material. This length is termed as h 0 hereafter and is computed using the following expression:

h 0 = T -l sin(β) (4.51)
Where T is the thickness of the honeycomb workpiece (T = 25(mm)) and l is the protection gap as presented in Figure 4.10.

• The loading obeys an increasing trend along the knife length. In other words, the distributed loading rises along the affected length(h 0 ) and meets its maximum value at the ending point of the knife. With this assumption, the model should still be capable of describing the knife deformation under a concentrated force at the tip as an extreme condition.

Among all the mathematical functions for describing the distribution condition that can satisfy all the above mentioned assumptions, choosing a general polynomial function would be a wiser option due to their vast usage and the simplicity of the algebraic operations. The distributed load function ω(x) of order n can be expressed with the following general expression:

ω(x) = ω 0 xn (4.52)
Where x is the dimensionless length variable with the following definition:

x = x -(h -h 0 ) h 0 (4.53)
And ω 0 is the maximum value of the distributed load met at the ending point of the knife. As the resultant force of the distributed load should be equal to the measured force (F) By adjusting a certain value to n, the value of ω 0 can be calculated from the following integral:

F = h h-h 0 ω(x)dx (4.54)
By changing the variable x to x in equation 4.54 and applying the effects on the limits of the integral and using the general expression explained in equation 4.52, we have:

h h-h 0 ω(x)dx = 1 0 ω 0 h 0 xn d x (4.55)
Leading to the relation between ω 0 , n and F: Displacement and deflection of the knife structure are computed using the governing equation. Since the loading is distributed only over one part of the knife, the torque expression is not unique all along the knife. Knife length can be separated into two parts, To proceed further, hereafter we use Left Part to refer to the zone of the knife which in not subjected to the external load (hh 0 > x > 0). Right Part on the other hand presents the loaded zone of the knife (h > x > hh 0 ).

ω 0 = (n + 1) F h 0 (4.56)
The bending-moment of the left part deals only with the placement of the equivalent force of the loading. Magnitude of the equivalent force is F and it is located on the centroid of the trapped area under a polynomial of order n. It can be easily shown that the x coordinate of the centroid of the area under a polynomial with general form of y(x) = ax n in the interval of x 1 > x > 0 is equal to x c = n+1 n+2 x 1 . Setting x 1 = h 0 and converting the origin of x from the beginning point of the external load to the fixed support, we have the below final format describing the bending-moment of the left part:

M 1 (x) = F * (h - h 0 n + 2 -x) x ∈ [0, h -h 0 ] (4.57)
For the right part, the bending-moment expression (M 2 (x)) was calculated by applying the equilibrium equation over the right part of a section at x position. Consider the right part of the cross section in Figure 4.48. To find the expression of M 2 (x), we use the following general form:

M 2 (x) = h x X × ω(X)dX (4.58) M 2 (x) = F n + 2 (h 0 + h -x) n+2 h n+1 0 -h 0 + (n + 2)(h + x) x ∈ [h -h 0 , h] (4.59)
Recalling the governing equation of the knife deformation:

d 2 z dx 2 = M(x) EI(x) (4.60)
Where E is the Young Modulus of the material and I(x) is the value of second moment of inertial of a cross-section at pint x in the vertical direction. We should note that in some experiments, we have noticed that the geometry of the cutting knives differed slightly case by case specially in the sensitive property of consistency of the knife thickness. In some cases, the variation of the thickness of knife structure (which is 2(mm) in the nominal value) could achieve 0.6(mm). In cases of such, the thickness has been modeled by a linear function of distance from the knife base (x) and integrated in the value of I(x). Due to the comparatively complex form of final format of this integral, we have adopted Simpson's rule to numerically evaluate the integral presented in equation 4.60 using 200 points in the domain.

To find the best value of n describing the experimental error values, we went through a minimization method. Er(n) is the error function of comparing the modelresulted of chamfering error with the experimental values for a certain n value:

Er(n) = ∑ N i=1 E (exp) (β) -E(β, f z (β)) 2 N (4.61)
In which N is the number of experimental data and E i(MP) is the slope of the knife tip (x = h) estimated by described compliance model of knife structure using distributed load of order n for machining feature corresponding to point i and E i(exp) is the experimentally evaluated chamfering error for the identical chamfering angle (β) value of point-index i. Figure 4.49 depicts the value of the described error function for wide range of n = 0 to n = 100 with the increments of 0.2. [START_REF] Krastev | Passing Through Jacobian Singularities in Motion Path Control of Redundant Robot Arms[END_REF].49 illustrates the error of the estimating model of the compliance behaviour of knife structure for a wide range. This range covers the complete uniformly distributed load (n = 0) to the ranges where the force can be practically assumed as a concentrated force (n = 100). The behaviour of the graphs for both criteria follows the same trend. By marching along the n axis in positive direction, the estimated error decrease to a certain point and then increases. This clearly indicates that the assumption of a distributed load instead of a concentrated force is valid.

Chapter 4. Robot-Tool Compliance Treatment In Machining Process

According to Figure 4.49, n = 7.8 best describes the chamfering angle error, while in terms of geometric error D1, this number changes to n = 16.4. A simple comparison between the error value (Er(n)) of the two mentioned n values in D1 criterion, declares that the difference does not exceed 0.1(mm) which is way lower than the precision of the measurement for D1. Thus hereafter, we use the value n = 7.8 for further steps. Figure 4.50 shows the predicted error values for several loading conditions (n) along with the experimental values. A noteworthy conclusion can be drawn form Figure 4.50 regarding the estimated deviation error using different values of load distribution order n. The estimated deviation error increases by increasing the value of n. According to the data presented in this Figure, the variation of this parameter can affect the estimated error drastically while the proportion of the highest and lowest estimated errors corresponding to the two extremities of this parameter differ up to 3 times.

Compensation Procedure

In previous sections, modeling of the knife structure was surveyed and experimentally validated. the loading condition of the machining operation on the knife was carried out using a non-uniformly distributed force. In this part, a compensation method is presented to overcome the geometric error in honeycomb cutting process with ultrasonic technology using triangular knife.

The compensation procedure proposed in this part is based on offline approach. The strategy wields the developed force model in knife's vertical direction to estimate the machining forces applied on the knife and consequently uses the compliance behaviour model generated previously to estimate the chamfering error based on the estimated forces for the required chamfering features. The process aims to propose a modified chamfering angle (β m ) for a corresponding desired chamfering angle of β. A detailed explanation of the proposed compensation strategy comes below.

More precisely, the logic behind the method is to achieve the answer of the following question: for a desired chamfering value of β, what is the corresponding chamfering angle β m , for which by programming the chamfering operation by this value, an actual chamfer with β will be gained?

Computing the value of β m corresponding to the desired value of β consists of an iterative algorithm. In the first step a chamfering angle is proposed (for the first iteration, desired value β is used). The cutting forces for the proposed chamfering value are estimated using the vertical machining force model. Deformation of the knife structure due to the estimated force is computed using the knife model which is termed as θ dev hereafter. This deformation can be estimated using either concentrated or distributed load condition. The estimated chamfering angle (β = β m + θ dev ) is then compared to the value of desired chamfering angle β. The chamfering error of each iteration er is defined as er = ββ . If the computed error of iteration (er) is higher than the predefined value of , the process continues with the updated value of β m computed as β m = β mer. Loop ends when the iteration error value falls lower than the value of . f Vertical (β m )

Knife Compliance Behaviour in Vertical Direction

Concentrated Force (Equation 4.22) Distributed load (n = 7.8) (Equation 4.59) In light of the explained procedure of the proposed compensation method and the flowchart presented in figure 4.51, it is evident that the elements involved in the compensation strategy such as force model and tool compliance behaviour are dependent to to the specific characters of a machining process. In other words, although the outline of the strategy is generic and can be adopted for a wide range of applications, but on the other hand, for each set of (tool, workpiece material, machining process), related blocks presented in Figure 4.51 must be developed based on the pre-executed experiments or simulations. using specific features of our case study, for hard honeycomb material, the offline compensation strategy was applied using developed vertical force model and knife compliance behaviour on a practical range of chamfering angle by setting the stopping criteria = 0.01 [deg]. Two graphs are presented in Figure 4.52. The difference between the two graphs are the load conditions applied on the knife (concentrated force assumption for n = ∞ and distributed load with the optimal distribution order n = 7.8 found in Figure 4.49). The geometry of the knife used for the graph having n = 7.8 is the modified geometry regarding the knife thickness as mentioned earlier. The trend of the both graphs indicates that the difference of the modified chamfering angle (β m ) and the desired chamfering angle (β) increases while decreasing the desired chamfering angle. This behaviour on the graph can be reasoned by the augmenting trend of the vertical machining force while chamfering angle increases which consequently increases the deviation of knife structure involved in the chamfering process.

θ dev er = β -β m -θ dev β m Accepted β m er > β m = β m -er
We should also note that since the distributed load assumption predicts lower chamfering error for a given vertical machining force compared to the one of concentrated force assumption, for an identical knife structure, the modified chamfering angle β m should have been closer to the value of the desired chamfering angle β. But the graph corresponding to distributed load presented in Figure 4.52 presents a contrasting result. This is due to the fact that the knife structure used for this load compensation, possesses a linearly decreasing thickness from 2.1(mm) in the base to 1.5(mm) in the tip. We believe that the additional compliance of knife due to the slight changes in the knife thickness is responsible for this behaviour.

To evaluate the efficiency of the proposed offline compensation procedure, experiments were conducted on the hard honeycomb material. Before presenting the results we should note that three different experimental sessions were executed. The presented data for machining force modeling presented in section 4.9 and machining errors brought in section 4.8 were resulted from Session 1. Two other sessions were also executed to verify the efficiency of the compensation procedure which here after, will be referred to as Session 2 and Session 3 depending on the time order of execution. In second experimental session, only the compensation procedure based on concentrated force model was verified. However on the third and last experimental session, both assumptions (concentrated and distributed load with n = 7.8 shown in Figure 4.52) were examined. Figure 4.53 presents the chamfering error in honeycomb cutting process on the hard material using compensation strategy along with the uncompensated error measured in the first session. An important notable point based on the data depicted in Figure 4.53 concerns the repeatability of the process. The machining features used for the evaluation of the compensation strategy based concentrated force assumption executed on sessions 2 and 3 were identical. But however, the resulted chamfering angular errors do not comply. We attribute the observed noncompliance to the non repeatability of this machining process. Although the subjected honeycomb material and machining features such as programmed chamfering angles and feed rate stayed unchanged in all three mentioned experimental sessions, several elements of the system were unavoidably changed such as robot programming software. It is also worth mentioning that the honeycomb material boosts its non-homogeneous characters which trivially affects machining forces while sensitive parameters such as the cells arrangement with respect to cutting knife change. To depict the non repeatable nature of the process, the vertical chamfering force applied on the knife and chamfering errors for uncompensated machining condition is presented in Figure 4.54. .54 clearly indicates that the machining forces in honeycomb chamfering process using ultrasonic technology was not repeatable during the different experimental sessions. We believe that repeatability of this important parameter can not be ensured unless a complete control over the various parameters (such as programming software, workpiece placement with respect to machining tool, cutting knife which comes with its geometrical manufacturing errors case by case, configuration of robotic arm and methodology of the operators designing the process) is in hand. Further investigations on this process and error improvement is out of the scope of this work. However, by comparing the chamfering errors resulted from different sessions with their corresponding vertical machining force, an evident similarity between the trend of both properties is vivid which can be served as a final evidence of the hypothesis of dominant importance of this force component in the observed geometrical manufacturing errors in honeycomb workpieces in ultrasonic chamfering process.

Conclusion

In this chapter, we proceeded to survey a robotic-based machining process in which the compliance behaviour of the machining tool could not be neglected compared to the one of the operating industrial robotic arm in the total observed manufacturing errors. The machining process adopted for this chapter was the chamfering process of Aramid honeycomb materials using ultrasonic technology based on the requirements and interests of Le Creneau company.

Firstly, a survey of recent research works was presented stating the stat of the art in the topic of robotic machining. Further on, a literature review was brought indicating the advantages of the honeycomb structure and materials for the engineering applications followed by the notable points regarding the challenges of manufacturing and simulation of the behaviour of materials of such type.

To proceed with our investigations on finding the main contributor to the machining errors in cutting honeycomb materials with ultrasonic cutting technology wielding a robot arm as the operator and triangular knife as the cutting tool, we were in need of the machining forces involved in the process. Due to the lack of an accurate machining force model in the literature, an empirical approach was used to generate a force model by executing experimental tests.

The generation of the machining force model was based on the variable control method and capturing experimental values of the exact machining forces. Data were then used in a curve-fitting process leading to expressions able to estimate machining forces in different directions defined in both tool and work-piece frame for a given chamfering angle for our two main focused honeycomb materials.

The chapter presented a novel strategy in measuring geometrical errors observed in the resulted workpieces by categorizing them into two main distance and angle errors. Due to the shortcomings and practical limits, the measurement procedure involved different measurement devices and on-site challenges. Our proposed solutions were presented to overcome the inconveniences.

In our first attempt, we simulated the compliance effect of robot structure. The results of the simulations on the robot compliance behaviour were inaccurate and incapable of explaining the measured error values in experimentally-resulted workpieces. The order of misalignment and displacement error in TCP orientation and positioning were negligible, therefore, the compliance behaviour of the cutting knife was integrated.

To integrate the compliance model of the cutting knife, we made use of the Structural Stiffness Matrix (SSM) approach. The stiffness matrix of the knife structure relating the applied force on knife tip to its displacement was then evaluated in three main directions of knife structure using either Euler-Bernoulli beam theory or Hooke's low depending on the nature of the deformation on each direction. the analysis revealed that the vertical direction of the knife structure is by far the dominant direction in terms of compliance behaviour. Therefore, the effect of the other directions was neglected. The compliance model of the knife was put to experiments in which the results proved its efficiency.

The integration of the compliance behaviour of the knife in the machining process was based on a novel approach proposed in this chapter. The logic behind the approach was simply to introduce the machining tool as an additional virtual joint and link. The virtual joint was chosen to be in the direction of the most compliant direction of the tool structure. According to this approach, the compliance behaviour Chapter 4. Robot-Tool Compliance Treatment In Machining Process of the entire robot-tool system was modeled using the equivalent 7 axes robot arm in which the corresponding stiffness value of the additional virtual joint was formulated using the equation governing the knife behaviour in the vertical direction.

Although due to the properties of our machining setup, the portion of contribution of the robot compliance behaviour was negligible compared to the one of tool structure, the proposed method stays with its generic form to serve for applications where compliance of both systems are in the same order.

In a further step, for more accurately modeling of the knife behaviour, the load applied on the knife side plane was modeled by a non-uniformly distributed force. The non-uniformly distributed force was expressed by a polynomial of order n, whose value was determined by applying the least square method on the experimental data.

the chapter continued by proposing a compensation method. The method was based on iterative computational procedure able to predict the initial chamfering value corresponding to a desired chamfering value taking into account the force model the consequently the deformation of the knife structure due to its corresponding machining force. The outcome of this compensation method has shown promising results.

Chapter 5

Optimal Workpiece Placement In Robotic Cells

Introduction

The percentage of the robot-involved machining process is augmenting rapid ally in different industrial sectors. The increase of demands for robot manipulators year after year is evident for the advantages of industrial robot arms in many applications [START_REF] Verl | Robots in machining[END_REF]. Industrial robots are known for their comparatively lower prices compared to the conventional apparatus while offering appealing advantages such as widen workspaces, reconfigurability, and dexterity. Nowadays, robot manipulators are playing roles in different machining processes such as milling, grinding, deburring and drilling to name a few [START_REF] Ji | Industrial robotic machining: a review[END_REF].

However, the number of disadvantages of robot manipulators is preventing their higher rate of employment and contribution to the precision manufacturing process. Inherent shortcomings of the robotic arms are among the main contributors to machining errors and inaccuracies. Low joint stiffness values and serial structures of industrial robotic arms introduce great challenges such as considerably higher order of compliance behaviour compared to the conventional machining devices, dynamic misbehavior responsible for undesired chatters [START_REF] Huynh | Modal updating of a 6-axis robot for milling application[END_REF], vibrations and increased machining loads.

Forgoing points explain the sharply increasing number of research items in the literature and researchers focusing their efforts on increasing robot accuracy to further improve their performance by reducing the deficiencies sourced from various error sources. Vastly different solutions have been proposed which differ in the approach and level of complexity based on specific requirements. These solutions span from simple off-line compensation methods for a defined machining load and conditions to costly error reduction strategies which usually involve additional sensors such as external encoders [START_REF] Klimchik | Robotic manipulators with double encoders: accuracy improvement based on advanced stiffness modeling and intelligent control[END_REF], piezoelectric [START_REF] Denkena | Adaptronic Systems in Robot Manufacturing[END_REF] and laser tracker [START_REF] Schneider | Improving robotic machining accuracy through experimental error investigation and modular compensation[END_REF]. [START_REF] Nelson | Locating assembly tasks in a manipulator's workspace[END_REF] opened a new chapter in the course of robotic research and pioneered a way which was a starting point for many upcoming researchers. To achieve a batter performance of a robotic application, they proposed to modify the configuration of the robot cell. The so-called workcell design expanded to optimize many performance indices of a robotic workcell such as cycle time [START_REF] Spensieri | Optimal robot placement for tasks execution[END_REF] and energy consumption (Bukata and Hanzálek, 2019) by further efforts of experts.

Chapter 5. Optimal Workpiece Placement In Robotic Cells

Workcell design approach entered the robot machining error reduction problem as a subclass of offline strategies. Researchers started to benefit from either the redundancy introduced by 6 axis robot arms dedicated to 5 axis machining tasks [START_REF] Xiong | Stiffness-based pose optimization of an industrial robot for five-axis milling[END_REF] or configuration-dependent properties of the robot such as local stiffness or eigenfrequencies [START_REF] Celikag | Chatter Suppression in Robotic Milling by Control of Configuration Dependent Dynamics[END_REF][START_REF] Chen | Posture-dependent stability prediction of a milling industrial robot based on inverse distance weighted method: Paper ID: 1222[END_REF][START_REF] Denkena | Design and optimization of a machining robot[END_REF] to improve machining accuracy.

According to the literature, one of the most efficient approaches to increase robot performance in a machining task is to modify the relative placement of the workpiece with respect to the robotic arm. According to [START_REF] Caro | Workpiece placement optimization for machining operations with a KUKA KR270-2 robot[END_REF], for a Kuka IR executing a milling operation on an aluminum block, the manufacturing error resulted from robot compliance behaviour is reduced about 80% in the optimal placement compared to the worst placement.

Optimizing the relative workpiece-robot placement is not limited only to workpiece displacement. In some applications such as arc welding over giant pipelines optimizing the base placement of the robot arm is more convenient [START_REF] Doan | Optimal robot placement with consideration of redundancy problem for wrist-partitioned 6R articulated robots[END_REF]. In applications where the robot's end-effector can not reach entire points on the workpiece with appropriate orientation with stationary location, researchers have turned their views to employ mobile platforms and consequently search for optimal paths and configurations for highly redundant systems of such type [START_REF] Subrin | Digital chain development for sanding application with a kinematically redundant robotic system[END_REF].

Problem Statement

In this section we tend to define different aspects of the statement of the concerned problem. Different steps of the procedure including problem definition and proposed solution are thoroughly reported in further sections.

Our effort in this work is dedicated to improve the quality of a trimming process as an ending machining process affected by the joint transmission periodic errors causing cyclic in the displacement of the robot end-effector and robot compliance behaviour. These errors are mainly cased by manufacturing and assembly errors of the gearbox resulting in geometrical imperfections of the transmission system from its desired state.

Several remedies have been offered in the literature to counter the joint transmission periodic effects [START_REF] Vakilinejad | Identification and Compensation of periodic gear transmission errors in Robot Manipulators[END_REF][START_REF] Olabi | Improving the accuracy of industrial robots by offline compensation of joints errors[END_REF]. But however, they are mainly based on using the mathematical error model to eliminate the error portion directly from the desired joint values. These solutions require postprocessing after each placement of the robot and probably command authority in robot controller level which may not be feasible or economically efficient for many applications.

In this work, by assuming that the mathematical model of the joint errors for the contributing joints as well as joint stiffness values are in hand, we propose an optimization procedure to reduce the undesired effect of these non geometrical sources on the machining quality by means of replacing the the workpiece in a close optimal position in robot workspace.

The workpiece adopted for this work consists of a circular path representing trimming trajectory of a circle with a diameter of 25(cm) made of Carbon fiber reinforced plastic (CFRP) due to their lower manufacturing tolerances required for these materials compered to their equivalent metal parts [START_REF] Denkena | Design and optimization of a machining robot[END_REF]. A Stäubli TX200 industrial robot arm equipped with a spindle is chosen to execute the task (Figure 5.1). As depicted in Figure 5.1, placement of the workpiece is fully defined using four variables. Three variables stand for Cartesian coordinate of the center of the workpiece (x c , y c , z c ) and the fourth variable represents the orientation of the circle plane along the vertical axis of universal plane. To sake simplicity, the universal frame is identical to robot base frame which is constructed by the following procedure: Z axis of the coordinate system is collinear to the axis of joint 1 (axis z 1 depicted in Figure 5.1) and x axis is the common normal of axes of joint 1 and 2 of the robot arm. The TCP frame with respect to the robot flange is also given as:

T TCP Flange =     cos( -3π 4 ) 0 sin( -3π 4 ) 300 0 1 0 0 -sin( -3π 4 ) 0 cos( -3π 4 ) 300 0 0 0 1    
(5.1)

Joint transmission Error

As mentioned previously in chapter 3, nonlinear joint transmission system errors are originated from the wear effect and eccentricity of the input and output shaft in transmission system of a joint. These imperfections cause periodic errors in the joint variable superimposed on the desired joint variables which results in flactuational displacement of TCP along a given path.

To integrate the effect of periodic errors of joint transmission systems, a mathematical model capable of estimating this error for a given desired joint values is needed. In [START_REF] Vakilinejad | Identification and Compensation of periodic gear transmission errors in Robot Manipulators[END_REF], we have proposed a method based on an optimal path generated in Cartesian space to identify and mathematically model the nonlinear transmission joint errors. Due to lack of experimental

k θ =    k θ1 • • • 0 . . . . . . . . . 0 • • • k θ6    (5.4)
And k c is the complementary stiffness matrix for a loaded robot situation introduced in [START_REF] Chen | Conservative congruence transformation for joint and Cartesian stiffness matrices of robotic hands and fingers[END_REF] This term is neglected in this work due to its minor effect as claimed in [START_REF] Xiong | Stiffness-based pose optimization of an industrial robot for five-axis milling[END_REF] The joint stiffness values used in this work are the ones achieved during elastostatic model generation stated in table 3.7.

Modeling of the forces involved in machining of CFRPs has been subjected in a vast number of research works lately. Researchers are putting efforts to generate models capable of estimating forces in different machining operations such as milling and cutting of these materials by taking into considerations different machining parameters like feed rate, tool geometry [START_REF] Chen | Study on cutting performance of diamond-coated rhombic milling cutter in machining carbon fiber composites[END_REF] and workpiece thickness [START_REF] Boudelier | Cutting force model for machining of CFRP laminatewith diamond abrasive cutter[END_REF].

According to th literature of this topic, the machining forces are mainly distributed in along in-plane normal and feed rate directions and the force component along the tool axis is neglected. This fact can be concluded from the experimental results presented in (Sheikh-Ahmad, [START_REF] Sheikh-Ahmad | Cutting force prediction in milling CFRPs with complex cutter geometries[END_REF][START_REF] Slamani | Influence of machining parameters on surface quality during high speed edge trimming of carbon fiber reinforced polymers[END_REF]. Therefore, this force component is eliminated from the machining interaction model in this work. Figure 5.3 depicts the loading condition of the CIRP trimming adopted in this work. Machining forces defined in the supporting table frame are as following:

-→ F st =   0 -F r cos(α) -F t sin(α) F r sin(α) -F t cos(α)   (5.5)
To transfer the force vector from supporting table frame to the robot base frame, we make use of the rotation matrix of θ about z axis:

-→ F = Rot z (θ) × -→ F st (5.6)
Which will lead to:

-→ F =   sin(θ) × (F r cos(α) + F t sin(α)) -cos(θ) × (F r cos(α) + F t sin(α)) F r sin(α) -F t cos(α)   (5.7)
As mentioned by Karpat, Bahtiyar, and De ger, 2012, the tangential and radial force components (F t and F r respectively) for carbon fiber reinforced polymers (CFRP) materials are dependent on different machining parameters. They have surveyed the effect of feed rate, tool rotation and fiber direction of the workpiece on the trimming machining forces. According to their experimentally driven values, the maximum tangential force component during one rotation of the cutting tool, varies between 100(N) to 140(N) for a feed rate of 350(mm/min) as a function of fiber direction. These values for the radial force component are 180(N) to 300(N). In this work, the maximum values for both components (F t = 140(N),F r = 300(N)) are employed to represent the trimming forces.

Cost Function Definition

In robot machining operations tool path deviations are not limited to the feed direction. Even in the processes where the load is mainly in the machining direction the tool may expose high value of deviations in side directions due to inequality of local stiffness in different directions. Furthermore, in considerable number of processes such as deburring and trimming the placement errors orthogonal to the feed direction are main contributors to the machining imperfections. Machining defeats for a milling operation shown in [START_REF] Ozkana | Milling Behavior Analysis of Carbon Fiber-Reinforced Polymer (CFRP) Composites[END_REF][START_REF] Yaşar | Experimental investigation on novel drilling strategy of CFRP laminates using variable feed rate[END_REF][START_REF] Yaşar | Experimental investigation on novel drilling strategy of CFRP laminates using variable feed rate[END_REF] are examples of this claim.

In this work we propose a cost function which is designed to minimize the sidedirectional machining errors independent to the concentrated error source. Figure 5.4 shows a schematic of a desired and an actual path resulted from robot related errors. Let us consider a group of two successive points used to describe the desired machining path P i and P i+1 shown in figure. Due to the described error sources, tool path is deviated so that P i and P i+1 are displaced to P i and P i+1 respectively. The error vector -→ e i is defined as the vector connecting the actual P i point to its desired origin:

-→ e i = ----→ P i -P i (5.8)

The normal vector to the desired path at point P i (termed as -→ d i ) is also computed by using the two adjacent desired points:

-→ d i = -----→ P i+1 -P i | -----→ P i+1 -P i | (5.9)
Using the plane normal vector -→ n , the orthogonal component of the error to the desired path in point P i (E i ) can be achieved as:

E i = -→ e i .( -→ n × -→ d i ) (5.10)
At this point different choices are available. Choices such as maximum error value (max{e i }) or mean error value (avg{e i }) are among interesting functions to define the quality of machining workpiece for a given set of error values along its machining path [START_REF] Xie | Posture Optimization Based on Both Joint Parameter Error and Stiffness for Robotic Milling[END_REF][START_REF] Caro | Workpiece placement optimization for machining operations with a KUKA KR270-2 robot[END_REF].

The errors observed in any machining process can mathematically be distinguished into two main categories: First is the average error value which is, in our application, a circle with a radius (r p ) differing from the desired circle radius (r) by the average error value: 5.11) Where: e = Σe i N (5.12)

r p = r + e (
The remaining part of the error can be viewed as a scattering around the average resulted path. In one hand, one can conclude from the literature that the effect of joint parameter errors has a lower impact on the tool displacement error compared to error sources such as compliance effect [START_REF] Xie | Posture Optimization Based on Both Joint Parameter Error and Stiffness for Robotic Milling[END_REF] and on the other hand, the effect of cyclic errors can more contribute to the scatterings along a defined path as shown in the experiments presented in [START_REF] Vakilinejad | Identification and Compensation of periodic gear transmission errors in Robot Manipulators[END_REF] In search of a suitable cost function imposing stronger constrain on cyclic errors as well as compliant behaviour error, we propose the variance of the error values (var{e i }) which thanks to its definition, eliminates the effect of average error value:

f c (x c , y c , z c , θ) = Var e x c , y c , z c , θ =
Σ(e ie) 2 N (5.13) Figure 5.5 depicts two scenarios of high and low values of the error variance.

In this work we have focused our effort on reducing the compliance behaviour and cyclic transmission system error of the joints. However the proposed cost function is defined generically and capable of implementation for errors rising from any other source. Modifications can be applied specially in the final step to better describe any arbitrary problem statement.

Path Definition

Undoubtedly finding seemingly optimal solutions will only result in an efficient practical outcome when the defined model in the research state is a strong representative of the actual problem. In our problem, the path definition addresses the systematic procedure we aim to take to mathematically describe the path presenting the geometrical lines and curves that the robot end-effector must follow to accomplish the trimming process.

To be able to completely distinguish a robot pose from another, one should define six parameters corresponding to the general six degree of freedom in Cartesian space. Since the transformation from the Cartesian space to joint variable space does not provide a explicit expression for a mathematically explicitly defined path, the remedy to provide the joint variable trajectories is to discretize the path with a reasonable distance gap in between adjacent points. One can easily conclude that the decision over the number of points describing the path (I.e. the distance gap defined as ∆d in Figure 5.6) is vividly a decisive decision. Describing the path by gazillion number of points will lead to a drastic increase of computation time which for an optimization process means less number of costfunction evaluation for potential solutions in a given time and thus lowers the chance of finding a better solution. On the other hand, describing robot duty with low number of points would clearly cause in inaccurate results and a poorly defined joint variable trajectory.

∆d ∆d P i-1 P i P i+1
To find a trade off in the number of points to describe the path of our machining problem we propose a simple procedure. Firstly, as a accurately defined path with distance gap of ∆d = 0.1(mm) cost function is evaluated for 50 different points widely spread in the work-space of the robot. These values are then stored and used as reference values for these locations in the workspace.

The distance gap (∆d)is then gradually increased by steps of 0.1(mm) to final value of 50(mm). At each step, cost function is evaluated and compared to the previously defined reference values. Maximum changes among different 50 locations in robot workspace is attributed as the discretization error of the corresponding distance gap. Results are exhibited in Figure 5.7.

Transmission System Error

Elastostatic Error As shown in figure 5.7, the discretization error increases by increasing the distance gap (∆d). Even though the increase of error was expected as a result of increasing the distance gap, the main conclusion drawn from the figure is the efficient value of distance gap. As shown in the figure, even by increasing the distance gap to 50(mm), maximum discretization error evolves to a value of 2%. Hereafter we pick a distance gap of 10(mm) for a compromise between the computation load and indispensable situations such as vicinity to joint limits and singular positions.

It is worth mentioning that the discretization constant (distance gap) described above is only used in the the optimization level. In the final stage, the robot trajectory developed by the machining software uses smaller gaps and fills the successive points by a pre-defined method (linear or constant evolution).

Constraint Implementation

Joint Limits

Among all the constraints and limiting conditions in robotics, respecting joint limits is of most importance. joint limits are imposed due to physical limits of robot structure and violating them may cause irreversible damage to the comparatively costly apparatus such as industrial robot arm. Although almost all of the commercially available robots come with external controllers taking extra care of respecting this limiting constrain, implementing this important constrain in the research level will ensure us of a practically feasible solutions offered by optimization procedure.

Since joint limit constrain is a binary condition one by nature, defining it as a binary condition is a conventional yet popular approach to this criteria as proposed in [START_REF] Caro | Workpiece placement optimization for machining operations with a KUKA KR270-2 robot[END_REF][START_REF] Xiong | Stiffness-based pose optimization of an industrial robot for five-axis milling[END_REF][START_REF] Ghariblu | Path Planning of Complex Pipe Joints Welding with Redundant Robotic Systems[END_REF][START_REF] Doan | Optimal robot placement with consideration of redundancy problem for wrist-partitioned 6R articulated robots[END_REF]. However, intelligent changes in defining joint limit constrains may not only appear to be more efficient in ensuring of respecting this constraint, but also increases the chance to better control robot behaviour in the vicinity of the joint limits as well as optimization procedures . Among inspiring mathematical presentation of joint limits, function proposed by [START_REF] Tsai | Workspace geometric characterization and manipulability of industrial robots[END_REF][START_REF] Marchand | Using the task function approach to avoid robot joint limits and kinematic singularities in visual servoing[END_REF][START_REF] Zghal | Efficient gradient projection optimization for manipulators with multiple degrees of redundancy[END_REF][START_REF] Liegeois | Automatic supervisory control of the configuration and behaviour of multibody mechanisms[END_REF][START_REF] Xiao | Redundancy and optimization of a 6R robot for five-axis milling applications: singularity, joint limits and collision[END_REF][START_REF] Ito | Kinematic Control of Redundant Manipulators for Admitting Joint Range of Motion Maximally[END_REF] are notable. In this work, we have adopted Tsai function to impose the joint limit constrain on the optimization problem:

P = 1 + exp -k × n ∏ j=1 (θ j -θ min,j ) × (θ max,j -θ j ) (θ max,j -θ min,j ) 2 (5.14)
According to the nature of pattern search optimization method performance, Tsai function offers two main advantages: The drastic increase of the function in close inner vicinity of the joint limit strongly penalizes the cost function while the considerable majority of solutions falling within the joint range stay fairly unchanged. Secondly, the steep curve created by Tsai function in neighbourhood and beyond joint limits enable pattern search to converge to the acceptable region even when the initial location of a population member falls out of the robot motion range.

Singularity Avoidance

Next important constrain in any robotic application is singularity avoidance. In robotics, trajectories are specified in either joint or operational space. In machining applications, operations (involving machining features, feed rates, positioning etc.) are generally defined and programmed in Cartesian space and requires further transformations to translate the desired point-time sequences from the operation space to joint coordinate to pass the proper commands to actuator control unit [START_REF] Siciliano | Robotics: modelling, planning and control[END_REF] Cartesian velocity requested by the operation is related to the joint angular velocity by the linear Jacobian transformation (J(Θ)) is used:

[ Θ] = J(Θ) -1 [ Ẋ]
(5.15)

Where [Θ]and[ Θ] are the 6 × 1 array of angular values and velocities respectively and [ Ẋ] is the time derivative of Cartesian components (position and orientation). Expression 5.15, vividly reveals that in certain combination of [Θ], where the determinant of matrix J(Θ) is null, the required angular velocity (of at least one joint) will rise beyond the capacity of joint actuator. Impeded control algorithm, increased force or torque imposed on actuators reduced mobility are also among inconveniences which may arise as consequences of singularities Donelan, 2010. In Xiao and[START_REF] Xiao | Redundancy and optimization of a 6R robot for five-axis milling applications: singularity, joint limits and collision[END_REF] Huan, have experimentally reported the maleffects of close singular configurations in surface quality of workpiece in a robot based milling operation.

Extensive research has been conducted to the topic of singularity in robot arm manipulators [START_REF] Wang | A study of the singular configurations of serial manipulators[END_REF][START_REF] Lipkin | Enumeration of singular configurations for robotic manipulators[END_REF][START_REF] Cheng | Study and resolution of singularities for a 6-DOF PUMA manipulator[END_REF]Wang et al., 2018a) and how to either avoid or tackle them in robotic applications [START_REF] Garcia | Maximal Singularity-Free Orientation Subregions Associated with Initial Parallel Manipulator Configuration[END_REF][START_REF] Krastev | Passing Through Jacobian Singularities in Motion Path Control of Redundant Robot Arms[END_REF]Wang et al., 2018a). Proposed methods to tackle singularities are mainly applicable in path planning by slightly deviating from singularityincluded path and require rather complex mathematical work and programming in actuator level [START_REF] Xu | Practical approaches to handle the singularities of a wrist-partitioned space manipulator[END_REF]. Aforementioned drawbacks clearly explain the necessity of avoiding singular configurations in optimization level.

Thanks to the wrist-partitioned design of Stäubli TX200 robotic arm, a singular configuration of robot is only met in the following three conditions [START_REF] Doan | Optimal robot placement with consideration of redundancy problem for wrist-partitioned 6R articulated robots[END_REF][START_REF] Hayes | Singular configurations of wrist-partitioned 6R serial robots: A geometric perspective for users[END_REF]):

• when rotation vector of joint 6 and 4 align. This singularity condition occurs when θ 5 = 0 and is known as Wrist singularity.

• In situation where the intersection point of wrist axes (4,5 and 6), lies on the spread plane of rotation vector of axes 2 and 3, robot encounters another type of singularity configurations termed as Elbow Singularity. For Stäubli TX200, θ 3 = 0 satisfies this condition.

• In robot postures where intersection point of wrist axes lies along the rotation vector of axis one, Shoulder Singularity takes place. Unlike the previous singularity conditions, configurations in which a wrist-partitioned robot undergoes shoulder singularity can not be specified by only one joint parameter.

To avoid the singular configurations in this step, likewise the joint limit condition we propose a penalization approach. In this approach, calculation of the Jacobian is avoided and thus a huge load of computation is relieved. A singularity coefficient (S) is proposed to penalize the results from the cost-function and is composed of three sub-coefficients S w , S e and S sh which represent the singularity closeness with respect to wrist, elbow and shoulder singular configurations respectively: S = S w × S e × S sh (5.16) Expressions proposed of the sub-coefficients S e and S w are straight forward and chosen in such way that their value will stay practically very close to 1 and dramatically rise to 10 only in vicinity of the singular conditions. Expressions for the mentioned coefficients are as following:

S w (x c , y c , z c , θ) = 1 + exp -k × min({|θ 5 i |}) -θ *
(5.17)

S e (x c , y c , z c , θ) = 1 + exp -k × min({|θ 3 i |}) -θ * (5.18)
Where the operator min indicates that the singularity criteria would be quantified on the most potential critical point of the path and value of δ * can be freely adjusted.

Proposing an expression representing the shoulder singularity criteria contains however further investigations. As presented in Figure 5.8 the distance between the intersection point of wrist axes (C) from the rotation vector of axis 1 (Or in other words the z axes of robot base frame), can be calculated using the following equation:

(θ 2 , θ 3 ) = l 1 + l 2 sin(θ 2 ) + l 3 sin(θ 3 + θ 2 ) (5.19)
It is evident that the distance r is only a function of θ 2 and θ 3 (r(θ 2 , θ 3 )). According to the description of shoulder singularity, the operator must avoid configurations where r is zero. In search for the sets of (θ 2 , θ 3 ) corresponding to the shoulder singularity configuration, the value of function r(θ 2 , θ 3 ) is computed for a practical working range of second and third joint (-120 ≤ θ 2 , θ 3 ≤ 120). Results are presented in the contour shown in Figure 5.9.

As shown in Figure, two zones corresponding to a null r and thus shoulder singularity, termed as z 1 and z 2 , are present. z 2 is related to the sets of configurations in which point C lies on the negative side of axes z and thus automatically removed from the feasible solutions by imposing the lower boundary of z component to zero for the optimization process. On the other hand, zone z 1 involves feasible solutions and thus requires further treatments. the shape created by zone z 1 in the r(q 2 , q 3 ) in Figure 5.9 is close to a straight line. The sets of values (θ 2 , θ 3 ) forming z 1 were extracted and passed through least square method to find the equation of the closest line (L) (values defined in [Deg]): L(θ 2 , θ 3 ) : θ 2 + 0.37 × θ 3 + 13.5 = 0 (5.20) The maximum error in modeling zone z 1 with a straight line (maximum distance of the points forming z 1 from line L) does not exceed 7.1 [Deg]. In light of above mentioned investments on expressing the shoulder singularity, the proper expression for penalizing the cost function with respect to this singularity type is based on the minimum distance of (θ 2 i , θ 3 i ) (for i = 1 : N where N is number of the points used for path definition). For a proposed placement of the workpiece using expression:

d i = |A × θ 2 i + B × θ 3 i + C| √ A 2 + B 2 (5.21) θ Z 1 Z 2 r > 0.2(m) r > 0.2(m) L θ FIGURE 5
.9: Contour of distance between wrist intersection point and vertical axis of first joint of the robot (r(θ 2 , θ 3 ))

In which A, B and C are constants of line L and hold the values of 1, 0.37 and 13. 5[deg] respectively. Finally, the Shoulder Singularity coefficient can be formulized as:

S sh (x c , y c , z c , θ) = 1 + expk × min({|d i |}) -7.1 (5.22) To avoid singular conditions in robotic applications, other approaches such as condition number presented in [START_REF] Khan | The kinetostatic optimization of robotic manipulators: the inverse and the direct problems[END_REF][START_REF] Khan | The kinetostatic optimization of robotic manipulators: the inverse and the direct problems[END_REF]tuned weights in Huo and[START_REF] Huo | The self-adaptation of weights for joint-limits and singularity avoidances of functionally redundant robotic-task[END_REF] are proposed in the literature. However these approaches are Jacobian based approaches. Computing jacobian matrix and, for some methods, its eigenvalues compared the approach presented in this work is more time consuming.

To compare the impact of presented singularity avoidance criteria implementation to the conventional methods, the predefined cost function was evaluated in 1000 points in the workspace along with their corresponding singularity coefficient using both strategies. Condition number presented in [START_REF] Khan | The kinetostatic optimization of robotic manipulators: the inverse and the direct problems[END_REF] was chosen as a jacobian based approach. The computation time for both trials are brought in Figure 5.10.

By having the values of penalizing coefficients related to the kinematic limits (P and s) of the robot, the final format of the cost function turns to be:

F c = P × S × f c
(5.23)

Optimization Method

In engineering, as well as many other fields, researchers have always faced problems in which the impact of various variables is inevitable. The daily progress and advancement of technology and stronger yet more complex modeling of physical phenomenon continuously increase dimensions of highly nonlinear problems where researchers can no more relay on conventional algebraic or linear methods to come up with best combinations, arrangements or designs in their problems.

Exploiting optimization methods to expose efficient solutions is common in different engineering applications. From optimizing the solutions in water resource planning Nicklow et al., 2009,design A brief glance on the workcell design optimization problem in the literature reveals that a hand full of optimization methods have been utilized for this goal. Among all, Genetic Algorithm (GA) [START_REF] Caro | Workpiece placement optimization for machining operations with a KUKA KR270-2 robot[END_REF] and particle swarm [START_REF] Doan | Optimal robot placement with consideration of redundancy problem for wrist-partitioned 6R articulated robots[END_REF] methods are the most popular ones due to the nonlinear nature of the problem. Their independency of the gradient of the cost function and their capability in dealing with highly constrained objectives are among notable advantages of these optimization methods.

However, in this work, we propose to approach the optimization problem with a modified free pattern search algorithm to discover its efficiency of this method in robotic workcell design problem.

The conventional pattern search introduced by Hooke and Jeeves, 1961 wields a simple yet strong logic for choosing the best direction for marching from an either randomly or arbitrary initial guess inside the design variables space.

In this optimization method, a set of N points (equal to the desired population size) is chosen either randomly or uniformly distributed in the n-dimensional space of variables. Then each point of the initial population successively undergo the following process: each point is initially set as the current point (X). A geometrical structure termed Grid [START_REF] Wen | Free Pattern Search for global optimization[END_REF] is formed around the current point. The group of points forming the grid is hereafter shown by ψ. Different strategies have been introduced for forming the grid.

In the conventional method, as well as later modified methods such as the one proposed in [START_REF] Wen | Free Pattern Search for global optimization[END_REF], the grid contains as many as two times of the number of variables (2n) of points. These points each vary only in the value of one variable in positive and negative direction with the distance of δ (the marching parameter) from the current point. Other strategies have been also introduced [START_REF] Lewis | Rank Ordering and Positive Bases in Pattern Search Algorithms[END_REF][START_REF] Dolan | Pattern search behavior in nonlinear optimization[END_REF].

The principal modification applied on the pattern search algorithm in this work is grid geometry. As mentioned earlier, in conventional grid constructing approach, each point of the grid differs only value of one component of its coordinate in the variable space from current point. However, the gradient of the cost function in a certain point may decrease in a direction other than orthogonal directions. Moreover, The gradient of the cost function may decrease in more than one direction. In this case, it is more probable to have a steeper degradation in a direction closer to the resultant of decreasing orthogonal directions.

Forgoing reasons convinced us to use a full grid which contains all the possible direction of changes in variables. As each variable has three different possibility of change (In increasing, constant and decreasing direction), grid consists of 3 n number of points. However, in a situation where all variables are constant, the corresponding point of the grid is equivalent to the current point (X) and thus eliminated from the grid. Thus the number of points in this case equals to 3 n -1 where n is number of variables. Full-sized grid format for a two variable space is depicted in Figure 5.11 where dashed-lined points are the additional points in this grid structure compared to the conventional grid structure. As explained in [START_REF] Wen | Free Pattern Search for global optimization[END_REF], a matrix containing the values of changes in variables termed as M is utilized to form the grid by adding each of its columns to the coordinates of the current point. In full-sized grid structure, matrix M contains equal 3 n -1 columns. The corresponding M matrix in this case is as follows:

M =      -δ 1 0 δ 1 • • • -δ 1 0 δ 1 -δ 2 -δ 2 -δ 2 • • • δ 2 δ 2 δ 2 . . . . . . . . . . . . . . . . . . . . . -δ n -δ n -δ n • • • δ n δ n δ n      3 n -1
(5.24)

All the points in the grid are then evaluated by the cost function and in case the best point of the grid offers a better solution than the current point, two points will be nominated to replace the current point. First point (X n1 ) is the best point of the grid (φ) and second point (X n2 ) is the reflection of the current point across the grid best point:

X n1 = φ X n2 = 2 × φ -X (5.25)
Trivially, the evaluation by cost function will decide whether to replace the current point for coming iteration. If the points in the grid do not offer a better alternative to the current point, the searching increment (δ) which has been set initially to ∆ 0 shrinks to a smaller value by the constant κ:

δ = κ × δ (5.26)
The algorithm continues to an instant when the value of δ is smaller than a predefined parameter ∆ min . At this point, the algorithm continues with the next point in the initial population. A detailed chart of the algorithm is depicted in figure 5.12.

Determination of Searching Parameter

One of the main advantages of using free pattern search algorithm compared to the heuristic algorithms such as GA, is the authority given to the operator to choose different properties related to search marching parameter (∆) throughout the optimization process. One should note that a large value of ∆ will speed up the optimization process but on the other hand it may miss comparatively smaller zones containing local minimums.

To determine a proper value of ∆ 0 , we require a basic knowledge on the effect of changes of each variable on the variation of cost function. In problems such the one focused in this work, a solution offering the best evaluation of the cost function may not necessarily be an efficient solution in practice. Practically-efficient solutions should also possess a certain level of robustness. Robust solutions in our application are the ones in which small variations in the design variables will not cause a drastic changes in the cost function evaluation.

The placement of the workpiece in the workspace of the robot is prone to the inaccuracies rising from either human or machine operators executing the task. Although in precise manufacturing processes, the transformation matrix between robot base frame and workpiece is updated using adequate sensors and probes before the task, however, this modification process does not offer any advantage in reducing the joint variable errors due to the misplacement of the workpiece.

To overcome the robustness problem, in this work we suggest to wisely determine the termination parameter ∆ min . By assigning this parameter, at minimum, to approximately the accuracy of the placement of work-piece, we tend to reduce the probability of the optimization method to stop at a local minimum trapped in a zone considerably smaller in size compared to the placement accuracy.

As mentioned earlier, in this work we try to optimize three positioning variables corresponding to the three coordinates of the center of work-piece (x c , y c , z c ) and one rotational parameter standing for the deviation of the work-piece surface in z direction of the robot base frame (θ). To sake simplicity, we assume that all three positioning variables behave in the same manner and thus only one set of values are assigned to the ∆ 0 and ∆ min . The fourth variable (θ) on the other hand, requires different treatments as explained later.

To observe the changes in the cost-function with respect to the positioning variables, function is evaluated throughout a plane at x = 700(mm) in the workspace of the robot using a grid with equally spaced points with distance of 1(cm) each for ranges of -1000(mm) to 1000(mm) and 0 to 1800(mm) in y and z directions respectively. Figure 5.13 presents the results.

Out of robot reachability

Out of robot reachability d = 100(mm) To determine the searching step parameter of the last variable (θ), same procedure is utilized. The effect of changes of frame rotation on the value of cost function is simulated for -90[deg] < θ < 90 [deg]. Center of the workpiece for observation of the frame rotation effect is located in [1000, 600, 500](mm). Figure 5.14 presents the evolution of the cost function with respect to the changes in frame rotation parameter. As shown in Figure 5.14, the sensitivity of the cost-function with respect to frame rotation (θ) is lower compared to the coordinating components. Therefore by setting ∆ min to 10[deg], the evolution of the cost function can be tracked. Figure 5.14 also shows that the effect of frame rotation is considerable in the robot related machining errors. For the mentioned location of the workpiece center, the variance of TCP positioning error can be reduced more than 50% by only rotating the workpiece frame from its initial posture (θ = 0) in the appropriate direction.

Out of Robot Reachability

Simulation Results

Previously explained optimization method was used to minimize the cost function defined earlier. The upper and lower limits used for the four design variables in optimization process was chosen to cover almost every reachable point of the front zone of the robot:

       0(mm) ≤ x c ≤ 2000(mm) -2000(mm) ≤ y c ≤ 2000(mm) 0(mm) ≤ z c ≤ 2000(mm) -90(deg) ≤ θ ≤ 90(deg) (5.27)
The optimization method was applied on the problem using both proposed and conventional grid structures to serve as a means of comparison. The termination criteria were: a) maximum cost function calls of 2600 times) A solution offering a value of cost function for less than 0.01(mm). Figure 5.15 illustrates the cost function value of the best found solution by the optimization method as function of function calls. log scale is adopted for y axis due to massive changes of cost function values. Another conclusion drawn form Figure 5.15 is regarding the performance of the proposed grid structure. One should note that for a 4-variable optimization problem, function calls of the proposed grid is 10 times more than the one of the conventional grid structure (80 times and 8 times respectively). Nevertheless, Figure 5.15 depicts that this computational cost of the full-sized grid has been compensated by its efficiency in converging to a proper solution. This demonstrates that the free pattern search method using the full-sized grid can converge in less number of grid constructions. However it is worth mentioning that a grid structure which exponentially increases the number of function calls with respect to the number of function variables for a single iteration, may only be efficient for problems with low number of variables.

Figure 5.15 also shows the benefit of the presented way of constraint implementation. The initial points for both optimization executions were out of robot reachability (cost function values 1) but the optimization process continued to reach a desirable solution.

To compare the compliance and transmission cyclic errors observed in the workpiece trimming trajectory, in the optimal placement of the work-piece and an intuitive placement, a point in the workspace of the robot is chosen to represent a probable intuitive placement of the workpiece. This point along with the optimized solutions are presented in table 5.1. An important conclusion deduced from Figure 5.16 is that although the the dimension of the cost function output in (mm), the practical feasibility of decision making over low order numbers as the variance of the error should not be judged by total error observed in the TCP displacement. It is vividly observed in Figure 5.16 that a workpiece placement having a 1.7(mm) of variance in normal errors along the path may include up to 4.9(mm) of deviation along the path which is much higher than the robot repeatability (about 0.05(mm)).

We should also mention that although the definition of cost-function is based on removing the constant displacement error (avg {e i }) from the error values, but due to the scattering trend of these errors, the defined cost-function was fully capable of proportionally explain the total error values. However, this cost-function may not be as efficient in the cases in which robot TCP undergoes a constant load (Such as gravitational forces).

Conclusion

Due to the inconveniences of the off-line compensation processes for improving machining accuracy by reducing robot-related error sources, in this chapter, we proposed an optimization procedure to find a close-optimal workcell design. Since the off-line compensations such as the one presented in chapter 3 require post-processes of the error models for a given workpiece location and usually the improvements are prone to inaccuracy of the error modelings, optimizing the workpiece location is of interest.

The workcell design in this chapter was approached by replacing the workpiece with respect to the robot base frame and thus dismantling of the robot structure from its initial placement is avoided. The targeted machining process was CFRP trimming Chapter 5. Optimal Workpiece Placement In Robotic Cells task and the workpiece was presented by a circular shape.

The goal was to reduce the errors originated from both elastostatic behaviour and periodic joint transmission errors. The proposed cost function was designed using the variance of the normal on the feed rate direction simulated from the previously generated error models. The kinematic constraints of the robot system namely joint limits and singularity avoidance were implemented using a penalization method. A novel way was proposed for implementing the singularity avoidance applicable on any wrist-partitioned robot structure which massively reduced the computation time of the optimization process. This method can be adopted specially for online path planning problems.

The optimization method adopted for this problem was a modification of the free pattern search algorithm in which a full-sized grid was employed to facilitate the decision making over the progress direction. The performance and efficiency of the proposed optimization procedure are validated by simulation results.

Chapter 6

Conclusion and Future Work

The concept of the manufacturing industry has to inevitably adapt to the fast-growing and demanding features of the recent signs of progress in other fields of engineering. The rapidly evolving complexity and dynamic nature of advancements in technology can only continue to thrive by leaning on a flexible production concept ready to reform and comply with short noticed changes at low cost. Industrial Robot (IR) arms have shown appealing advantages to lead this revolutionary breakthrough.

Although IRs offer vast merits such as extended workspaces and reconfigurability, however, they come with shortcoming required for further treatments and compensations to meet the strict standard accuracy demanded in high-tech sectors.

The contribution of each non geometrical error source was examined in chapter 2. The examinations were executed in loaded and unloaded conditions on a case study robot (Stäubli TX200). results showed that errors resulted from the compliance behaviour of robot structure are dominant. This error source was divided into two categories namely self-gravity effect and external load effect. The former category imposed higher values of errors compared to the latter.

Along with compliance error source, inaccuracies originated from periodic errors of joint transmission systems and robot thermal behaviour were evaluated. According to the experimental results, the contribution of each error source to the non geometrical errors for a quasi-static loaded condition was: 51%, 36%, 9% and 4% for self-gravity effect, external load effect, nonlinear joint transmission errors and thermal behaviour respectively for the test conditions reported in this chapter.

Chapter 3 was dedicated to identification and compensation of two min contributors to non geometrical errors namely compliance behaviour and joint transmission errors. To overcome the compliance behaviour, an elastostatic model was developed capable of accounting the torques resulted from three main torque imposition sources: gravity (termed as self gravity effect in the thesis), external loads and selfbalancing system attached to the robot. The stiffness values of robot joints were determined experimentally. The performance of the generated elastostatic model was validated experimentally. Although experiments revealed an acceptable efficiency of the model, however, the joint stiffness values obtained from the robot supplier showed better performance.

To compensate for the effect of nonlinear joint transmission errors, a mathematical model of this error type was required. A novel method was proposed for the identification of this error source based on constructing an optimal Cartesian path. The path was designed to enable the operator to identify the error model for all Chapter 6. Conclusion and Future Work targeted joints simultaneously. Therefore the main advantage of this method was its considerably shorter execution time compared to the available methods in the literature. An offline compensation procedure was presented with which the path accuracy for a straight line increased by 40% according to the obtained experimental data.

Chapter 4 deals with a practical case of a robotic based machining process. In this chapter, a machining setup composed of a robotic arm, spindle and machining tool designed for cutting Aramid honeycomb material with the aid of ultrasonic technology was adopted as a case study. The main concern in this machining process was the presence of a high level of undesired geometrical errors observed in resulted workpieces. Geometrical machining errors were measured using a practical approach presented in the chapter.

The problem was initially approached by developing a force model based on experimental force values for the chamfering process of honeycomb materials. The elastostatic behaviour of the robot structure was simulated using empirically driven tool-workpiece interaction forces. Mismatching simulation results with the actual error values inspired us to include the compliance behaviour of the triangular knife (machining tool) in the total robot-tool compliance model.

To integrate the tool compliance into the process, tool stiffness was modeled firstly using the structural stiffness matrix (SSM) method. The model was verified by experimental measurements. The results showed that the compliance effect of the knife structure mainly comes from the vertical direction. The entire robot-tool system was then modeled by a 7 axis robot assuming the vertical direction of the knife a virtual joint. The behaviour of this new system was analyzed and results showed that the effect of the knife compliance was more than 8 times greater than the one of the robot structure.

To decrease the machining errors, different loading conditions applied to the machining tool was examined. Simulations showed that a non-uniformly distributed load along the longitudinal axis of the knife can best describe the angular chamfering errors measured during the experimental sessions. Using the compliance model of the knife structure and machining forces model in this process, an offline compensation method was proposed based on an iterative algorithm. In this compensation method, the algorithm was designed to output a secondary chamfering angle corresponding to a certain desired chamfering angle. The experience proved the efficiency of the compensation method which resulted in a reduction between 40% to 90% of the geometrical errors. However, the chamfering machining process was observed to possess a non-repetitive nature in terms of machining forces.

To overcome the inconveniences encountered in post-process stages in off-line compensation approaches, in chapter 5 an optimization method was proposed to minimize the non geometrical errors of an IR by replacing the workpiece in a closeoptimal configuration. The targeted machining task was trimming of a circular Carbon Fiber Reinforced Plastic (CFRP) workpiece. The subjected error sources were robot compliance behaviour and joint transmission system errors.

A cost function was proposed based on the variance of the components of displacement errors normal to the feed rate direction. Among the contributions of this chapter, the constraint implementation approach is notable. Joint limits and singularity avoidance were implemented in the function using penalizing strategy. A straight line was fitted to the sets of values of joints 2 and 3 for which the robot faces shoulder singularity configurations. This line was implemented in the singularity avoidance criteria to avoid jacobian computations. The method of constraint implementation drastically increased the computational efficiency of the optimization process.

A modified pattern search algorithm was employed to execute the optimization process. A full-sized grid was adopted in the algorithm to facilitate the decision over the progression direction for a proposed initial point. Simulation results showed that a great portion of non geometrical errors originated from the robotic arm can be eliminated by replacing the workpiece in a close optimal location. By obtaining the above-mentioned conclusions from throughout this thesis, some proposals for the continuation of research in this area worth mentioning. To increase the accuracy of the elastostatic model of the robot, adopting a stiffness identification method in which the deformation of robot links are eliminated in the measurement process (such as dynamic stiffness identification procedure) can serve better.

As concluded in chapter 4 that the ultrasonic knife structure causes the dominant effect on the geometrical errors observed in the honeycomb workpieces, an efficient approach is to enhance the compliance behavior of the tool. Later can be achieved by focusing on redesigning the knife geometry considering the amplitude of the dynamic stress along the knife body. High dynamic stress would result in weak fatigue performance of the tool.

The new approach proposed in chapter 5 for implementing the singular avoidance criteria in the optimization problem consists of three coefficients representing three different singularity conditions in a wrist portioned robotic arm. The feasibility of using the presented way to develop a novel method for describing the manipulability of a wrist oriented robot can be an interesting research subject.

Le chapitre 2 examine les différentes sources potentielles d'imprécision dans les processus d'usinage impliquant des robots. Une analyse documentaire du sujet est préparée. Les principaux facteurs d'erreurs de positionnement du robot sont mentionnés. La part de contribution pour trois sources d'erreurs non géométriques, à savoir le comportement de conformité, l'imprécision induite thermiquement et les erreurs de transmission non linéaires des boîtes de vitesses robotisées, est évaluée pour notre étude de cas. L'évaluation est réalisée sur la base de données expérimentales pour le bras robotique Stäubli TX200 chargé et déchargé.

En raison de l'importance des sources d'erreurs observées expérimentalement, deux sources d'erreurs non géométriques, le comportement de conformité de la structure du robot et les erreurs non linéaires des systèmes de transmission du robot ont été choisies pour des recherches plus approfondies au chapitre 3. Le comportement de conformité a d'abord été ciblé. Pour compenser les erreurs résultant de cette source d'erreur, un modèle élastostatique du robot faisant l'objet de l'étude de cas est développé, impliquant différents aspects dont l'effet d'auto-gravité de la structure du robot, l'impact des charges externes et le mécanisme d'auto-équilibrage attaché au robot. En tant qu'élément vital du modèle élastostatique, les valeurs de rigidité des articulations du robot sont évaluées de manière expérimentale. Le modèle est ensuite utilisé pour compenser la trajectoire du robot pour une condition de charge donnée. Des valeurs expérimentales sont présentées pour montrer l'efficacité de la performance du modèle élastostatique.

Le reste du chapitre traite des efforts consacrés à l'identification et à la réduction des erreurs de transmission non linéaires des articulations du robot. Dans cette partie, les erreurs périodiques d'un robot d'étude de cas en laboratoire sont identifiées à l'aide d'un nouveau chemin proposé dans l'espace cartésien, ce qui réduit considérablement le temps du processus d'identification des erreurs par rapport aux méthodes proposées dans la littérature. L'erreur identifiée est utilisée pour une procédure de compensation hors ligne. La performance et l'efficacité des procédures d'identification et de compensation sont validées sur la base des résultats expérimentaux.

L'une des principales applications du Creneau industriel étant la découpe par ultrasons de matériaux en nid d'abeille aramide, le chapitre 4 donne un aperçu des caractéristiques spécifiques de ce processus d'usinage exécuté par le bras robotique Stäubli TX200. Afin d'examiner en profondeur le comportement du robot dans les conditions d'usinage, les forces d'interaction de l'usinage et les valeurs d'erreur géométrique ont été quantifiées lors de plusieurs sessions de tests. Un modèle de force est généré pour estimer les forces d'usinage pour des caractéristiques données dans le processus de chanfreinage en nid d'abeille. Ce modèle est ensuite utilisé pour simuler les déviations et le déplacement du point central d'outil du robot (TCP) dans un processus de chanfreinage d'un matériau alvéolaire.

En raison de la valeur élevée du comportement de conformité de l'outil d'usinage dans ce processus d'usinage, une nouvelle méthode est proposée pour contrer l'effet de conformité pour un système conforme à un robot-outil. Le modèle de conformité de la lame triangulaire à ultrasons, qui sert d'outil d'usinage dans cette application, est développé selon la méthode SSM (Structural stiffness matrix). Le modèle est examiné sur place avec des valeurs expérimentales. La direction la plus conforme de l'outil est modélisée comme un joint virtuel de la structure du robot. L'effet d'une charge externe sur ce nouveau robot 7 axes est simulé à l'aide de la rigidité torsionnelle équivalente de l'outil et des valeurs de rigidité des articulations du robot obtenues expérimentalement.

En raison de la grande influence de la structure de l'outil d'usinage, la condition de charge appliquée à l'outil est modélisée en utilisant une charge non uniformément répartie sur toute sa longueur. Une procédure hors ligne pour diminuer les erreurs géométriques observées dans les pièces en nid d'abeilles est proposée sur la base d'un algorithme itératif. La procédure utilise la force d'usinage générée et le modèle de conformité de l'outil. Les valeurs d'erreur compensées et non compensées sont indiquées.

Le chapitre 5 présente une stratégie pour optimiser le placement d'une pièce en fibre de carbone soumise à un processus de coupe. Les critères de ce problème d'optimisation sont de réduire l'effet du comportement de conformité et des erreurs de transmission non linéaires dans la pièce finale. Une nouvelle fonction de coût basée sur la variance des erreurs induites est proposée. Pour éviter des configurations singulières du robot, une nouvelle approche est introduite basée sur la pénalisation de la fonction de coût en fonction de la distance à laquelle se trouve une configuration proposée, à partir de trois conditions singulières différentes qu'un robot divisé au poignet peut rencontrer.

La méthode d'optimisation utilisée dans ce chapitre est un algorithme de recherche de modèle libre modifié avec une construction de grille généralisée. Cette grille aide l'algorithme d'optimisation à prendre de meilleures décisions sur le choix de la meilleure direction à prendre pour aller de l'avant. Les résultats de simulation sont présentés pour montrer la performance de la fonction de coût proposée et de la méthode d'optimisation. La figure A.2 présente les grandes lignes de cette thèse.

A.2 Chapitre 2

Les inexactitudes des robots peuvent provenir de différentes sources. Selon quelques travaux de recherche tels que [START_REF] Kim | Robotic Machining: A Review of Recent Progress[END_REF][START_REF] Schneider | Experimental investigation of sources of error in robot machining[END_REF][START_REF] Gong | Nongeometric error identification and compensation for robotic system by inverse calibration[END_REF] ces sources peuvent être classées en deux catégories différentes qui sont à savoir Geometrical et Non geometrical erreurs sources :

• Erreurs géométriques

Ces erreurs proviennent des tolérances des processus de fabrication, d'assemblage et d'installation des éléments du robot. Ils sont présents à la fois dans les maillons et dans les articulations, donc ils existent toujours indépendamment de l'état de fonctionnement du robot. Plutôt que les tolérances de fabrication qui entraînent une déviation des pièces robotisées fabriquées par rapport à leurs valeurs nominales, les désalignements des articulations et des axes sont également classés dans cette zone.

• Erreurs non géométriques D'autres types d'erreurs qui ne peuvent pas être compensées par une modification permanente du modèle géométrique du robot sont appelées erreurs non géométriques. La caractéristique la plus importante de ces erreurs est leur dépendance à d'autres facteurs tels que la charge appliquée et la régularité de la trajectoire en provoquant le comportement de conformité du robot, des facteurs environnementaux tels que les variations de température dans le temps provoquant la déformation des liaisons du robot ou des erreurs liées à la configuration comme les irrégularités dans les variables communes.

La réduction des erreurs géométriques du robot est généralement basée sur la mise à jour du modèle mathématique servant de cinématique avant en utilisant des décalages obtenus expérimentalement au lieu des valeurs nominales du manipulateur du robot Nubiola and Bonev, 2013. Les méthodes d'étalonnage géométrique peuvent différer de diverses manières. Le modèle mathématique adopté pour la transformation des trames de base à la fin (modèle géométrique direct) peut varier de la mise à jour des paramètres Denavit-Hartenberg conventionnels (également appelés paramètres DH) à la transformation exponentielle dans Lie Theory utilisé dans [START_REF] Fu | Analysis of unified error model and simulated parameters calibration for robotic machining based on Lie theory[END_REF]. Cependant, l'idée derrière toutes les approches proposées est de comparer le placement réel de l'effecteur final du robot avec le placement souhaité calculé à partir du modèle géométrique nominal. La différence entre les valeurs souhaitées et les valeurs réelles est utilisée pour mettre à jour en permanence le modèle géométrique du robot afin de surmonter les erreurs géométriques. L'analyse documentaire présentée par [START_REF] Wu | Geometric calibration of industrial robots using enhanced partial pose measurements and design of experiments[END_REF] peut servir de guide pour les lecteurs inclinés dans le domaine de l'étalonnage géométrique.

Bien que l'impact des erreurs géométriques du robot soit indéniable dans les applications basées sur la robotique, un coup d'oeil sur la littérature révèle que cette source d'erreur est principalement considérée comme le principal facteur contribuant à une situation de déchargement du robot. [START_REF] Schneider | Experimental investigation of sources of error in robot machining[END_REF] A la lumière des points mentionnés ci-dessus, dans ce chapitre de la thèse, l'effet des erreurs non géométriques est étudié pour révéler leur effet et leur contribution dans le mauvais positionnement du point central de l'outil robot (TCP).

Tout au long de ce chapitre, nous avons étudié l'effet des erreurs non géométriques sur la pose IR destinée à l'usinage. Le comportement de conformité de la structure du robot a été étudié. Il a été mentionné que le comportement de conformité pour les tâches de robot statiques et quasi-statiques peut être divisé en deux catégories différentes, à savoir l'effet d'auto-gravité qui est dû aux couples imposés aux articulations du robot pour résister au poids des maillons du robot et l'effet de charge externe qui est la contribution de la charge appliquée sur TCP aux flexions des articulations.

Au cours des expériences réalisées sur le comportement de conformité du robot, nous avons constaté des irrégularités cycliques qui ont été réparties tout au long de la trajectoire souhaitée. Ces erreurs cycliques étaient indépendantes de la trajectoire externe. En attribuant ces erreurs au désalignement et à l'usure des engrenages dans les transmissions articulées, nous avons concentré notre attention sur ce type d'erreur. Ce chapitre a ensuite été suivi d'une revue de la littérature sur les sources de cette erreur et plus généralement sur les imperfections du comportement du réseau de transport.

Pour observer l'effet des variations de température le long de la structure du robot et sa précision due à l'exécution des tâches, l'appareil du robot a été programmé pour déplacer continuellement ses articulations à vitesse maximale et l'effet de la chaleur produite par les actionneurs du robot a été rapporté en termes de changements de température en plusieurs points du corps et a augmenté la précision due à la dilatation thermique des liaisons. D'après l'analyse des données recueillies dans les simulations et expériences présentées dans ce chapitre, pour un bras de robot industriel Stäubli TX200 exécutant une tâche quasi statique avec une charge verticale externe de 500(N) (ce qui n'est pas une hypothèse irréaliste alors qu'une broche ordinaire pèse environ 40(kg)), la portion de contribution de chaque source d'erreur non géométrique étudiée est indiquée dans le graphique circulaire illustré à la La conclusion tirée de la Figure 2.10 est l'influence vivante du comportement de conformité sur le déplacement de pose du robot qui représente environ 87% (Autogravité et effet de charge externe). Selon la figure, l'effet du comportement de conformité est 7 fois supérieur à celui des autres sources combinées.

A.3 Chapitre 3

L'adoption de bras robotisés pour l'usinage présente à la fois des avantages et des défis. L'un des principaux défis des IRs utilisés pour les tâches d'usinage est leur imprécision par rapport à leurs alternatives conventionnelles. Comme nous l'avons conclu au chapitre 2, les principales sources d'erreurs non géométriques, par ordre d'importance de leur contribution, sont le comportement de conformité (résultant à la fois des effets de l'auto-gravité et des charges externes), les erreurs cycliques attribuées aux systèmes de transmission communs et les erreurs induites par la chaleur.

En raison de l'importance de ces sources d'erreur, ce chapitre est axé sur l'identification et la compensation des deux premières sources d'erreur non géométriques mentionnées, à savoir Compliance behaviour et Erreurs cycliques du système de transmission par articulation. Ces sources sont responsables d'environ 96% des erreurs de positionnement du robot TCP dans les conditions expérimentales décrites au chapitre 2.

A.3.1 Modélisation et compensation du comportement de conformité du robot

La structure en série des IRs, bien qu'elle offre une magnifique polyvalence pour différentes applications, induit un inconvénient majeur pour la précision du robot. L'imprécision due à la flexibilité des articulations et des maillons est à l'origine d'une grande imprécision de pose due à l'autogravité et aux forces d'usinage et reste le principal obstacle à l'utilisation généralisée des robots pour des applications d'usinage de précision.

Cette section est dédiée à la génération d'un modèle mathématique capable de calculer le couple, l'écart angulaire des joints et le déplacement du TCP dû à l'autogravité et aux forces externes basées sur la méthodologie VJM. Le modèle est subdivisé en trois parties : la première est l'explication de l'effet de la gravité sur chaque articulation pour une version robotisée de base vers le bas (montée au sol). La deuxième partie présente la modélisation du système d'équilibrage sur l'articulation numéro deux du robot et la dernière partie calcule le couple induit par rapport aux forces externes.

Effet de gravité de la structure du robot on peut calculer le couple appliqué sur chaque axe d'articulation causé par les forces gravitationnelles alors que la base du robot a été montée au sol en utilisant un algorithme récursif à partir du joint 6 et en appliquant l'équation d'équilibre sur chaque liaison. La valeur du couple imposé au joint i est déterminée par l'expression suivante :

   T i = T i+1 + -→ d i × m i -→ g + -→ r i × ∑ N j=i+1 m j -→ g , (i < 6) T i = -→ d i × m i -→ g , (i = 6) (A.1)
Où T i est le vecteur de couple appliqué sur l'articulation i, -→ d i et -→ r i sont des vecteurs reliant l'origine de la trame i au centre de gravité du lien i (CG i ) et l'origine de la trame i + 1. Et m i représente la masse du lien i. La figure A.4 représente les paramètres décrits pour calculer le couple imposé sur l'articulation 2 pour une configuration de robot arbitraire due à la force gravitationnelle.

Système d'équilibrage

Le robot Stäubli TX200 utilise un système d'équilibrage armé sur l'articulation 2 pour réduire l'effet du couple gravitationnel appliqué sur cette articulation. Ce système consiste en un ressort fixé à un point fixe de la liaison 1 à une extrémité et près de l'extrémité de la liaison 2 à l'autre extrémité (Figure A.5).

Pour observer l'effet de ce système d'équilibrage sur les performances du robot, nous avons simulé un mouvement du robot dans lequel l'articulation 2 tourne de -90(Deg) à +90(Deg) tout en ayant d'autres articulations à leur position zéro constamment sans charge externe appliquée sur l'effecteur final. La figure A.6 montre la différence de couple appliqué sur l'articulation 2 pour deux cas d'articulation non armée et armée avec système d'équilibrage en lignes bleue et rouge respectivement. Comme on peut le voir, le système d'équilibrage compense massivement l'effet d'auto-gravité du robot sur le couple appliqué sur l'articulation 2.

Effet d'un couple externe

Pour tenir compte de l'effet du couple externe résultant des forces d'usinage, il faut convertir le couple exercé en couples résultants appliqués aux articulations du robot. Bien que ces coupures provoquent à la fois des forces et des couples sur les articulations et les maillons du robot, mais en raison de l'hypothèse de la méthode VJM,

A.3.2 Erreurs cycliques du système de transmission articulaire

Il existe différentes approches pour surmonter les erreurs communes. Les performances des robots peuvent être améliorées en utilisant soit des moteurs à entraînement direct (comme dans certains robots SCARA), soit des engrenages de haute précision caractérisés par un très faible jeu, un rendement élevé, une faible inertie, une faible friction, une rigidité élevée et un faible poids ; soit en plaçant des codeurs haute résolution à la sortie des engrenages. Cependant, ces solutions ne peuvent pas être appliquées simultanément et augmenteraient considérablement le coût de fabrication d'un robot industriel [START_REF] Slamani | Characterization and experimental evaluation of gear transmission errors in an industrial robot[END_REF] Dans le laboratoire LISPEN, l'expérience de travail sur l'identification et la compensation du comportement des engrenages non linéaires remonte à 2012 à partir des travaux de Olabi et al., 2012. Bien que les auteurs n'aient pas présenté de résultats pratiques dans le document mentionné, le résultat final de la compensation des défaites non linéaires du système de transmission a été présenté dans Olabi, 2011. Là, l'auteur a identifié les défaites d'engrenages des articulations 2 et 3 d'un bras de robot Stäulbi RX 170BH et a compensé la trajectoire articulaire liée à une ligne droite dans le système de coordonnées cartésien. Dans cette partie, nous proposons une nouvelle méthode pour éviter le problème de la consommation de temps. Dans cette méthode, nous avons proposé une nouvelle approche pour définir la trajectoire du robot pour l'identification des erreurs de transmission : l'erreur de toutes les articulations souhaitées peut être identifiée en même temps, ce qui peut réduire considérablement le temps d'arrêt du robot. Nous ciblons la plage acceptable des fluctuations causées par les engrenages pour qu'elles se situent dans les limites de la répétabilité du manipulateur du robot.

Dans le cas général des robots 6 axes, en raison de l'effet de levier et des rapports de démultiplication plus élevés, l'effet des trois premières articulations est supposé être dominant dans l'erreur de positionnement du robot. Puisque la recherche de la meilleure trajectoire dans l'espace cartésien conduit éventuellement à des incréments variables de points discrets dans l'espace d'articulation, la trajectoire proposée est le résultat de la discrétion d'un certain intervalle de valeurs communes dans l'espace articulaire. La trajectoire est composée d'un certain nombre de poses dans lesquelles le robot doit s'arrêter et la coordonnée de l'effecteur final est mesurée. L'équation A.2 exprime le format général des variables communes de la trajectoire :

θ i = [θ 1 (0) + i * ∆θ 1 , θ 2 (0) + i * ∆θ 2 , θ 3 (0) + i * ∆θ 3 , 0, 0, 0] (A.2)
La figure A.8 représente le chemin cartésien résultant des variables communes mentionnées du robot TCP définies dans le système de coordonnées de base du robot. Pour annuler l'effet des autres sources d'erreur, on soustrait le polynôme ajusté du total des valeurs d'erreur. La partie résiduelle est ce que nous prétendons qui est essentiellement causée par le comportement non linéaire des engrenages de transmission. En raison de la nature périodique de cette erreur, nous avons choisi d'ajuster la somme des fonctions sinusiennes pour exprimer l'erreur pour une articulation donnée j :

∆θ nj = N ∑ k=1 a k sin(b k θ j + c k ) (A.3)
Pour une trajectoire désirée, principalement, les propriétés cartésiennes désirées devraient être converties dans l'espace articulaire en utilisant la cinématique inverse du robot et seulement dans cet espace on peut estimer et éliminer par la suite les erreurs non linéaires de la trajectoire induites par les engins par cette approche proposée. L'élimination de l'erreur n'est que la soustraction de l'erreur estimée des variables communes souhaitées.

La validité et la performance de la méthode proposée ont été examinées par l'exécution d'un test simple.nous avons programmé le robot pour qu'il suive une droite de 25(cm) avec le point de départ de X s = [500, 50, 0](mm) et le point final de X e = [750, 50, 0](mm). La trajectoire a été générée point par point avec des incréments de 0,2(mm) menant à 1250 points. Le robot subit une pause de 3 secondes sur chaque point pour permettre au traqueur d'enregistrer la position de l'effecteur final et aussi pour éviter toute vibration de la structure du robot. Les valeurs des coordonnées TCP résultantes de ce chemin sont présentées dans la Figure A.9 pour les variables communes compensées et non compensées.

A.3.3 Conclusion

Dans ce chapitre, des efforts ont été consacrés à l'identification et à la compensation du comportement de conformité de la structure du robot et des systèmes de transmission à articulation cyclique comme principaux facteurs contribuant aux erreurs de pose non géométrique du robot. En raison de l'influence considérablement plus grande du comportement de conformité sur l'inexactitude du robot, cette source a été étudiée au préalable dans la section 3.2.

Pour modéliser le comportement de conformité, nous avons utilisé l'approche VJM (Virtual Joint Method) pour générer un modèle élastostatique en fonction de sa performance convaincante présentée dans la littérature existante. Le modèle élastostatique généré implique à la fois l'auto-gravité et l'effet de force externe. La section 3.2.3 a été attribuée à l'effet d'auto-gravité dans laquelle la procédure de mise en oeuvre du modèle de masse du robot (Liens poids et coordonnées des centres de gravité), le couple gravitationnel imposé aux articulations du robot pour une configuration donnée et l'effet du système de balancement attaché au robot ont été expliqués en détail.

En raison de la nécessité de valeurs de rigidité articulaire dans le modèle élastostatique, des expériences ont été menées pour déterminer ces valeurs. La section 3.2.5 présente les procédures d'imposition des forces et les dispositifs de mesure utilisés pour chaque articulation. L'effet non linéaire du système d'équilibrage sur la rigidité de l'articulation 2 a également été modélisé dans cette section à l'aide d'un ressort de torsion à rigidité variable. Cependant, d'après les enquêtes, son effet était mineur par rapport à la rigidité de l'articulation 2 et a donc été négligé pour d'autres considérations.

Pour vérifier le modèle élastostatique généré, des expériences ont été menées en utilisant l'option de compensation de charge offerte par le fournisseur du bras robotique Stäubli company. Ce module de compensation étant conçu pour contrer l'auto-gravité et les dispositifs supplémentaires fixés à la bride du robot (Broche ou générateur d'ondes ultrasonores par exemple), il a permis de valider notre modèle élastostatique. Le processus de vérification a révélé la validité du modèle élastostatique. Cependant, les valeurs de rigidité obtenues dans la procédure expérimentale ont été remplacées par la paire fournie par Stäubli en raison des simulations plus précises utilisant cette dernière.

Dans la section 3.3, les erreurs cycliques du système de transmission étaient concentrées. Une nouvelle méthode a été proposée pour modéliser mathématiquement cette erreur sur la base d'un chemin optimal pour réduire le temps de l'identification par rapport à la méthode existante dans la littérature. Une méthode de compensation hors ligne a été proposée sur la base de l'estimation et de l'élimination des erreurs cycliques pour une trajectoire requise dans l'espace de travail du robot. Des expériences sur un bras robotique de laboratoire (Stäubli TX90) ont montré que le procédé d'identification et de compensation proposé était capable de réduire l'amplitude de ces fluctuations par deux fois la répétabilité du robot. Cette méthode peut être utilisée pour augmenter la précision de la trajectoire des manipulateurs du robot, en particulier ceux qui sont utilisés pour des applications d'usinage où les déviations de trajectoire provoquent des imperfections de qualité indésirables. L'angle de chanfreinage en régime permanent est appelé l'angle que la rampe génère avec la surface inférieure dans la zone de coupe en ligne droite (β dans la figure A.11). Ce paramètre étant notre principale préoccupation en raison de l'application du Creneau, cet angle a été mesuré à l'aide d'un appareil et d'un procédé plus précis. Nous avons utilisé un dispositif à bras pour mesurer cet angle. A chaque fois, la pièce était fixée sur une table en marbre dont la planéité est de l'ordre du micron. Les pièces ont été montées sur la table en utilisant deux étaux comme supports. En touchant la table en marbre avec l'effecteur final du Faro, nous avons construit le plan de table. Ce processus a été automatisé par AutoDesk Logiciels( Figure A.12). Pour des applications telles que la simulation du comportement de conformité de la structure du robot dans des conditions de chargement du processus d'usinage par ultrasons en nid d'abeilles, on peut utiliser soit le vecteur de force défini initialement dans le système de coordonnées Kistler (en supposant que la transformation entre le Kistler et le bâti du robot soit, de manière plus pratique, des forces d'usinage définies dans le bâti de l'outil qui est équivalent au bâti TCP du robot).

A.4 Chapitre 4 A.4.1 Erreurs d'usinage dans le processus de chanfreinage en nid d'abeilles

L'expression de force générique proposée pour n'importe quelle direction du cadre de la pièce et du cadre de l'outil est supposée avoir le format suivant : F c cos(β) + F s sin(β) + F 0 (A.4) Pour trouver les valeurs appropriées de F c ,F s et F 0 Least Square, la méthode a été appliquée aux données expérimentales obtenues dans la procédure de mesure de force décrite. Les expressions suivantes estiment la force d'usinage pour un angle de chanfreinage donné de β exprimée dans le système de coordonnées de la table Kistler présenté dans la Comme nous utilisons le modèle élastostatique du bras du robot développé au chapitre 3, les simulations peuvent prendre en compte les trois effets suivants : 1. Effet gravitationnel des éléments du robot (effet auto-gravité), 2. masse de l'ensemble Les résultats de la simulation étaient inexacts pour ce qui est de décrire les erreurs d'usinage prévues. Dans la section suivante, nous allons générer un modèle décrivant le comportement de conformité de l'outil d'usinage (qui est la lame de coupe à ultrasons) à intégrer dans une étape ultérieure dans le comportement de conformité total du système robot-outil.

A.4.4 Comportement de conformité des couteaux

Pour générer un modèle de conformité du couteau, nous utilisons l'approche SSM (Matrice de rigidité structurelle). Cette matrice est capable d'estimer la déformation du couteau pour une condition de charge donnée. La matrice de rigidité structurelle (SSM) est définie comme la matrice de rigidité reliant le déplacement de la pointe du couteau à la charge appliquée sur le couteau. La condition de chargement exprimée dans la Figure A.17 est utilisée pour générer le SSM. Comme on peut facilement conclure en comparant différents éléments du déplacement de la pointe du couteau, le comportement de conformité du couteau dans la direction z (direction verticale de la structure du couteau) est de loin supérieur au reste. Le déplacement dans le sens longitudinal est inférieur à un micron. Le déplacement de la direction latérale est environ 20 fois inférieur à la direction verticale. Il faut noter que non seulement les valeurs de rigidité de ces deux directions sont considérablement plus grandes que la direction verticale, mais aussi les forces de chanfreinage dans ces directions sont toujours plus petites que la composante verticale.

A.4.5 Comportement de conformité des robots et couteaux

L'approche pour vérifier l'effet de chaque système (couteau et structure du robot), est d'utiliser une articulation virtuelle avec une rigidité rotationnelle équivalente pour représenter l'effet de conformité de l'outil fixé au robot. Ainsi, la modélisation de la conformité de la structure du robot à 7 axes résultante serait un format généralisé de celui appliqué sur le robot à 6 axes. dépend de la configuration, le comportement de conformité de la structure du robotoutil sera mesuré dans une configuration dans laquelle le couple maximal est imposé à l'articulation du robot (Figure A.18). Le déplacement du TCP, comme mentionné précédemment dans le chapitre 3, peut être calculé par l'expression suivante :

δP t ω t = (J 6×7 × [K] -1 × J T 7×6 ) × F (A.9)
Le remplacement des valeurs expérimentales de la composante de conformité robot-couteau dans la direction z, révèle le pourcentage suivant de contribution pour la structure du robot et l'élément couteau:

δz tcp f z = 561.3 × k -1 2 + 201.4 × k -1 3 + 38.3 × k -1 5 11%
+ 0.49 × k -1 eq 89% (A.10)

Les résultats présentés par l'équation A.10 indiquent que l'effet du comportement de conformité de la structure du couteau est dominant par rapport à l'augmentation de plus de 8 fois de la portion de la structure du robot.

A.4.6 Modélisation de la condition de chargement

Parmi toutes les fonctions mathématiques pour décrire la condition de distribution qui peut satisfaire toutes les hypothèses mentionnées ci-dessus, le choix d'une fonction polynomiale générale serait une option plus sage en raison de leur vaste utilisation et de la simplicité des opérations algébriques. La fonction de charge distribuée ω(x) de l'ordre n peut être exprimée par l'expression générale suivante :

A.4.7 Procédure de Compensation

La procédure d'indemnisation proposée dans cette partie est basée sur une approche hors ligne. La stratégie utilise le modèle de force développé dans la direction verticale du couteau pour estimer les forces d'usinage appliquées sur le couteau et utilise par conséquent le modèle de comportement de conformité généré précédemment pour estimer l'erreur de chanfreinage basée sur les forces estimées pour les caractéristiques de chanfreinage requises. Le processus vise à proposer un angle de chanfreinage modifié (β m ) pour un angle de chanfreinage souhaité correspondant de β. Le calcul de la valeur de β m correspondant à la valeur souhaitée de β consiste en un algorithme itératif. Dans la première étape, un angle de chanfreinage est proposé (pour la première itération, la valeur souhaitée β est utilisée). Les forces de coupe pour la valeur de chanfreinage proposée sont estimées à l'aide du modèle de force d'usinage vertical. La déformation de la structure du couteau due à la force estimée est calculée en utilisant le modèle de couteau qui est appelé θ dev ci-après. Cette déformation peut être estimée en utilisant des conditions de charge concentrée ou répartie. L'angle de chanfreinage estimé (β = β m + θ dev ) est alors comparé à la valeur de l'angle de chanfreinage souhaité β. L'erreur de chanfreinage de chaque itération er est définie comme er = ββ . Si l'erreur calculée de l'itération (er) est supérieure à la valeur prédéfinie de , le processus continue avec la valeur actualisée de β m calculée comme β m = β mer. La boucle se termine lorsque la valeur de l'erreur d'itération est inférieure à la valeur de . La figure A.21 présente l'erreur de chanfreinage dans le processus de coupe en nid d'abeille sur le matériau dur en utilisant la stratégie de compensation avec l'erreur non compensée mesurée dans la première session.

A.4.8 Conclusion

Dans ce chapitre, nous avons procédé à l'étude d'un processus d'usinage robotisé dans lequel le comportement de conformité de l'outil d'usinage ne pouvait être négligé par rapport à celui du bras robotique industriel opérationnel dans le total des erreurs de fabrication observées. Le procédé d'usinage adopté pour ce chapitre était le chanfreinage de matériaux en nid d'abeille aramide utilisant la technologie ultrasonore basée sur les exigences et les intérêts de la société Le Creneau.

Pour poursuivre nos recherches afin de trouver le principal responsable des erreurs d'usinage lors de la découpe de matériaux en nid d'abeilles à l'aide de la technologie de découpe par ultrasons avec un bras de robot comme opérateur et un couteau triangulaire comme outil de coupe, nous avions besoin des forces d'usinage impliquées dans ce processus. En raison de l'absence d'un modèle précis des forces d'usinage dans la littérature, une approche empirique a été utilisée pour générer un Lors de notre première tentative, nous avons simulé l'effet de conformité de la structure du robot. Les résultats des simulations sur le comportement de conformité du robot étaient imprécis et incapables d'expliquer les valeurs d'erreur mesurées dans les pièces obtenues expérimentalement. L'ordre de désalignement et d'erreur de déplacement dans l'orientation et le positionnement TCP étant négligeable, le comportement de conformité de la lame de coupe a été intégré.

Dans une étape ultérieure, pour une modélisation plus précise du comportement du couteau, la charge appliquée sur le plan latéral du couteau a été modélisée par une force non uniformément répartie. La force non uniformément répartie était exprimée par un polynôme d'ordre n, dont la valeur était déterminée en appliquant la méthode des moindres carrés aux données expérimentales. le chapitre s'est poursuivi par la proposition d'une méthode de compensation. La méthode était basée sur une procédure de calcul itérative capable de prédire la valeur initiale du chanfreinage correspondant à une valeur de chanfreinage désirée en tenant compte du modèle de force et par conséquent de la déformation de la structure de la lame due à sa force d'usinage correspondante. Les résultats de cette méthode de compensation sont prometteurs.

A.5 Chapitre 5

Le pourcentage d'usinage robotisé augmente rapidement dans les différents secteurs industriels. L'augmentation de la demande de manipulateurs de robots année après année est évidente pour les avantages des bras de robots industriels dans de nombreuses applications [START_REF] Verl | Robots in machining[END_REF]. Les robots industriels sont connus pour leurs prix comparativement plus bas par rapport aux appareils conventionnels tout en offrant des avantages attrayants tels que des espaces de travail plus vastes, la reconfigurabilité et la dextérité. De nos jours, les manipulateurs de robots jouent un rôle dans différents processus d'usinage tels que le fraisage, le meulage, l'ébavurage et le perçage pour ne citer que quelques exemples [START_REF] Ji | Industrial robotic machining: a review[END_REF].

Cependant, le nombre d'inconvénients des manipulateurs de robots empêche leur taux d'utilisation plus élevé et leur contribution au processus de fabrication de précision. Les défauts inhérents aux bras robotiques sont parmi les principaux facteurs d'erreurs et d'imprécisions d'usinage. Les faibles valeurs de rigidité des articulations et les structures en série des bras robotiques industriels présentent de grands défis, tels qu'un comportement de conformité considérablement plus élevé par rapport aux dispositifs d'usinage conventionnels, un mauvais comportement dynamique responsable de bavardages indésirables [START_REF] Huynh | Modal updating of a 6-axis robot for milling application[END_REF], des vibrations et des charges de travail accrues. [START_REF] Nelson | Locating assembly tasks in a manipulator's workspace[END_REF] a ouvert un nouveau chapitre dans le cours de la recherche robotique et a ouvert la voie à une voie qui a été un point de départ pour de nombreux futurs chercheurs. Pour obtenir une performance de batterie d'une application robotique, ils ont proposé de modifier la configuration de la cellule robotisée. La soi-disant conception de cellule de travail s'est développée pour optimiser de nombreux indices de performance d'une cellule de travail robotique tels que le temps de cycle [START_REF] Spensieri | Optimal robot placement for tasks execution[END_REF] et la consommation énergétique [START_REF] Libor | Optimizing energy consumption of robotic cells by a Branch & Bound algorithm[END_REF] par des efforts supplémentaires des experts. L'approche de conception de la cellule de travail est entrée dans le problème de réduction des erreurs d'usinage du robot en tant que sous-catégorie des stratégies hors ligne. Les chercheurs ont commencé à bénéficier soit de la redondance introduite par les bras robotiques 6 axes dédiés aux tâches d'usinage 5 axes [START_REF] Xiong | Stiffness-based pose optimization of an industrial robot for five-axis milling[END_REF], soit des propriétés du robot dépendantes de la configuration comme la rigidité locale ou les fréquences propres [START_REF] Celikag | Chatter Suppression in Robotic Milling by Control of Configuration Dependent Dynamics[END_REF][START_REF] Chen | Posture-dependent stability prediction of a milling industrial robot based on inverse distance weighted method: Paper ID: 1222[END_REF][START_REF] Denkena | Design and optimization of a machining robot[END_REF] pour améliorer l'usinage.

Selon la littérature, l'une des approches les plus efficaces pour augmenter les performances du robot dans une tâche d'usinage est de modifier le positionnement relatif de la pièce par rapport au bras robotique. D'après Caro et al., 2013, pour un Kuka IR exécutant une opération de fraisage sur un bloc d'aluminium, l'erreur de fabrication résultant du comportement de conformité du robot est réduite d'environ 80 % dans le placement optimal par rapport au pire placement.

A.5.1 Énoncé du problème

Notre effort dans ce travail est consacré à l'amélioration de la qualité d'un processus d'ébavurage en tant que processus d'usinage final affecté par les erreurs périodiques de transmission conjointe causant cyclique dans le déplacement de l'effecteur final du robot et le comportement de conformité du robot. Ces erreurs sont principalement dues à des erreurs de fabrication et d'assemblage de la boîte de vitesses qui entraînent des imperfections géométriques du système de transmission dans l'état souhaité (Figure A.22). Selon les valeurs en Karpat, Bahtiyar, and De ger, 2012 entraînées expérimentalement, la composante de force tangentielle maximale pendant une rotation de l'outil de coupe varie entre 100(N) et 140(N) pour une avance de 350(mm/min) en fonction de la direction des fibres. Ces valeurs pour la composante de force radiale sont de 180(N) à 300(N). Dans ce travail, les valeurs maximales des deux composantes (F t = 140(N),F r = 300(N)) sont utilisées pour représenter les forces de coupe.

A.5.4 Méthode d'optimisation

Un bref coup d'oeil sur le problème de l'optimisation de la conception des cellules de travail dans la littérature révèle qu'une main pleine de méthodes d'optimisation a été utilisée pour cet objectif. Parmi toutes les méthodes, l'algorithme génétique (GA) [START_REF] Caro | Workpiece placement optimization for machining operations with a KUKA KR270-2 robot[END_REF] et l'essaim de particules [START_REF] Doan | Optimal robot placement with consideration of redundancy problem for wrist-partitioned 6R articulated robots[END_REF] sont les plus populaires en raison de la nature non linéaire du problème. Cependant, dans ce travail, nous proposons d'aborder le problème d'optimisation avec un algorithme de recherche de modèle libre modifié pour découvrir l'efficacité de cette méthode dans la conception de cellules de travail robotiques.

La principale modification appliquée à l'algorithme de recherche de motif dans ce travail est la géométrie de la grille. Comme mentionné précédemment, dans l'approche conventionnelle de construction de grille, chaque point de la grille diffère seulement la valeur d'une composante de sa coordonnée dans l'espace variable du point courant. Comme chaque variable a trois possibilités différentes de changement (dans le sens croissant, constant et décroissant), une grille pleine grandeur (employée dans ce travail) se compose de 3 n nombre de points.

A.5.5 Résultats de la simulation

La méthode d'optimisation expliquée précédemment a été utilisée pour minimiser la fonction de coût définie précédemment. Les limites supérieure et inférieure utilisées pour les quatre variables de conception dans le processus d'optimisation ont été choisies pour couvrir presque tous les points accessibles de la zone avant du robot.

Pour comparer les erreurs cycliques de conformité et de transmission observées dans la trajectoire de coupe de la pièce, dans le placement optimal de la pièce et dans un placement intuitif, un point dans l'espace de travail du robot est choisi pour représenter un placement intuitif probable de la pièce.Figure A.25 montre un déplacement exagéré de l'outil d'usinage dû à une erreur de transmission cyclique le long de la circonférence de la trajectoire de coupe.

A.5.6 Conclusion

La conception de la cellule de travail dans ce chapitre a été abordée en remplaçant la pièce à usiner par rapport au châssis de base du robot, ce qui permet d'éviter le démontage de la structure du robot dès son placement initial. Le processus d'usinage ciblé était une tâche d'ébavurage en CFRP et la pièce à usiner était présentée par une forme circulaire. L'objectif était de réduire les erreurs dues à la fois au comportement élastostatique et aux erreurs périodiques de transmission articulaire. La fonction de coût proposée a été conçue en utilisant la variance de la normale sur la direction de l'avance simulée à partir des modèles d'erreur générés précédemment. Les contraintes cinématiques du système robotique, à savoir les limites articulaires et l'évitement de la singularité, ont été mises en oeuvre par une méthode de pénalisation. Une nouvelle méthode a été proposée pour mettre en oeuvre l'évitement de la singularité applicable à toute structure de robot partitionné au poignet. La méthode d'optimisation adoptée pour ce problème était une modification de l'algorithme de recherche de modèle libre dans lequel une grille pleine grandeur était utilisée pour faciliter la prise de décision sur la direction du progrès. La performance et l'efficacité de la procédure d'optimisation proposée sont validées par les résultats de simulation.

A.6 Conclusion de la thèse

Le concept de l'industrie manufacturière doit inévitablement s'adapter à la croissance rapide et exigeante des signes récents de progrès dans d'autres domaines de l'ingénierie. L'évolution rapide de la complexité et de la nature dynamique des progrès technologiques ne peut continuer à prospérer qu'en s'appuyant sur un concept de production flexible, prêt à se réformer et à s'adapter à de courts changements observés à faible coût. Les bras de robots industriels (IR) ont démontré des avantages attrayants pour mener cette percée révolutionnaire.

Bien que les infrarouges offrent d'énormes avantages tels que des espaces de travail étendus et la reconfigurabilité, ils présentent toutefois des inconvénients pour les traitements et compensations ultérieurs, afin de répondre aux normes strictes de précision exigées dans les secteurs de haute technologie.

La contribution de chaque source d'erreur non géométrique a été examinée au chapitre 2. Les examens ont été effectués en conditions chargées et déchargées sur un robot d'étude de cas (Stäubli TX200). Les résultats ont montré que les erreurs résultant du comportement de conformité de la structure du robot sont dominantes. Cette source d'erreur a été divisée en deux catégories, à savoir l'effet d'auto-gravité et l'effet de charge externe. La première catégorie impose des valeurs d'erreurs plus élevées que la seconde.

En plus de la source d'erreurs de conformité, des inexactitudes provenant d'erreurs périodiques des systèmes de transmission conjoints et du comportement thermique des robots ont été évaluées. D'après les résultats expérimentaux, la contribution de chaque source d'erreur aux erreurs non géométriques pour une condition de charge quasi statique était : 51%, 36%, 9% et 4% pour l'effet de gravité, l'effet de charge externe, les erreurs de transmission non linéaire des joints et le comportement thermique respectivement pour les conditions d'essai décrites dans ce chapitre.

Le chapitre 3 était consacré à l'identification et à la compensation de deux min contributeurs aux erreurs non géométriques, à savoir le comportement de conformité et les erreurs de transmission communes. Pour surmonter le comportement de conformité, un modèle élastostatique capable de comptabiliser les couples a été développé à partir de trois sources principales d'imposition de couple : la gravité (appelée effet auto-gravité dans la thèse), les charges externes et le système autoéquilibrant attaché au robot. Les valeurs de rigidité des articulations du robot ont été déterminées expérimentalement. La performance du modèle élastostatique généré a été validée expérimentalement. Bien que les expériences aient révélé une efficacité acceptable du modèle, les valeurs de rigidité des articulations obtenues auprès du fournisseur du robot ont donné de meilleurs résultats.

Pour compenser l'effet des erreurs non linéaires de transmission articulaire, un modèle mathématique de ce type d'erreur était nécessaire. Une nouvelle méthode a été proposée pour identifier cette source d'erreur basée sur la construction d'un chemin cartésien optimal. Le trajet a été conçu pour permettre à l'opérateur d'identifier simultanément le modèle d'erreur pour toutes les articulations ciblées. Par conséquent, le principal avantage de cette méthode était son temps d'exécution considérablement plus court par rapport aux méthodes disponibles dans la littérature. Une procédure de compensation hors ligne a été présentée avec laquelle la précision de la trajectoire d'une ligne droite a augmenté de 40 % en fonction des données expérimentales obtenues. Le problème a d'abord été abordé en développant un modèle de force basé sur des valeurs expérimentales de force pour le processus de chanfreinage des matériaux en nid d'abeille. Le comportement élastostatique de la structure du robot a été simulé à l'aide de forces d'interaction outil-pièce à commande empirique. L'inadéquation des résultats de simulation avec les valeurs d'erreur réelles nous a incités à inclure le comportement de conformité de la lame triangulaire (outil d'usinage) dans le modèle de conformité de l'ensemble du robot-outil.

Pour intégrer la conformité de l'outil dans le processus, la rigidité de l'outil a d'abord été modélisée en utilisant la méthode SSM (Structural stiffness matrix). Le l'adoption d'une méthode d'identification de la rigidité dans laquelle les déformations des liaisons du robot sont éliminées dans le processus de mesure (telle que la procédure d'identification de la rigidité dynamique) peut mieux servir.

Comme il a été conclu dans le chapitre 4 que la structure du couteau ultrasonique est à l'origine de l'effet dominant sur les erreurs géométriques observées dans les pièces en nid d'abeille, une approche efficace consiste à améliorer le comportement de conformité de l'outil. Plus tard, il est possible de se concentrer sur la reconception de la géométrie du couteau en considérant l'amplitude de la contrainte dynamique le long du corps du couteau. Une contrainte dynamique élevée entraînerait une faible performance de fatigue de l'outil.

La nouvelle approche proposée au chapitre 5 pour mettre en oeuvre les critères d'évitement de singularité dans le problème d'optimisation consiste en trois coefficients représentant trois conditions de singularité différentes dans un bras robotisé à poignet portionné. La faisabilité de l'utilisation de la méthode présentée pour développer une nouvelle méthode pour décrire la manipulabilité d'un robot orienté vers le poignet peut être un sujet de recherche intéressant.

  A) A Conventional Gantry Structured CNC machine (B) Stäubli TX200 Robot Arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Thesis Positioning and Configuration . . . . . . . . . . . . . . . . . . . 2.1 Classification of Mechanical Robot Pose Errors . . . . . . . . . . . . . . 2.2 Configuration of a robot joint armed with external (secondary) encoder (Möller et al., 2017) . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Robot configuration and results for compliance behaviour of Stäubli TX200 robot arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Cyclic TCP position error along a straight line for loaded and unloaded conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Location of the chosen points for temperature measurements . . . . . . 2.6 (A) Laser tracker used for position measurements. (B) Thermometer utilized for temperature measurements . . . . . . . . . . . . . . . . . . 2.7 Temperature rising trend for different points on robot armor as function of warming time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Location of the chosen spots for repeatability measurements with respect to robot coordinate system . . . . . . . . . . . . . . . . . . . . . . 2.9 Evolution of position Repeatability of robot as function of Warming Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Contribution of each non geometrical error source to the robot pose misplacement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Stäubli VS DH frames of Stäubli TX200 robot arm . . . . . . . . . . . . 3.2 Computation of Self gravity torque imposed on joint 2 . . . . . . . . . . 3.3 schematic of balancing system attached to joint 2 of Stäubli TX200 robot arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Force imposition of balancing system on joint 2 . . . . . . . . . . . . . . 3.5 Effect of Balancing system attached to joint 2 on the imposed torque . . 3.6 Imposed torque on joints due to external load . . . . . . . . . . . . . . . 3.7 Illustration of measurement details of axis 1 stiffness identification using Faro arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Force imposition method for stiffness identification of joint 1 . . . . . . 3.9 Circle construction for joint 1 stiffness identification . . . . . . . . . . . 3.10 Measured and Fitted data for Stiffness identification of axis 1 . . . . . . 3.11 Configuration of the robot and force imposition system for joint 2 stiffness identification process . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Experimental data and fitted model to estimate stiffness of joint 2 . . . 3.13 Stiffness configuration of joint 2 . . . . . . . . . . . . . . . . . . . . . . . 3.14 Evolution of the non linear stiffness of robot balancing system attached to joint 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 Configuration of the robot and force imposition system for joint 3 stiffness identification process . . . . . . . . . . . . . . . . . . . . . . . . . . x 3.16 Experimental data and fitted model to estimate stiffness of joint 3 (Left) and 2 (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Robot Configuration and introduced frame attached to link 4 for joint 5 stiffness identification . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Measuring Process of Joint 5 Stiffness ID using Faro Arm . . . . . . . . 3.19 Plane construction for Stiffness ID of Joint 5 . . . . . . . . . . . . . . . . 3.20 experimental vs fitted data for joint 5 stiffness identification in different configurations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Illustration of Desired, Actual and Compensated path for compensation algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.22 Measurement Setup for Elastostatic Model verification . . . . . . . . . 3.23 TCP path for different test conditions of Elastostatic Model verification 3.24 Effect of backlash in the experimental path . . . . . . . . . . . . . . . . 3.25 Experimental and fitted values for test 1 and 2 . . . . . . . . . . . . . . 3.26 Experimental error values E for the entire rectangle on the left. Torques resulted from external force applied on different joint of robot arm on the right . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.27 Contour of estimation error value Es for a range of K 2 and K 3 (Dimension of the ES is (mm)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.28 Experimental Vs simulated positioning error values of TCP for the rectangle experiment using different joint stiffness values . . . . . . . . 3.29 Effect of compensation of Nonlinear behavior of transmission in a straight line in Olabi, 2011. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.30 Proposed Trajectory for Gear Transmission identification in Olabi et al., 2012 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.31 Diagram of the Proposed methodology for gear error Identification . . 3.32 Experimental setup and Measurement device for identification of gear transmission errors executed in LISPEN laboratory . . . . . . . . . . . 3.33 TCP path of proposed robot movement for gear error identification . . 3.34 Total joint error values (∆θ j ) in blue along with ∆θ j or part of error in Red for joint j (j : 1, 2, 3) . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.35 Periodic Gear Errors(∆θ jn ) (in Blue) and Fitted Sinus Functions (in Red) 3.36 Gear Error Compensation Procedure . . . . . . . . . . . . . . . . . . . 3.37 Experimental results for Compensated (in red) and Uncompensated path (in blue) for a straight line . . . . . . . . . . . . . . . . . . . . . . . 4.1 General Schematic of a robot-involved machining compliance behaviour . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Different elements combined in the ultrasonic cutting process in Le Creneau Industriel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Different Cutting Tools Offered by Le Creneau. Ultrasonic Disks are shown at the top and Ultrasonic Knives are at the bottom of the figure 4.4 Typical Honeycomb Panel End Design[START_REF] Soovere | Random vibration analysis of stiffened honeycomb panels with bevelededges[END_REF] . . . . . . . . . . 4.5 Overall configuration of a Chamfer . . . . . . . . . . . . . . . . . . . . . 4.6 General schematic of chamfering operation of Honeycomb work-piece 4.7 Example of a workpiece after four chamfering process on each side . . 4.8 A resulted workpiece from the test session . . . . . . . . . . . . . . . . 4.9 General schematic of chamfering operation of Honeycomb work-piece 4.10 General schematic of chamfering operation of Honeycomb work-piece 4.11 Penetration of knife in the protective moquette . . . . . . . . . . . . . . 4.12 Placement of the force sensor regarding robot base . . . . . . . . . . . .

xi 4 .

 4 13 Coordinate system of the Kistler table on the left [Note that the XY plane coincides the surface of the Kistler ]. On the right: An illustration of a work-piece mounted on the Kistler . . . . . . . . . . . . . . . 4.14 Hard and Soft honeycomb materials used in the tests (on the left and right respectively ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Chamfering test configuration and orientation of the end-effector (Cutting Knife) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Measuring the entrance angle using a Goniometer . . . . . . . . . . . . 4.17 Installation of the Workpiece for measuring the chamfering angle in Steady State zone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.18 Chamfering Angle error in both Soft and Hard Honeycomb materials for different desired chamfering angle values . . . . . . . . . . . . . . . 4.19 Displacement Errors as function of chamfering Angle of hard honeycomb material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.20 Ultrasonic Cutting Process Forces for a Chamfering operation for β = 50[deg] in hard honeycomb material . . . . . . . . . . . . . . . . . . . . 4.21 Machining Forces in Hard and Soft Honeycomb chamfering process defined in Kistler sensor coordinate system . . . . . . . . . . . . . . . . 4.22 Knife frame (in red) with respect to the Kistler frame (in black) . . . . . 4.23 Machining Forces in Hard and Soft Honeycomb chamfering process defined in Knife coordinate system . . . . . . . . . . . . . . . . . . . . .4.24 Placement of the Nida Work-piece regarding to the robot . . . . . . . . 4.25 Machining head assembly for the bench mark experiment . . . . . . . 4.26 Dimensions of the cutting knife placement and orientation with respect to the robot flange frame . . . . . . . . . . . . . . . . . . . . . . . . 4.27 Evolution of CAM generated Euler angle components for TCP frame orientation during the chamfering process . . . . . . . . . . . . . . . . . 4.28 Knife and TCP frame rotation regarding Spindle rotation angle . . . . 4.29 Joint Variables for the bench mark test . . . . . . . . . . . . . . . . . . . 4.30 Cutting Forces of honeycomb material on each edge in presence of robot base frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.31 TCP Displacement (on the left) and deviation (on the right) resulted form simulation analysis for the described three loading conditions . . 4.32 Joint Deviation in NO compensation mode on the left. Software generated displacement schema on the right . . . . . . . . . . . . . . . . . . 4.33 Measurement Strategy for the effect of applied force on the deflection in the supporting table . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.34 Dimensions of cutting knife involved in ultrasonic chamfering process 4.35 Illustration of different knife structure directions . . . . . . . . . . . . . 4.36 Deformation of the beam model under a concentrated force in vertical direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.37 Cross Section of the actual and modeled Knife geometry for vertical loading condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.38 Moment imposed by a concentrated vertical force . . . . . . . . . . . . 4.39 Deformation of the beam model under a concentrated force in Lateral direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.40 Cross section of the knife for lateral bending deformation . . . . . . . . 4.41 Deformation of the beam model under a concentrated force in Longitudinal direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.42 Axial stress distribution of an element along the longitudinal axis . . . 4.43 Experimental setup along with loading and measurement schematic .

xii 4 .

 4 44 Experimentally driven with model predicted deformation of the knife structure under several vertical concentrated force condition . . . . . . 4.45 Virtual joint model of the knife structure for vertical compliance effect 4.46 Developed 7 axes robot structure for robot-tool compliance analysis . . 4.47 Different Load distribution conditions . . . . . . . . . . . . . . . . . . . 4.48 Bending moment in an arbitrary cross section in the right part . . . . . 4.49 Value of Error function (Er(n)) for different Loading Order (n) values . 4.50 Deviation and displacement error estimation of the knife structure using different loading order parameter (n) . . . . . . . . . . . . . . . . . 4.51 Offline Chamfering Angle Compensation Logic . . . . . . . . . . . . . 4.52 Modified chamfering angle values (β m ) for a range of desired chamfering angle (β) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.53 Uncompensated and compesated chamfering angle error values . . . . 4.54 Vertical Machining Forces and Chamfering angle errors for the three experimental sessions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Schematic of workcell for trimming process . . . . . . . . . . . . . . . 5.2 Implementation of Joint Transmission System errors in the optimization problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Illustration of Machining Forces . . . . . . . . . . . . . . . . . . . . . . . 5.4 Illustration of desired and actual path for Cost Function definition . . 5.5 Examples of two actual paths with high error variance (A) and low error variance (B) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Representation of discretization parameter (∆d) used for path definition124 5.7 Effect of discretization parameter (∆d) on the cost function evaluation error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Defining variable r for shoulder singularity avoidance . . . . . . . . . . 5.9 Contour of distance between wrist intersection point and vertical axis of first joint of the robot (r(θ 2 , θ 3 )) . . . . . . . . . . . . . . . . . . . . . . 5.10 Comparison of computation time in a jacobian based method and proposed method for singularity avoidance . . . . . . . . . . . . . . . . . . 5.11 Illustration of conventional and full-sized grid structure for free pattern search optimization method . . . . . . . . . . . . . . . . . . . . . . 5.12 Optimization Flow chart . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.13 Evaluation of Cost function for a plane in workspace . . . . . . . . . . 5.14 Effect of Frame rotation in evaluation of cost function . . . . . . . . . . 5.15 Final Result Of Work-piece Placement . . . . . . . . . . . . . . . . . . . 5.16 Final Result Of Work-piece Placement . . . . . . . . . . . . . . . . . . . A.1 (A) Une machine CNC conventionnelle à portique structuré (B) Stäubli TX200 Bras de robot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Positionnement et configuration de la thèse . . . . . . . . . . . . . . . . A.3 Contribution de chaque source d'erreur non géométrique à l'erreur de positionnement du robot . . . . . . . . . . . . . . . . . . . . . . . . . . . A.4 Calcul du couple d'auto-gravité imposé à l'articulation 2 . . . . . . . . A.5 schematic of balancing system attached to joint 2 of Stäubli TX200 robot arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.6 Effect of Balancing system attached to joint 2 on the imposed torque . . A.7 Couple de serrage imposé aux articulations sous l'effet d'une charge extérieure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii A.8 Chemin TCP du mouvement du robot proposé pour l'identification d'erreur d'engrenage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.9 Résultats expérimentaux pour les trajectoires compensées (en rouge) et non compensées (en bleu) pour une ligne droite . . . . . . . . . . . . A.10 Exemple d'une pièce après quatre chanfreins de chaque côté . . . . . . A.11 Schéma général de l'opération de chanfreinage de la pièce en nid d'abeilles157 A.12 Montage de la pièce à usiner pour la mesure de l'angle de chanfreinage dans la zone en régime permanent . . . . . . . . . . . . . . . . . . . . . A.13 Erreur d'angle de chanfreinage dans les matériaux en nid d'abeilles souples et durs pour différentes valeurs d'angle de chanfreinage souhaitées158 A.14 Cadre du couteau (en rouge) par rapport au cadre Kistler (en noir) . . A.15 Machining head assembly for the bench mark experiment . . . . . . . A.16 Le déplacement TCP (à gauche) et l'écart (à droite) ont résulté de l'analyse de simulation pour les trois conditions de charge décrites . . A.17 Illustration des différentes directions de la structure des couteaux . . . A.18 Développement d'une structure robotisée 7 axes pour l'analyse de conformité des robots-outils . . . . . . . . . . . . . . . . . . . . . . . . . A.19 Différentes conditions de répartition de la charge . . . . . . . . . . . . . A.20 Estimation des déviations et des erreurs de déplacement de la structure de la lame à l'aide de différents paramètres d'ordre de chargement (n) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.21 Valeurs d'erreur d'angle de chanfreinage non compensées et compensées . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.22 Schéma de la cellule de travail pour le processus d'ébavurage . . . . . A.23 Illustration des forces d'usinage . . . . . . . . . . . . . . . . . . . . . . . A.24 Illustration du chemin souhaité et réel pour la définition de la fonction de coût . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.25 Résultat final du placement de la pièce à usiner . . . . . . . . . . . . . . xv List of Tables 2.1 Coordinates of Repeatability Measurement . . . . . . . . . . . . . . . . 3.1 Least Square Circle Fitting Properties for different Point Groups . . . . 3.2 Values of the stiffness of joint 5 for different configurations . . . . . . . 3.3 Joint Stiffness Values used in Elastostatic Model (Defined in N.m/rad) 3.4 Details of the designed tests for elastostatic model verification (a) . . . 3.5 Details of the designed tests for elastostatic model verification (b) . . . 3.6 Values of D 1 and D 2 for each experimental test path . . . . . . . . . . . 3.7 Joint Stiffness Values used in Elastostatic Model (Defined in N.m/rad) 3.8 Error Properties of Joint 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Error Properties of Joint 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Error Properties of Joint 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Different characteristics of the examined Honeycomb Materials . . . . 4.2 Material and Geometrical properties of Knife Structure . . . . . . . . . 5.1 different workpiece placement coordinates . . . . . . . . . . . . . . . . xvii

  FIGURE 1.1: (A) A Conventional Gantry Structured CNC machine (B)Stäubli TX200 Robot Arm

  FIGURE 1.2: Thesis Positioning and Configuration

Chapter 2 .

 2 Error Sources In Robotic-Based Machining joints are armed with external encoders implemented in control loop of robot joints (Figure 2.2).

FIGURE 2 . 2 :

 22 FIGURE 2.2: Configuration of a robot joint armed with external (secondary) encoder[START_REF] Möller | Machining of large scaled CFRP-Parts with mobile CNC-based robotic system in aerospace industry[END_REF]

FIGURE 2 . 3 :

 23 FIGURE 2.3: Robot configuration and results for compliance behaviour of Stäubli TX200 robot arm

FIGURE 2 . 4 :

 24 FIGURE 2.4: Cyclic TCP position error along a straight line for loaded and unloaded conditions

FIGURE 2 . 5 :

 25 FIGURE 2.5: Location of the chosen points for temperature measurements

  FIGURE 2.6: (A) Laser tracker used for position measurements. (B)Thermometer utilized for temperature measurements

FIGURE 2 . 8 :

 28 FIGURE 2.8: Location of the chosen spots for repeatability measurements with respect to robot coordinate system

FIGURE 2 .

 2 FIGURE 2.10: Contribution of each non geometrical error source to the robot pose misplacement

FIGURE 3

 3 FIGURE 3.1: Stäubli VS DH frames of Stäubli TX200 robot arm

FIGURE 3 . 2 :

 32 FIGURE 3.2: Computation of Self gravity torque imposed on joint 2

FIGURE 3

 3 FIGURE 3.3: schematic of balancing system attached to joint 2 of Stäubli TX200 robot arm

FIGURE 3 . 5 :

 35 FIGURE 3.5: Effect of Balancing system attached to joint 2 on the imposed torque

FIGURE 3 . 6 :

 36 FIGURE 3.6: Imposed torque on joints due to external load

FIGURE 3 . 7 :

 37 FIGURE 3.7: Illustration of measurement details of axis 1 stiffness identification using Faro arm

FIGURE 3 . 8 :

 38 FIGURE 3.8: Force imposition method for stiffness identification of joint 1

FIGURE 3 . 9 :

 39 FIGURE 3.9: Circle construction for joint 1 stiffness identification

FIGURE 3 .

 3 FIGURE 3.10: Measured and Fitted data for Stiffness identification of axis 1

FIGURE 3 .

 3 FIGURE 3.11: Configuration of the robot and force imposition system for joint 2 stiffness identification process

FIGURE 3 .

 3 FIGURE 3.12: Experimental data and fitted model to estimate stiffness of joint 2

FIGURE 3 .

 3 FIGURE 3.13: Stiffness configuration of joint 2

FIGURE 3 .

 3 FIGURE 3.16: Experimental data and fitted model to estimate stiffness of joint 3 (Left) and 2 (right)

FIGURE 3 .

 3 FIGURE 3.18: Measuring Process of Joint 5 Stiffness ID using Faro Arm

FIGURE 3 .

 3 FIGURE 3.20: experimental vs fitted data for joint 5 stiffness identification in different configurations

  1. Using generated elastostatic model with the experimental stiffness values 2. Using generated elastostatic model with the Stäubli set of stiffness values 3. Using load compensation feature of robot controller which uses the compensation method of Stäubli and Stäubli set of stiffness values

FIGURE 3 .

 3 FIGURE 3.21: Illustration of Desired, Actual and Compensated path for compensation algorithm

FIGURE 3 .

 3 FIGURE 3.22: Measurement Setup for Elastostatic Model verificationTABLE 3.6: Values of D 1 and D 2 for each experimental test path Test No. D 1 (mm) D 2 (mm) 1 1.4 0.5 2 2.4 1.2 3 0.8 0.3 4 0.2 0.1 5 0.9 0.3 6 0.9 0.5 7 0 0.3 8 0.5 0.3 9 0 0.3 10 1.5 0.7

FIGURE 3 .

 3 FIGURE 3.23: TCP path for different test conditions of Elastostatic Model verification

FIGURE 3 .

 3 FIGURE 3.25: Experimental and fitted values for test 1 and 2

FIGURE 3 .

 3 FIGURE 3.26: Experimental error values E for the entire rectangle on the left. Torques resulted from external force applied on different joint of robot arm on the right

FIGURE 3 .

 3 FIGURE 3.27: Contour of estimation error value Es for a range of K 2 and K 3 (Dimension of the ES is (mm))

  FIGURE 3.28: Experimental Vs simulated positioning error values of TCP for the rectangle experiment using different joint stiffness values

FIGURE 3 .

 3 FIGURE 3.29: Effect of compensation of Nonlinear behavior of transmission in a straight line in Olabi, 2011.

  FIGURE 3.30: Proposed Trajectory for Gear Transmission identification in Olabi et al., 2012

  FIGURE 3.31: Diagram of the Proposed methodology for gear error Identification

FIGURE 3 .

 3 FIGURE 3.32: Experimental setup and Measurement device for identification of gear transmission errors executed in LISPEN laboratory

  3.3. Cyclic Joint Transmission System Errors55Above mentioned values were determined due to the constraint of movement in the environment where robot was located and also reachability of the reflector target by the laser beam.

Figure 3 .

 3 Figure 3.33 depicts the Cartesian path resulting from the mentioned joint variables of the robot TCP defined in robot's base coordinate system.

FIGURE 3 .

 3 FIGURE 3.33: TCP path of proposed robot movement for gear error identification

FIGURE 3 .

 3 FIGURE 3.35: Periodic Gear Errors(∆θ jn ) (in Blue) and Fitted Sinus Functions (in Red)

  FIGURE 3.36: Gear Error Compensation Procedure

Figure 4 .

 4 Figure 4.1 presents different elements and properties affecting compliance errors in a robot-involved machining process with their corresponding section attributed in this chapter.

FIGURE 4 .

 4 FIGURE 4.2: Different elements combined in the ultrasonic cutting process in Le Creneau Industriel

Figure 4 .

 4 3 shows different cutting tools utilized by Le Creneau company.

FIGURE 4 .

 4 FIGURE 4.3: Different Cutting Tools Offered by Le Creneau. Ultrasonic Disks are shown at the top and Ultrasonic Knives are at the bottom of the figure

FIGURE 4 . 4 :

 44 FIGURE 4.4: Typical Honeycomb Panel End Design (Soovere, 1986)

FIGURE 4 . 5 :

 45 FIGURE 4.5: Overall configuration of a Chamfer

FIGURE 4 . 6 :

 46 FIGURE 4.6: General schematic of chamfering operation of Honeycomb work-piece

Figure 4 .

 4 Figure 4.7 is the final state of an ultrasonic cutting process in Le Creneau. The machining process of this work-piece combined of four chamfering operation for the four sides. The meteorology observations using Faro arm device reveal an error of about 4 [deg] for a desired 30 [deg] angle of chamfering.

FIGURE 4 .

 4 FIGURE 4.7: Example of a workpiece after four chamfering process on each side

Figure 4 .

 4 8 depicts a typical resulted workpiece.

FIGURE 4 .

 4 FIGURE 4.8: A resulted workpiece from the test session

FIGURE 4 .

 4 FIGURE 4.9: General schematic of chamfering operation of Honeycomb work-piece

  FIGURE 4.10: General schematic of chamfering operation of Honeycomb work-piece

FIGURE 4 .

 4 FIGURE 4.11: Penetration of knife in the protective moquette

kg m 3 FIGURE 4 .

 34 FIGURE 4.13: Coordinate system of the Kistler table on the left [Note that the XY plane coincides the surface of the Kistler ]. On the right: An illustration of a work-piece mounted on the Kistler

FIGURE 4 .

 4 FIGURE 4.14: Hard and Soft honeycomb materials used in the tests (on the left and right respectively )

  Chamfering angle (β) varied by the increments of 10 • (deg) from 80 • (deg) to 40 • (deg). The feed rate and direction of the cutting were 2( m min ) and positive Y direction in the sensor frame shown in Figure 4.13 and the lead angle was always set to 20 [deg].

Figure 4 .

 4 15 illustrates the configuration of the honeycomb workpiece placement and cutting knife orientation before execution of the cutting process.

FIGURE 4 .

 4 FIGURE 4.15: Chamfering test configuration and orientation of the end-effector (Cutting Knife)

FIGURE 4 .

 4 FIGURE 4.16: Measuring the entrance angle using a Goniometer

FIGURE 4 .

 4 FIGURE 4.17: Installation of the Workpiece for measuring the chamfering angle in Steady State zone

FIGURE 4 .

 4 FIGURE 4.18: Chamfering Angle error in both Soft and Hard Honeycomb materials for different desired chamfering angle values

FIGURE 4 .

 4 FIGURE 4.19: Displacement Errors as function of chamfering Angle of hard honeycomb material

FIGURE 4 .

 4 FIGURE 4.20: Ultrasonic Cutting Process Forces for a Chamfering operation for β = 50[deg] in hard honeycomb material

Figure 4

 4 Figure 4.20 reveals several notable characteristics of the honeycomb ultrasonic operation forces using knife tool. There are two zones in which all the components of machining forces are null. These zones correspond to the time before the first knife-material contact and after complete exit of the knife. Knife entrance zone depicted in the figure, corresponds to the time window between the first contact of the knife with honeycomb work-piece and the instance when knife has completely penetrated the material. In this zone, the active knife length increases by time and so do

  256.5 sin(β) + 63.1 cos(β) -314.9(N) f y (β) = -314.8 sin(β) -70.4 cos(β) + 348.8(N) f z (β) = 176.0 sin(β) -10.8 cos(β) -167.4

Figure 4 .

 4 Figure 4.21 illustrates different components of the modeled machining forces along with the experimental force values defined in Kistler coordinate system.

FIGURE 4 .

 4 FIGURE 4.21: Machining Forces in Hard and Soft Honeycomb chamfering process defined in Kistler sensor coordinate system

FIGURE 4 .

 4 FIGURE 4.22: Knife frame (in red) with respect to the Kistler frame (in black)

FIGURE 4 .

 4 FIGURE 4.23: Machining Forces in Hard and Soft Honeycomb chamfering process defined in Knife coordinate system

FIGURE 4 .

 4 FIGURE 4.24: Placement of the Nida Work-piece regarding to the robot

FIGURE 4 .

 4 FIGURE 4.25: Machining head assembly for the bench mark experiment

FIGURE 4 .

 4 FIGURE 4.26: Dimensions of the cutting knife placement and orientation with respect to the robot flange frame

  Following

4. 10 .

 10 Simulation of Robot Behaviour in honeycomb Machining Process 85 comes the transformation matrix between TCP and 6th frame:

FIGURE 4 .

 4 FIGURE 4.27: Evolution of CAM generated Euler angle components for TCP frame orientation during the chamfering process

Figure 4 .

 4 Figure 4.28 indicates that due to the data generated by the software, the term γθ s remains constant during the operation. This figure also shows the relative configuration of the TCP and Knife frames.

Figure 4 .

 4 Figure 4.28 depicts that during this operation spindle rotates continuously to keep the orientation of the knife on a desired constant configuration.

FIGURE 4 .

 4 FIGURE 4.28: Knife and TCP frame rotation regarding Spindle rotation angle

FIGURE 4 .

 4 FIGURE 4.29: Joint Variables for the bench mark test

FIGURE 4 .

 4 FIGURE 4.30: Cutting Forces of honeycomb material on each edge in presence of robot base frame

Figure 4 .FIGURE 4 .

 44 Figure 4.31 shows the results obtained from the simulations based on the three defined conditions for TCP displacement and deviation.

  8[deg] to describe an error of 4 [deg] in the chamfering process of the mentioned honeycomb workpiece. However, the angular deviations of the robot joints resulted from the simulations for the extreme condition (C), indicate drastically lower values (Figure 4.32).

FIGURE 4 .

 4 FIGURE 4.32: Joint Deviation in NO compensation mode on the left. Software generated displacement schema on the right

FIGURE 4 .

 4 FIGURE 4.33: Measurement Strategy for the effect of applied force on the deflection in the supporting table

FIGURE 4 .

 4 FIGURE 4.34: Dimensions of cutting knife involved in ultrasonic chamfering process

FIGURE 4 .

 4 FIGURE 4.35: Illustration of different knife structure directions

FIGURE 4 .

 4 FIGURE 4.36: Deformation of the beam model under a concentrated force in vertical direction

FIGURE 4 .

 4 FIGURE 4.37: Cross Section of the actual and modeled Knife geometry for vertical loading condition

FIGURE 4 .

 4 FIGURE 4.38: Moment imposed by a concentrated vertical force

FIGURE 4 .

 4 FIGURE 4.39: Deformation of the beam model under a concentrated force in Lateral direction

FIGURE 4 .

 4 FIGURE 4.40: Cross section of the knife for lateral bending deformation

FIGURE 4 .

 4 FIGURE 4.41: Deformation of the beam model under a concentrated force in Longitudinal direction

FIGURE 4 .

 4 FIGURE 4.42: Axial stress distribution of an element along the longitudinal axis

FIGURE 4 .

 4 FIGURE 4.43: Experimental setup along with loading and measurement schematic

FIGURE 4 .

 4 FIGURE 4.44: Experimentally driven with model predicted deformation of the knife structure under several vertical concentrated force condition

FIGURE 4 .

 4 FIGURE 4.45: Virtual joint model of the knife structure for vertical compliance effect

FIGURE 4 .

 4 FIGURE 4.46: Developed 7 axes robot structure for robot-tool compliance analysis

Figure 4 .

 4 Figure 4.47 depicts the loading condition for different n values.

FIGURE 4 .

 4 FIGURE 4.47: Different Load distribution conditions

FIGURE 4 .

 4 FIGURE 4.48: Bending moment in an arbitrary cross section in the right part

FIGURE 4 .

 4 FIGURE 4.49: Value of Error function (Er(n)) for different Loading Order (n) values

Figure

  Figure 4.49 illustrates the error of the estimating model of the compliance behaviour of knife structure for a wide range. This range covers the complete uniformly distributed load (n = 0) to the ranges where the force can be practically assumed as a concentrated force (n = 100). The behaviour of the graphs for both criteria follows the same trend. By marching along the n axis in positive direction, the estimated error decrease to a certain point and then increases. This clearly indicates that the assumption of a distributed load instead of a concentrated force is valid.

FIGURE 4 .

 4 FIGURE 4.50: Deviation and displacement error estimation of the knife structure using different loading order parameter (n)

Figure 4 .

 4 51 illustrates the chart of the above explained procedure.er < Vertical Force Model(Equation 4.5)

FIGURE 4 .

 4 FIGURE 4.51: Offline Chamfering Angle Compensation Logic

Figure 4 .

 4 52 depicts the values of β m as a function of desired chamfering angle β.

FIGURE 4 .

 4 FIGURE 4.52: Modified chamfering angle values (β m ) for a range of desired chamfering angle (β)

FIGURE 4 .Figure 4 .

 44 FIGURE 4.53: Uncompensated and compesated chamfering angle error values

FIGURE 4 .

 4 FIGURE 4.54: Vertical Machining Forces and Chamfering angle errors for the three experimental sessions

  FIGURE 5.1: Schematic of workcell for trimming process

FIGURE 5 . 3 :

 53 FIGURE 5.3: Illustration of Machining Forces

FIGURE 5 . 4 :

 54 FIGURE 5.4: Illustration of desired and actual path for Cost Function definition

FIGURE 5 . 5 :

 55 FIGURE 5.5: Examples of two actual paths with high error variance (A) and low error variance (B)

FIGURE 5 . 6 :

 56 FIGURE 5.6: Representation of discretization parameter (∆d) used for path definition

FIGURE 5 . 7 :

 57 FIGURE 5.7: Effect of discretization parameter (∆d) on the cost function evaluation error

FIGURE 5

 5 FIGURE 5.8: Defining variable r for shoulder singularity avoidance

  r

FIGURE 5 .

 5 FIGURE 5.10: Comparison of computation time in a jacobian based method and proposed method for singularity avoidance

  of hydropower plants Yildiz and Vrugt, 2019 and shape configuration in wave energy converters Esmaeilzadeh and Alam, 2019 to applications such as scheduling problems Sun et al., 2019, vibration control of vehicles He et al., 2019 and wind farm layout arrangement Baker et al., 2019 are among highly diverse domains of engineering prospering from contributions of such methods.

FIGURE 5 .

 5 FIGURE 5.11: Illustration of conventional and full-sized grid structure for free pattern search optimization method

FIGURE 5 .

 5 FIGURE 5.13: Evaluation of Cost function for a plane in workspace

FIGURE 5 .

 5 FIGURE 5.14: Effect of Frame rotation in evaluation of cost function

Figure 5 .

 5 Figure 5.16 shows an exaggerated displacement of the machining tool due to cyclic transmission error along the circumference of the trimming trajectory.

  FIGURE A.2: Positionnement et configuration de la thèse

  Figure A.3.

FIGURE A. 3 :

 3 FIGURE A.3: Contribution de chaque source d'erreur non géométrique à l'erreur de positionnement du robot

FIGURE A. 4 :

 4 FIGURE A.4: Calcul du couple d'auto-gravité imposé à l'articulation 2

FIGURE A. 5 :

 5 FIGURE A.5: schematic of balancing system attached to joint 2 of Stäubli TX200 robot arm

FIGURE A. 7 :

 7 FIGURE A.7: Couple de serrage imposé aux articulations sous l'effet d'une charge extérieure

FIGURE A. 8 :

 8 FIGURE A.8: Chemin TCP du mouvement du robot proposé pour l'identification d'erreur d'engrenage

FIGURE A. 9 :

 9 FIGURE A.9: Résultats expérimentaux pour les trajectoires compensées (en rouge) et non compensées (en bleu) pour une ligne droite

Figure

  Figure A.10 est l'état final d'un procédé de découpe par ultrasons au Creneau. Le processus d'usinage de cette pièce combine quatre opérations de chanfreinage sur les quatre faces. Les observations météorologiques utilisant le bras Faro révèlent une erreur d'environ 4[deg] pour un angle de chanfreinage souhaité de 30[deg].

FIGURE A. 10 :

 10 FIGURE A.10: Exemple d'une pièce après quatre chanfreins de chaque côté

FIGURE A. 12 :

 12 FIGURE A.12: Montage de la pièce à usiner pour la mesure de l'angle de chanfreinage dans la zone en régime permanent

  FIGURE A.13: Erreur d'angle de chanfreinage dans les matériaux en nid d'abeilles souples et durs pour différentes valeurs d'angle de chanfreinage souhaitées

  FIGURE A.14: Cadre du couteau (en rouge) par rapport au cadre Kistler (en noir)

FIGURE A. 16 :

 16 FIGURE A.16: Le déplacement TCP (à gauche) et l'écart (à droite) ont résulté de l'analyse de simulation pour les trois conditions de charge décrites

FIGURE A. 17 :

 17 FIGURE A.17: Illustration des différentes directions de la structure des couteaux

  FIGURE A.18: Développement d'une structure robotisée 7 axes pour l'analyse de conformité des robots-outils

  FIGURE A.19: Différentes conditions de répartition de la charge

FIGURE A. 20 :

 20 FIGURE A.20: Estimation des déviations et des erreurs de déplacement de la structure de la lame à l'aide de différents paramètres d'ordre de chargement (n)

  FIGURE A.21: Valeurs d'erreur d'angle de chanfreinage non compensées et compensées

  FIGURE A.22: Schéma de la cellule de travail pour le processus d'ébavurage

FIGURE A. 23 :

 23 FIGURE A.23: Illustration des forces d'usinage

∆

  FIGURE A.25: Résultat final du placement de la pièce à usiner

Le chapitre 4

 4 traite d'un cas pratique d'usinage robotisé. Dans ce chapitre, une installation d'usinage composée d'un bras robotisé, d'une broche et d'un outil d'usinage conçu pour couper des matériaux en nid d'abeille à l'aide de la technologie ultrasonore a été adoptée comme étude de cas. La principale préoccupation dans ce processus d'usinage était la présence d'un niveau élevé d'erreurs géométriques indésirables observées dans les pièces résultantes. Les erreurs géométriques d'usinage ont été mesurées à l'aide d'une approche pratique présentée dans le chapitre.

SSM Structural Stiffness Matrix CFRP Carbon Fiber Reinforced Plastic Chapter 1

	IR	Industrial Robot
	TCP	Tool Center Point
	LT	Laser Tracker
	FK	Forward Kinematics
	IK	Inverse Kinematics
	FEA	Finite Element Analysis
	MSA Matrix Sturactural Analysis
	VJM	Virtual Joint Method

TABLE 3 .

 3

	1: Least Square Circle Fitting Properties for different Point
		Groups
	Point Group No. R (mm) Standard Deviation (mm)
	1	1424.0	0.024
	2	1333.6	0.023
	3	1194.2	0.040

TABLE 3 .

 3 2: Values of the stiffness of joint 5 for different configurations

	Conf No. θ 5 ( o deg) k 5 (N.m/rad)
	1	31.65	2.2651e5
	2	-40	2.5938e5
	3	-60	1.67788e5

  Consider a rectangle with following edge coordinates defined in the robot base frame:

	Point X(mm) Y(mm) Z(mm)
	P 1	900	-600	445
	P 2	1900	-600	445
	P 3	1900	600	445
	P 4	900	600	445

TABLE 3 .

 3

	Test Number	6	7	8	9	10
	Actual Load (N)	500	500	0	500	0
	Compensation Method	ESM	ESM	Stäubli Stäubli Stäubli
	Compensation Load (N)	500	500	0 (SG)	500	500
	Stiffness Values Source Experimental Stäubli Stäubli Stäubli Stäubli

4: Details of the designed tests for elastostatic model verification (a)

TABLE 3 .

 3

5: Details of the designed tests for elastostatic model verification (b)

TABLE 3

 3

		.8: Error Properties of Joint 1
	k	a k (rad)	b k	c k
	1 3.838e-05	66.8	2.844
	2	1.65e-05	89.07	-1.567
	3 1.314e-05 44.54	-3.046
	4 8.299e-06 200.4	1.474
	5 7.979e-06 111.3	2.508
	6 7.675e-06 22.27 -0.9643
	7 5.072e-06 133.6	-1.505

TABLE 3

 3

		.9: Error Properties of Joint 2
	k	a k (rad)	b k	c k
	1 3.002e-05	64.2	0.8516
	2 2.122e-05 33.92 -1.009
	3 1.493e-05 127.8	3.624
	4 1.186e-05 29.25	1.55
	5 1.316e-05 68.55	3.401
	6	9.42e-06	198.4	2.081
	7 7.494e-06 59.52 -1.898

TABLE 4 .

 4

	kg m 3 ) Manufacturer

1: Different characteristics of the examined Honeycomb Materials Given Name Aramid Polymer Type Density (

TABLE 4 .

 4 2: Material and Geometrical properties of Knife Structure

	Parameter	h	t	E	b 0
	Value	7 × 10 -2 (m)

2 × 10 -3 (m) 55 × 10 10 (Pa) 12.75 × 10 -3 (m)

TABLE 5 .

 5

		1: different workpiece placement coordinates
	Point Origin	x c (mm) y c (mm) z c (mm) θ[deg] Cost value (mm)
	Proposed Method	572.2	480.1	518.5	-49.6	0.187
	Intuitive Placement	1500	200	1200	0	1.697

  rapportent un résultat contrasté de leurs recherches sur un bras robotique exécutant une tâche d'usinage. De plus, les progrès continus dans les méthodes de fabrication offrent des éléments plus précis avec des tolérances géométriques plus faibles promettant des valeurs d'erreurs géométriques réduites.Une classification générique des erreurs non géométriques est apparemment impossible en raison du grand nombre de sources qui répondent à sa définition. De plus, ces sources ne se limitent pas aux aspects mécaniques d'un système de robot industriel (IR). D'autres sources provenant d'aspects électriques et de programmation participent également, à un certain niveau, aux erreurs non géométriques. L'effet de bande passante limitée mentionné dans[START_REF] Kim | Robotic Machining: A Review of Recent Progress[END_REF] est un exemple d'erreur non géométrique possédant une nature autre que mécanique.

(A) (B)

Linear Variable Differential Transformer

FIGURE 5.12: Optimization Flow chart

Une cellule robotisée d'usinage contient au moins trois éléments[START_REF] Olabi | Amélioration de la précision des robots industriels pour des applications d'usinage à grande vitesse[END_REF]: Bras de robot, outil d'usinage et pièce à usiner. La conception de la cellule robotisée dépend fortement du processus d'usinage ciblé. Certains procédés tels que le fraisage de blocs métalliques présentent des forces d'usinage élevées et il convient donc d'être particulièrement attentif à améliorer la précision de l'usinage grâce à un traitement de comportement de conformité du robot. Cependant, dans des applications de ce type, le comportement de conformité des outils d'usinage, par exemple les

Sim (A) Sim (B) Sim (C)

En ayant la rigidité torsionnelle équivalente de l'outil d'usinage, qui est maintenant considéré comme le 7ème maillon du robot, le comportement de conformité de la structure du robot-outil sera modélisé selon la même approche qu'un robot 6 axes ordinaire. Comme le comportement de conformité de la structure du robot

Dans ce chapitre, nous aborderons une procédure de compensation des erreurs induites par le comportement de conformité d'un processus d'usinage robotisé dans lequel l'effet de conformité de l'outil d'usinage ne peut être négligé par rapport à celui de la structure du robot. Bien que la documentation sur les procédés d'usinage pour lesquels la conformité de l'outil n'est pas un problème (comme le perçage et le fraisage) soit abondante, le but de ce chapitre est d'introduire un schéma pour contrer l'effet de conformité de l'outil d'usinage ainsi que de la structure du robot dans les procédés d'usinage combinant un outil conforme.

Le procédé de chanfreinage sur nid d'abeilles a été adopté comme procédé d'usinage d'étude de cas, car ce procédé d'usinage est livré avec un outil conforme et conforme aux objectifs du partenaire industriel de ce travail Le Creneau Industriel.

Chapter 5. Optimal Workpiece Placement In Robotic Cells model for our case study robot Stäubli TX200, we make use of the model generated for Stäubli TX90 described in section 3.3.4. The implementation of this error type in the optimization method consists of three main steps. For a given placement (x c , y c , z c ) and frame rotation θ of the workpiece, the desired joint variables are computed using the inverse geometrical model of the robot. Secondly the desired joint variables are fed into the error model to compute the corresponding transmission errors. Eventually, the error values are added to the desired joint variables to calculate for the actual placement of the TCP using robot's forward geometrical model. Figure 5.2 illustrates the described procedure.

Elastostatic Error Source

The compliance behaviour of serial robots can be modelled using different approaches.

As explained in [START_REF] Klimchik | Stiffness modeling for perfect and non-perfect parallel manipulators under internal and external loadings[END_REF], these approaches can be named and ordered by their complexity into the following three categories: Virtual Joint Method (VJM), Matrix Structural Analysis (MSA) and Finite Elements Analysis (FEA). It is evident that these methods differ in the accuracy of the robot compliance behaviour prediction based on the level of detail involved in the method, but on the other hand the computation and effort taken for high level of details may not be cost-effective for certain ranges of machining processes. Due to this fact, even recent research works such as Theissen, Laspas, and Archenti, 2019, mainly adopt the VJM method for elastostatic modeling of robot structure.

The VJM method considers robot links as rigid structures undergoing no deformation due to the external load and replaces robot joints with linear rotational springs. To calculate the displacement of the robot end effector (∆X) to the external couple (F), the following relationship is given by [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF]:

Where:

In equation 5.3, matrix J is the jacobian matrix of robot manipulator relating the Cartesian velocity screw of the end-effector defined in the base frame to the angular velocities of the joint space variables. k θ is the diagonal matrix involving the joint stiffness values as follows:

Appendix A

Résumé en français

A.1 Introduction

Les progrès réalisés dans divers domaines de l'ingénierie tels que l'électronique et la mécanique ont élargi les horizons des applications de l'ingénierie. Les robots sont l'un des exemples parfaits de dispositifs rendus possibles par des réalisations simultanées sur des sujets autrefois non relatifs. Avec les progrès technologiques et scientifiques, le mouvement du cinéma et de la littérature de fiction à la fin des années 70 a pris le pouvoir comme deuxième pilier pour imaginer l'avenir des vies humaines affectées par les robots.

En bref, bien que l'avenir dépeint dans lequel les humains se distinguent encore de leurs robots émotionnels et structurés de façon identique, les robots occupent progressivement les fonctions et les positions d'experts humains dans de nombreuses applications.

Les bras de robots industriels (IR) sont un dérivé de robots conçus et fortement impliqués dans différentes applications industrielles. Cependant, les applications dans lesquelles les robots industriels étaient impliqués, ont été soumises à des changements. Au départ, les RI n'étaient utilisés que pour des applications telles que le pick and place et la peinture. Ces applications sont connues pour leur exigence de répétabilité au lieu de précision. Auparavant, les exigences des applications industrielles étaient relativement immuables dans le temps. Les aspects géométriques et fondamentaux des pièces ont rarement fait l'objet de mises à jour. Par conséquent, la versatilité et l'adaptabilité des lignes de production n'ont pas été considérées comme une propriété vitale. D'autre part, aujourd'hui, grâce à des progrès significatifs dans des domaines tels que les logiciels de simulation numérique, la capacité de calcul des ordinateurs et la science des matériaux, combinés à l'évolution rapide du marché dans des secteurs tels que l'aérospatiale et le transport, l'industrie manufacturière est en pleine révolution pour proposer des concepts et des configurations dont la principale caractéristique est la capacité à évoluer.

A la lumière de cette révolution industrielle, une attention particulière a été portée aux RI en raison de caractères notables en ligne avec les objectifs des besoins de fabrication récents. Ils disposent d'un espace de travail relativement grand, facile à programmer avec un poids plus léger par rapport aux machines d'usinage conventionnelles et moins coûteux. Les imperfections observées dans les pièces résultantes d'un processus d'usinage robotisé proviennent de différentes sources. Du côté du robot, ces erreurs, comme mentionné dans les chapitres précédents, proviennent de sources géométriques et non géométriques. Mais comme indiqué précédemment, le comportement de conformité du robot est la principale cause de mauvais positionnement de l'outil pour un bras de robot chargé pour un IR lourd. De plus, dans certaines applications d'usinage, l'effet de la conformité des outils d'usinage ne peut être négligé dans les erreurs de fabrication. d'usinage supplémentaire et 3. Les forces d'usinage. l'écart angulaire dans les articulations du robot (∆Θ) imposé par un ensemble de valeurs de couple appliquées aux articulations du robot (τ t ), peut être calculé en utilisant l'expression suivante :

le déplacement TCP peut être obtenu en comparant le placement TCP résultant des variables articulaires réelles avec celui résultant des variables articulaires désirées en utilisant le modèle géométrique avant du bras du robot.

La figure A.16 montre les résultats obtenus à partir des simulations basées sur les trois conditions définies pour le déplacement et l'écart TCP.

A.5.2 Définition de la fonction de coût

Dans ce travail, nous proposons une fonction de coût conçue pour minimiser les erreurs d'usinage side-directionnelles indépendantes de la source d'erreur concentrée. La figure A.24 montre un schéma d'un chemin souhaité et d'un chemin réel résultant d'erreurs liées au robot. A la recherche d'une fonction de coût appropriée imposant une contrainte plus forte sur les erreurs cycliques ainsi que sur les erreurs de comportement conforme, nous proposons la variance des valeurs d'erreur (var{e i }) qui, grâce à sa définition, élimine l'effet de la valeur moyenne des erreurs :

Σ(e ie) 2 N (A.12)

A.5.3 Implémentation des contraintes

Parmi toutes les contraintes et les conditions limites en robotique, le respect des limites articulaires est de la plus haute importance. Les limites articulaires sont imposées en raison des limites physiques de la structure du robot et leur violation peut causer des dommages irréversibles aux appareils relativement coûteux tels que les bras robotiques industriels. Dans ce travail, nous avons adopté la fonction Tsai pour imposer la contrainte de limite commune sur le problème d'optimisation.

La prochaine contrainte importante dans toute application robotique est l'évitement de la singularité. Grâce à la conception du bras robotique Stäubli TX200 à cloison de poignet, une configuration singulière du robot n'est atteinte que dans les trois conditions suivantes [START_REF] Doan | Optimal robot placement with consideration of redundancy problem for wrist-partitioned 6R articulated robots[END_REF][START_REF] Hayes | Singular configurations of wrist-partitioned 6R serial robots: A geometric perspective for users[END_REF]: Singularité du poignet, Singularité du coude et Singularité de l'épaule.

Pour éviter les configurations singulières dans cette étape, nous proposons également une approche de pénalisation de la condition de limite commune. Dans cette approche, le calcul du jacobien est évité et donc une énorme charge de calcul est allégée. Un coefficient de singularité (S) est proposé pour pénaliser les résultats de la fonction de coût.

Appendix A. Résumé en français modèle a été vérifié par des mesures expérimentales. Les résultats ont montré que l'effet de conformité de la structure du couteau provient principalement de la direction verticale. L'ensemble du système robot-outil a ensuite été modélisé par un robot à 7 axes prenant la direction verticale du couteau comme une articulation virtuelle. Le comportement de ce nouveau système a été analysé et les résultats ont montré que l'effet de la compliance du couteau était plus de 8 fois supérieur à celui de la structure du robot.

Afin de réduire les erreurs d'usinage, les différentes conditions de charge appliquées à l'outil d'usinage ont été examinées. Des simulations ont montré qu'une charge non uniformément répartie le long de l'axe longitudinal du couteau peut le mieux décrire les erreurs angulaires de chanfreinage mesurées pendant les séances expérimentales. En utilisant le modèle de conformité de la structure de la lame et le modèle des forces d'usinage dans ce processus, une méthode de compensation hors ligne basée sur un algorithme itératif a été proposée. Dans cette méthode de compensation, l'algorithme a été conçu pour fournir un angle de chanfreinage secondaire correspondant à un certain angle de chanfreinage souhaité. L'expérience a prouvé l'efficacité de la méthode de compensation qui a permis de réduire de 40% à 90 % les erreurs géométriques. Cependant, on a observé que le processus d'usinage par chanfreinage possédait un caractère non répétitif en termes de forces d'usinage.

Pour surmonter les inconvénients rencontrés lors des étapes de post-processus dans les approches de compensation hors ligne, au chapitre 5 une méthode d'optimisation a été proposée pour minimiser les erreurs non géométriques d'un IR en remplaçant la pièce dans une configuration proche-optimale. La tâche d'usinage ciblée était l'ébavurage d'une pièce circulaire en plastique renforcé de fibres de carbone (CFRP). Les sources d'erreur soumises étaient le comportement de conformité du robot et les erreurs du système de transmission commun.

Une fonction de coût a été proposée sur la base de la variance des composantes des erreurs de déplacement normales à la direction de l'avance. Parmi les contributions de ce chapitre, l'approche de mise en oeuvre des contraintes est remarquable. Des limites articulaires et l'évitement de la singularité ont été mis en oeuvre dans la fonction à l'aide d'une stratégie pénalisante. Une ligne droite a été ajustée aux ensembles de valeurs des articulations 2 et 3 pour lesquelles le robot fait face à des configurations de singularité d'épaulement. Cette ligne a été implémentée dans les critères d'évitement de singularité pour éviter les calculs jacobiens. La méthode d'implémentation des contraintes a considérablement augmenté l'efficacité informatique du processus d'optimisation.

Un algorithme de recherche de modèle modifié a été utilisé pour exécuter le processus d'optimisation. Une grille pleine grandeur a été adoptée dans l'algorithme pour faciliter la décision quant à la direction de progression pour un point initial proposé. Les résultats de la simulation ont montré qu'une grande partie des erreurs non géométriques provenant du bras robotique peuvent être éliminées en replaçant la pièce à usiner dans un endroit proche et optimal.

En obtenant les conclusions mentionnées ci-dessus tout au long de cette thèse, certaines propositions pour la poursuite de la recherche dans ce domaine méritent d'être mentionnées. Pour augmenter la précision du modèle élastostatique du robot, The involvement of industrial robotic arms in different manufacturing applications is going through an ever-changing era. These devices which were once only used in applications based on their acceptable repeatability are now taking place of dexterous human operators or rigid machining devices in high precision manufacturing tasks. However, the inherited shortcomings of these structures require further investigations to ensure acceptable performance. In this thesis, efforts were dedicated firstly to observe the contribution of different error sources in robotic involved manufacturing. The main contributors to the robot inaccuracy were investigated for identification and compensation processes. An elastostatic model of an industrial robot is generated. A novel method in identifying joint transmission errors is addressed capable of massively reducing the identification time. To enhance the workpiece quality in robotic-based honeycomb ultrasonic machining process, a new approach to integrate tool compliance behaviour with the one of robot structure is proposed. A detailed procedure of developing the machining force model and geometrical error measurement is also presented for this machining process. To reduce the effect of non-geometrical error sources, an optimization process is executed for optimal placement of workpiece in robot workspace.

Keywords : Robotic Manufacturing, Non Geometrical Errors, Ultrasonic cutting, Honeycomb Material, Workcell Optimization, singularity avoidance