Miguel Ortiz

David Bétaille
email: david.betaille@ifsttar.fr

François Peyret

Ola Martin Lykkja

Svend-Peder Oseth

Assessment of End-To-End performances of a GNSSbased Road User Charging System Based on Ifsttar-GEOLOC and Q-Free cooperation

Keywords: GNSS positioning, intelligent transportation systems, road user charging, error modelling I

The paper deals with the issues that have to be faced when designing a national Road User Charging (RUC) in Norway based upon GNSS and with an innovative methodology especially developed to address these issues proposed by CEN-CENELEC TC5/WG1 and supported by the SaPPART COST Action. The paper presents a case study applying the methodology to a specific RUC algorithm developed by the Norwegian company Q-Free. A dedicated GNSS position error model has been produced by Ifsttar GEOLOC for this purpose.

INTRODUCTION

A. Context

Global Navigation Satellite Systems (GNSS) have a very high potential in the development of Intelligent Transport Systems (ITS), Personal Mobility and associated services. This has been widely demonstrated through the use of GPS in supporting the provision of ITS services such as personal navigation, fleet management, cooperative traffic monitoring and more recently: Road User Charging (RUC), Pay-As-You-Drive insurance, emergency call (eCall), tracking and tracing of dangerous good, Advanced Driver Assistance Systems (ADAS), etc. [START_REF]GNSS Market Report, issue 4[END_REF].

Given the principle of GNSS positioning, performance is highly influenced by the conditions of the operational environment. Therefore, GNSS integrators and users are facing two major challenges: the difficult problem of estimating the expected performance of the service when using GNSS, and the lack of standards and certification references on positioning performance, which are necessary to guide their choices [START_REF]Satellite Positioning Performance Assessment for Road Transport[END_REF] [START_REF]SaPPART White paper: Better use of Global Navigation Satellite Systems for safer and greener transport[END_REF].

B. Road User Charging with GNSS

As alternative to conventional physical-gantries-based RUC systems, GNSS-based RUC systems offer numerous advantages, the most important one being their flexibility. In Europe, these systems are being adopted by more and more countries. Slovakia, Germany, Switzerland and Hungary have successfully implemented GNSS-based tolling, Belgium and Russia have launched similar projects implementing GNSSbased schemes and France, Finland, Bulgaria, Denmark, The Netherlands and Lithuania have all declared their interest in GNSS-based schemes [START_REF]GNSS Market Report, issue 4[END_REF].

But, RUC is very demanding in terms of false detections which are prohibited to avoid the risk of rejection of the system by the end-user. For instance, for the French system, the acceptable false detection rate of the system had been fixed at 10-6 with a confidence level of 5%, which could be reached for example if there were no missed detection at all during 3.106 consecutive charging events all over the network. To be capable of such performance level, the detection of the charging event (the vehicle crossing the charging point) has to be extremely reliable, which can be quite tricky in difficult environments.

The detection performance of the system at the charging point depends on the performance of both the GNSS-based positioning terminal (GBPT) and the RUC detection algorithm in the given charging points environment. In challenging environments, such as urban locations, where poor GNSS performances are usual, a very smart detection algorithm can overcome these poor GNSS performances and ensure a high performance for the whole system. This means that, performance assessment of such system needs tests applied to the whole system, with a perfect respect of the environmental conditions [START_REF]SaPPART White paper: Better use of Global Navigation Satellite Systems for safer and greener transport[END_REF].

Figure 1 illustrates the 2 main components which impact the end-to-end (E2E) service of the system.

C. Assessment of Road User Charging system performances

The execution of realistic test campaigns for assessing such performance as a false detection rate of 10-6 faces a difficult issue: how to carry out millions of tests in real or realistic environments ? Obviously, this is not possible to be done in real full-scale conditions, and laboratory-only conditions would not allow the respect of the realistic environment, so, a combined approach has to be followed.

CEN-CENELEC TC5 WG1, which has addressed this issue in his 1st standard [START_REF]Space -Use of GNSSbased positioning for road Intelligent Transport Systems (ITS) -Part 1: Definitions and system engineering procedures for the establishment and assessment of performances[END_REF], proposes a general method, applicable to any ITS system based on GNSS, called "Sensitivity analysis". This method is based upon field tests in real conditions of the GBPT to identify a PVT error model that will be used to generate automatically a high number of synthetic degraded trajectories that will be processed by the application module. From the outputs of the application module, using the E2E performance metrics, will be assessed the E2E performances of the system, depending from performances of both the GBPT and the application itself.

This method offers the advantage to multiply field tests executed in real operational conditions, the only ones capable of capturing the reality of the physical phenomena, in order to run a high number of trials, necessary to assess performances which are generally expressed by very low probabilities.

The sensitivity of the system to the performance of the GBPT can be analysed by increasing artificially step by step the amplitude of the errors of the PVT error model until the moment when the target E2E performance is no longer fulfilled. To summarize it, the Sensitivity analysis consists in 6 main steps:

Safety margin analysis

More details can be found in the draft standard itself [START_REF]Space -Use of GNSSbased positioning for road Intelligent Transport Systems (ITS) -Part 1: Definitions and system engineering procedures for the establishment and assessment of performances[END_REF]. This methodology was already followed, under an early form, by the test laboratory which was in charge of assessing the false detection rate performance of the French RUC system [START_REF] Duchâteau | Challenges in GNSS Road User Charging and French Implementation Experience[END_REF]. The following sections will be devoted to the application of this methodology to the RUC algorithm proposed by Q-Free Company.

The three first steps have already been achieved and already reported in previous papers [START_REF] Peyret | How GNSS performance standardisation can support the deployment of critical its applications[END_REF]. The next section will recall the main principles used for the development of the error models and the following ones will describe in more details how the degraded trajectories have been generated and used for assessing the performances of the RUC algorithm.

II. ERROR MODELLING AND SIMULATION OF DEGRADED

TRAJECTORIES

The error model is based upon the decomposition of the horizontal position error vector into its two local polar components: radius (or modulus) R and angle φ relative to the local Frenet frame, tangent to trajectory of the vehicle. Figure 2 shows these two local polar components. The 2 components of the error data acquired from field tests that have been chosen as seeds of the synthetic degraded trajectories are respectively modelled by: -For the radius, a filtered piecewise constant signal where a new step appears with the probability p (a Bernoulli trial at every sample); in case of a step, a new amplitude R for the radius is drawn from a Laplace distribution and a Cauchy additive noise, which is accumulated over each step, is added to it:

R[n+1] = { + C r [n] (1)
where L r [n] is a realization of a truncated Laplace distribution (positive values only) giving the height of the new step, and C r [n] is a realization of a truncated Cauchy distribution.

The accumulation of the Cauchy noise does not operate over the whole signal, but it is accumulated only during one step. It is reset when there is a new draw from the Laplace distribution.

In a final step, the radius signal is filtered using an autoregressive filter which relates the expected error signal with the simulated signal.

-

For the angle, a random walk in which the additional noise is a realisation of a Cauchy distribution:

R[n] with probability (1-p) L r [n] with probability

p φ[n+1] = φ[n]+S φ [n] (2)
where S φ [n] is a realization of a truncated Cauchy distribution. The dispersion of the angle distribution (parameter γ) will be modulated with the speed of the vehicle, so that when stopping, the error may vary in amplitude but not in direction (avoiding that way a possible large rotating error).

The final angle error signal is also low-passed filtered.

The identification of the parameters of the Laplace distribution of the radius is realized with [START_REF] Hutter | Exact Bayesian regression of piecewise constant functions[END_REF]. This method gives us the parameters of the distribution for the step amplitudes and also the instants when these steps occur, from which we can estimate p. An ARMAX identification is used to identify the linear filter for the radius [START_REF] Lennart | System identification: theory for the user[END_REF]. The identification of the parameters of the Cauchy distribution of the angle is realized by fitting the experimental PDF [START_REF] Bloch | A note on the estimation of the location parameter of the Cauchy distribution[END_REF].

Once the models of the radius and the angle components of the error have been correctly identified and tuned, the generation of the along-track and cross-track components of the error is straightforward:

{ (3)
where w a and w c are normalisation constants used for the situations where the error vector does not distribute equally in both directions [START_REF] Groves | Intelligent Urban Positioning using Multi-Constellation GNSS with 3D Mapping and NLOS Signal Detection[END_REF].

A. GNSS data collection

The GNSS position error modelling, as mentioned in item 3 of the Sensitivity analysis, shall rely on various full scale experiments reproducing operational scenarii. In the frame of this study, a test vehicle is equipped with several low-cost GNSS receivers from uBlox, ST, and SiRF, among which a low cost receiver uBlox LEA6-T is considered. The equipped vehicle is shown on Figure 3. The test devices were placed behind the windshield.

A reference trajectory is needed, not only for the error modelling step of the work, but also later on as a basis for cloning. The reference trajectory, supposed to be errorless, will actually be degraded N times, applying N different radius and angle errors, so that N clones of the initial uBlox trajectory will be able to be built in order to be fed into the RUC process. The synthetic errors are built by drawing random values in each of the probabilistic laws that are used in the radius and angle models.

On the roof one can see the antenna of the reference NovAtel CPOS receiver. This receiver is hybridized with an IMU and post-processed by the kinematic Precise Point Positioning (PPP) Terrapos software. The data set is composed of 28 trajectories from Frankfurt, Germany, each of them being approximately 5 km long (see Figure 4). They were performed at the average speed of 15 km/h. Despite it is without tunnels and covered sections, the driven route goes through deep urban canyons, as visible in Figure 5.

X a = w a R[n] cos (φ) X c = w c R[n] sin (φ)
large error peaks of 400 meters for a duration of up to 30 seconds (see Figure 6).

B. RUC virtual gantries and performance metrics

The RUC application is based on virtual gantries (VG) defined as geo-objects specified by latitude and longitude. They comprise a tolling point (pink circle on Figure 7), a central polyline, and cross-sections (ABCD) with left and right extensions of the road (dashed line), and outer tolerances to the left and right (solid line), in meters. The tolerances can be made very wide if there is no neighboring road. A virtual gantry is oriented. When a vehicle passes through the gantry, the intersection point between the GNSS trajectory and the crossing lines are identified and the distance from the center point to the intersection point is calculated. This distance is compared to the road width and the tolerance and gives a passage score, with a maximum value at 1 if the intersection coincide with the center point, and going to zero at the outer tolerance point. If the intersection is outside the tolerance limits, the passage of this crossing line is ignored. After the passage through the VG, the total score is calculated and if the total score is above a predefined limit and more than p out of q crossing lines are crossed and the crossing lines are crossed in correct order, then a VG passage is assessed. The score can be used as an indication of VG passage confidence.

The arrangement makes the passage detection independent of driving speed and GNSS sampling rate. The method also spreads the decision out on a large geographical area (several hundreds of meters, even kilometers in rural road scenario, with multiple crossing lines) thus reducing the effects of the correlation between individual GNSS samples. The VG can be applied to curved roads and also in (or close to) intersections. This makes the method also suitable for urban use. The tolerance and width information are set individually for each crossing line based on local geographical and road topology properties taking into account typical GNSS performance at that location.

A first set of virtual gantries have been defined manually all along the itinerary repeated in Frankfurt (in blue). These virtual gantries should be detected and they will make it possible a measurement of the missed detections of the process. Additional virtual gantries (in orange) were added in the vicinity of the first set. The aim of this second set is to create false detections and make it possible a measurement of those. Some of these additional virtual gantries belong to other streets, but some are also totally fake. They are located in such a way to trigger false detections. Thus the false detection rates seen in this article does not reflect at all those of a real system. ^^^ From Ola's txt in BLUE • blue virtual gantry: virtual gantry in the original QFREE database • orange virtual gantry: virtual gantry added specially for false detection purpose. •CCER Correct Charge Event Recognition = (NbCorrectDetections for VGi) / (RefNbCorrectDetection for VGi)

• FPCER False Positive Charge Event Recognition = (NbFalseDetections for VGi) / N For instance here is the result of the assessment of the RUC algorithm with 17,000 cloned trajectories, with a certain set of error modeling parameters: The speed is normalized by the average speed during the test campaign in Frankfurt (approximately 15 km/h) and it is used to modulate the parameter γ of the Cauchy law (i.e. the dispersion of the distribution).

One simulation result (for the radius) is shown in Figure 11 along with the true error radius of the best trajectory (25), using the default parameters. Another result, using the maximum upper bound of 400 m, is shown in Figure 12 along with the true error radius of the worst trajectory (21).

IV. ANALYSIS OF THE RESULTS OBTAINED

A. Fixing the number N of clones

The first parameter to set is the number N of simulated trajectories we have to produce in order to be confident that computed metrics are representative of the SA set. For instance, N=1 will not produce representative metrics: there is a good probability that CCR would be 100% and OVR 0%. By contrast, N= 10,000,000 will produce much more representative metrics ; but it will take too much time to run the whole process on a computer. For information, it takes about 1,5 h of computation for N=10,000 on a quite robust laptop. In order to estimate the minimum N value to be used, a run with N=30,000 has been conducted.

0 0.5 1 1.5 1.7 2 2.5 3
Number of simulated trajectories (N) The red line corresponds to the mean of the last third part (from 20,000 to 30,000); the red dash lines correspond to +10% and -10% of the first oscillation. As we can observe, from N=17,000 oscillations stay inside the 10% bounds. This value of N=17,000 will be used for the next phase of RUC assessment.

B. Sensitivity analysis

This sensitivity analysis (SA) has focused on varying the upper bound of the Laplace law which drives the radius error. The two following figures illustrate the random draws for 5 different sets of SA. Afterward, the notation SA(11-36-x) means SA set with all radius and angle parameters set to values of Table I and II, excepted for the maximum of radius Laplace parameter equal to x (cf. Table III). As our original reference trajectory on which cloned radius and angle errors had been added is a 1138 samples trajectory, the Laplace random draw (Figure 16) is not so close to the ideal plot (Figure 15). However, we could observe that some high errors should be draw with low occurrence. 85,000 trajectories have been injected in the RUC algorithm -17,000 for each SA -with the pre-defined set of virtual gantries (Figure 8): We can observe an exponential trend of the evolution of CCR and OCR over the Laplace upper bound (Figures 17 and18). In order to confirm this trend (Figures 19 and20), other runs have been realized with different sets of SA (from 11-36-22 to 11-36-400): Note 1: For CCR the 3 first points have been discarded for exponential fit. Note 2: OCR results are not representative of a realistic set of virtual gantries. They are only representative of the fake virtual gantries that had been added to the original ones. Note 3: CCR results are representative of a very dense network of virtual gantries. Here, each large error makes a missed detection. We could imagine that a low density network will be more resilient to large position error.

The exponential trends are confirmed. These asymptotic trends could be explained by the Laplace parameter we have chosen to make varying trajectories. The average of such a random draw with the upper bound varying is indeed nearby 40 m.

C. Conclusions

This study has allowed assessing a real Road User Charging (RUC) application. Cloned (or simulated) trajectories have been built from a model issued from a u-Blox-6T GPS receiver. It is based on Radius generation on one hand, and Angle generation on the other hand. Radius and Angle generation are driven by different random laws, and their parameters are the ones defined in Table 1. This kind of model fits quite well real trajectories with low errors, and has more difficulties to fit large error trajectories.

Assessment of the RUC algorithm could have been done by injecting trajectories increasingly degraded. Up to 18 sets of 17,000 trajectories have been injected. These parametric tests have been based on the increase of the maximum Laplace error. The following global metrics have been computed: Correct Charging Rate (CCR) and Over Charging Rate (OCR). Metrics relative to each virtual gantry have also been computed: this allows identifying a possible virtual gantry which has been badly designed. As a matter of fact, all false detections have been generated by virtual gantries that have been intentionally designed for that purpose (virtual gantries added to create false alarms). Furthermore, assessment of the RUC application has shown that CCR and OCR reach an asymptotic value when the upper bound of the Laplace law is increasing. In our study, these asymptotic values themselves make no general sense because the use case is here totally artificial and nonrepresentative of a real RUC system, but the existence of these asymptotes tends to prove that the couple {virtual gantries + RUC algorithm} presents minimum performances in front of a given receiver. These performances are actually guaranteedand in particular over-charging cannot exceed the maximum rate estimate -even if the experienced trajectories are all very bad ones.

Figure 1 :

 1 Figure 1: The 2 main components of a RUC system impacting the E2E service

1 .

 1 Definition of the operational scenario and of the test protocol 2. GNSS-based positioning terminal (GBPT) field tests execution 3. PVT error model identification 4. Generation of degraded simulated trajectories 5. E2E performances assessment

Figure 2 :

 2 Figure 2: The local Frenet frame and the two polar components of the position error vector

  <<< until here, in GREEN, is what I added to the draft >>> from here, in YELLOW, is what I took from the draft III. APPLICATION OF THE DEGRADED TRAJECTORIES TO ROAD USER CHARGING

Figure 3 :

 3 Figure 3: test vehicle from Q-Free

Figure 4 :

 4 Figure 4: Driven route in Frankfurt recorded by the reference receiver

Figure 5 :

 5 Figure 5: Same route in Frankfurt in a 3D building model

Figure 6

 6 Figure 6 displays the extreme performances of the uBlox receiver observed throughout these 28 revolutions downtown Frankfurt. Most trajectories typically show positioning error time series close to the best one, but 2 trajectories contain

Figure 6 :

 6 Figure 6: Position error of best (25) and worst (21) trajectories

Figure 7 :

 7 Figure 7: Description of a virtual gantry

Figure 8 :

 8 Figure 8: Virtual gantries designed in our experimental area in Frankfurt: the reference trajectory crosses through 20 virtual gantries In Figure 8 and in the following Figures 9 and 10, this color code has been applied: • green line: reference trajectory • blue line: real trajectory of uBlox • magenta line: cloned trajectory of uBlox

Figure 9 :Figure 10 :

 910 Figure 9: An example of missed virtual gantry (magenta line)

1

 1 NbCorrectDetections for VGi = RefNbCorrectDetection for VGi -Nb of missed detections for VGi RefNbCorrectDetection for VGi = N for all VG because each of them has been crossed only once C. Tuning of the simulation process A first set of error modeling parameters (denoted the default set) has been estimated with the best trajectory (25) using the Hutter method implemented by Monsifrot [6] [7] [12]. Afterward the upper bound of the truncation of the Laplace law is the only parameter we make varying, from 22 m to 400 m. Filter coef b { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1 }

Figure 11 :Figure 12 :

 1112 Figure 11: Time evolution of simulated and best target radii (trajectory 25)

  Figure 15: SA random draw for 100,000 samples trajectory

  Figure 17:

  Figure 18:

  Figure 19:

  Figure 20:

  Metrics for each Virtual Gantry with at least one missed or one false #VGi : Virtual Gantry Id #MISSED : Nb of missed dectections for VGi

	#Output RUC algorithm
	total missed number : 9
	total false number : 480
	#Global metrics
	Correct Charging Rate (CCR) : 99.99735%
	Over Charging Rate (OCR) : 2.82353%
	##CCER
	#FALSE : Nb of false dectections for VGi
	#FPCER
	VGi MISSED CCER	FALSE FPCER
	30279 0 100.00000% 210 1.23529%
	30275 0 100.00000% 145 0.85294%
	30273 0 100.00000% 59 0.34706%
	30283 0 100.00000% 27 0.15882%
	30276 0 100.00000% 39 0.22941%
	30113 4	99.97647% 0	0.00000%
	30112 3	99.98235% 0	0.00000%
	30099 1	99.99412% 0	0.00000%
	30111 1	99.99412% 0	0.00000%

TABLE III

 III

				.		VARYING RADIUS UPPER BOUND AND SIMULATION RESULTS
			SAi			Total		CCR	Total	OCR
						missed	(%)	false	(%)
						number		number
	SA0(11-36-22)		9			99.99	480	2.82
	SA1(11-36-44)		2062		99.39	7505	44.15
	SA2(11-36-66)		26631	92.16	13651	80.30
	SA3(11-36-88)		47105	86.15	17360	102.12
	SA4(11-36-176)	72672	78.63	21265	125.09
		120
		115
	Correct Charging Rate (CCR) (%)	85 90 95 100 105 110
		80
		75
		20	40	60	80	100	120	140	160	180
				max of Laplace for radius(m)

Evolution of CCR metric over max Laplace bound

	data
	fit using exponential

Acknowledgment

This study has been carried out in the frame of SaPPART COST Action, in particular during a Short Term Scientific Mission (STSM) executed by Miguel Ortiz from Ifsttar at Q-Free, using data sets acquired in the frame of the Norwegian SAVE project, dedicated to GNSS-based Road User Charging issues.


Philippe Bussetta

Narges Dialami

Michele Chiumenti

Carlos Agelet De Saracibar

Miguel Cervera

Romain Boman

Jean-Philippe Philippe Ponthot
email: jp.ponthot@ulg.ac.be

3D numerical models using a fluid or a solid formulation of FSW processes with a non-cylindrical pin

Keywords: Friction stir welding (FSW), Finite element method, Remeshing, Arbitrary Lagrangian Eulerian (ALE) formalism

3D numerical models using a fluid or a solid formulation of FSW processes with a non-cylindrical pin

Background

Friction stir welding (FSW) is a relatively recent welding process, which was developed at the Welding Institute (UK) and patented in 1991 [START_REF] Thomas | Friction stir welding for the transportation industries[END_REF]. FSW is a solid-state joining process. It means that during welding the materials to be joined are not melted. The joining is constituted by mechanical intermixing of the two materials. A rotating non-consumable tool is inserted between the two work-pieces and displaced along the welding direction (see Fig. 1). The tool is composed of two parts: a pin and a shoulder. The pin is introduced into the welded joint to mix deeply the two materials together. The aim of the shoulder is to contain the material around the pin. The part of the welding joint where the velocity of the tool and the advancing velocity add up is named the advancing side. The other part, where the two velocities are in opposite directions, is called the retreating side. The friction between the rotating tool and the work-pieces as well as the plastic deformation in the neighbourhood of the tool increase the temperature in the welded zone and thus soften the materials. But, during the process, the temperature is always smaller than the melting temperature of the materials. So, the heat-affected zone is smaller and the quality of the welding is higher than in more classical welding processes. In spite of the important number of applications of FSW, the phenomena happening during welding are still not very well understood. Therefore, the investigations on this process and especially regarding numerical simulations are still very active [START_REF] De Saracibar | Computational modeling and sub-grid scale stabilization of incompressibility and convection in the numerical simulation of friction stir welding processes[END_REF][START_REF] Feulvarch | A simple and robust moving mesh technique for the finite element simulation of friction stir welding[END_REF][START_REF] Dialami | Material flow visualization in friction stir welding via particle tracing[END_REF][START_REF] Assidi | Friction model for friction stir welding process simulation: calibrations from welding experiments[END_REF][START_REF] Santiago | Modelado termo-mecánico del proceso de friction stir welding utilizando la geometría real de la herramienta[END_REF][START_REF] Buffa | A new friction stir welding based technique for corner fillet joints: experimental and numerical study[END_REF][START_REF] Heurtier | Mechanical and thermal modeling of friction stir welding[END_REF]. From the numerical modelling point of view, a bibliography review is presented by He et al. [START_REF] He | A review of numerical analysis of friction stir welding[END_REF], by Dialami et al. [START_REF] Dialami | An apropos kinematic framework for the numerical modelling of friction stir welding[END_REF] and by Bussetta et al. [START_REF] Bussetta | Comparison of a fluid and a solid approach for the numerical simulation of friction stir welding with a non-cylindrical pin[END_REF]. Recently, the FSW process has been used to weld different material types, such as steel in Toumpis et al. [START_REF] Toumpis | Thermomechanical deformation behaviour of dh36 steel during friction stir welding by experimental validation and modelling[END_REF], magnesium in He et al. [START_REF] He | A multi-scale model for description of strain localization in friction stir welded magnesium alloy[END_REF] even polymers in Simões and Rodrigues [START_REF] Simoes | Material flow and thermo-mechanical conditions during friction stir welding of polymers: literature review, experimental results and empirical analysis[END_REF]. The mechanical intermixing induces very high strains in the material in the neighbourhood of the tool. Using a classical Lagrangian formulation as it is generally the case in solid mechanics would thus inevitably lead to mesh distortion. Consequently, classical numerical simulation techniques have to be extended in order to track the correct material deformations. One of the possibilities is to use the Arbitrary Lagrangian Eulerian (ALE) formulation [START_REF] Chiumenti | Numerical modeling of friction stir welding processes[END_REF][START_REF] Donea | Arbitrary Lagrangian Eulerian methods[END_REF][START_REF] Boman | Efficient ALE mesh management for 3D quasi-eulerian problems[END_REF][START_REF] Boman | Enhanced ALE data transfer strategy for explicit and implicit thermomechanical simulations of high-speed processes[END_REF]. This formulation is used to keep the mesh motion under control regardless of the real material displacements. The ALE formulation is also used to maintain a good mesh quality during the computation.

This paper deals with the extension to 3D of the works exposed in [START_REF] Bussetta | Comparison of a fluid and a solid approach for the numerical simulation of friction stir welding with a non-cylindrical pin[END_REF]. This article presents and compares two different 3D numerical approaches of the FSW process. The first model is based on a solid approach written in terms of nodal positions and nodal temperatures. The second model is based on a fluid approach written in terms of the velocity, the pressure and the temperature fields. Both models use advanced numerical techniques such as remeshing and the ALE formulation. 2D models are useful to test and easily compare both numerical formulations. Nevertheless, the FSW process is a fully 3D thermomechanical process. The effect of the shoulder and the thermal boundary conditions have a great influence on the FSW process, but they cannot be considered in the 2D models.

This article is split into three parts. First, both 3D numerical approaches are presented. Secondly, the 3D models are compared and the differences with the 2D models are exposed. Finally, some conclusions and explanations about the differences between the results of both models are made.

3D numerical modelling of FSW process

To model this welding process, the displacement of the tool is split into an advancing movement (actually assigned to the work-pieces but, in the opposite direction) and a rotation (imposed to the tool). In other words, the centre of the pin is fixed and a constant velocity is imposed to the plates (see Fig. 2). The tool is described by a classical Lagrangian mesh. Due to high deformations in the neighbourhood of the tool, the use of a Lagrangian formalism would lead very quickly to mesh entanglement. Then, the plates are modelled using the ALE formulation. On top of this, the ALE formulation allows the model to take into account tools with no rotational symmetry. In relation with the distance from the rotation axis of the tool, the plates are split in three zones. In the closest zone around the tool (red region in Fig. 3), the mesh has the same rotational speed as the tool. In the model, this region is limited by the value of the distance from the rotation axis of the tool equal to the value of the radius of the shoulder. In the furthest zone from the tool, the grey zone in Fig. 3, the mesh is fixed and the material is flowing through the mesh. The last zone is a transition zone, white region numbered 2 in Fig. 3. This zone connects the meshes of both other zones. Therefore, the quality of the [START_REF] Thomas | Friction stir welding for the transportation industries[END_REF], the transition zone corresponds to the white region (2), and the Eulerian formulation is applied on the grey region [START_REF] Feulvarch | A simple and robust moving mesh technique for the finite element simulation of friction stir welding[END_REF] mesh in this zone decreases with the simulation time. The numerical techniques used to overcome this problem are explained in the section named "The transition zone".

Thermomechanical formulation

The numerical models presented here are based on the finite element method. In this paper, two numerical formulations are compared (for more information about these formulations see [START_REF] Bussetta | Comparison of a fluid and a solid approach for the numerical simulation of friction stir welding with a non-cylindrical pin[END_REF]). The first one is based on a solid mechanics approach. It is written in terms of nodal positions and temperatures. The second one is based on a fluid mechanics approach. The equilibrium is written as a function of nodal velocities, pressures and temperatures.

Solid approach

The position and temperature fields are computed at each node of the elements. The mesh is composed of 21,980 linear hexahedral elements. The stresses and the internal variables are computed at each quadrature point of the element (8 Gauss points). To overcome the locking phenomenon, the pressure is considered constant over the element and computed only at a central quadrature point. The thermomechanical equations are split into a mechanical part and a thermal part. At each time step, the mechanical equations are first solved using a constant temperature field. This temperature field is the one obtained at the previous increment. Then, the thermal equations are solved on the frozen resulting geometrical configuration that has just been obtained.

Fluid approach

The fluid approach is based on a stabilized mixed linear temperature-velocity-pressure finite element formulation. This formulation is stabilized adopting the Orthogonal Sub-Grid Scale method (OSS) [START_REF] De Saracibar | On the orthogonal subgrid scale pressure stabilization of finite deformation J2 plasticity[END_REF][START_REF] Cervera | Mixed linear/linear simplicial elements for incompressible elasticity and plasticity[END_REF][START_REF] Chiumenti | A stabilized formulation for incompressible plasticity using linear triangles and tetrahedra[END_REF] to solve both the pressure instability induced by the incompressibility constraint and the instabilities coming from the convective term. A mesh of 74,127 linear tetrahedral elements is used for the domain discretisation. The velocity, the pressure and the temperature fields are computed at each node of the elements. The deviatoric stresses and the other internal variables are computed at each quadrature point of the element. Finally, the coupled thermomechanical problem is solved by means of a staggered time-marching scheme where the thermal and mechanical sub-problems are solved sequentially, within the framework of the classical fractional step methods [START_REF] De Saracibar | On the formulation of coupled thermoplastic problems with phasechange[END_REF][START_REF] Cervera | Thermo-mechanical analysis of industrial solidification processes[END_REF].

The transition zone

Solid approach

In the solid approach, the transition zone is a ring with a finite thickness (region 2 in Fig. 3). In this region, the evolution of the rotational speed of the mesh, which differs from the material velocity, is linearly interpolated between the ALE region and the Eulerian zone. As the mesh distortion grows with time, a remeshing operation is periodically required. The remeshing operation can be divided into two steps. First, a better-suited mesh, called the new mesh, is created. In this case, the simple geometry of this region allows an easy generation of the new hexahedral mesh. Then, to carry on the computation over this new mesh, the data are transferred from the old mesh to the new one (for more informations about the data transfer see [START_REF] Bussetta | Efficient 3D data transfer operators based on numerical integration[END_REF]).

Fluid approach

In the fluid model the transition zone (region 2 in Fig. 3) is limited to a tube of zero wall thickness. Each node of this zone is duplicated. One node is linked to the ALE region (numbered 1) and the other one to the Eulerian region (numbered 3). The coupling between both regions is performed using a specific node-to-node link approach. In this case the ALE mesh would slide precisely from one Eulerian interface node to the next one at each time step.

Thermomechanical constitutive model

In both models, the constitutive model of the tool is thermo-rigid. So, no mechanical fields are computed over this material. However, from the point of view of the thermal equations, the tool has a classical thermal behaviour as far as heat conduction is concerned. In addition, the material behaviour of the plates is modelled as thermo-visco-plastic using a Norton-Hoff constitutive model: where m and μ are the strain rate sensitivity and viscosity parameters respectively. Both are temperature dependent. S is the Cauchy stress deviator tensor and D is the deviatoric part of the strain-rate tensor.

S = 2μ(T ) D √ 3 2 3 D : D m(T )-1 (1
In the FSW process, the heat is mostly generated by the mechanical dissipation, which is computed as a function of the plastic strain rate and the deviatoric stresses as:

D mech = γ S : D (2)
where γ ≈ 0.9 is the fraction of the total plastic energy converted into heat.

Solid approach

In the solid model, the value of the variation of the pressure (dp) is computed thanks to the variation of the volume (dV ) and the bulk modulus (K ): dp = KdV . In addition, with the solid approach, it is possible to replace the Norton-Hoff constitutive model with a thermo-elasto-visco-plastic one, see e.g. [START_REF] Ponthot | Unified stress update algorithms for the numerical simulation of large deformation elasto-plastic and elasto-viscoplastic processes[END_REF]. With this kind of constitutive model, it is possible to compute the residual stresses.

Fluid approach

In the fluid model, the material is assumed to be incompressible and this constraint is incorporated into the equations to be solved. Fig. 8 Evolution of the pressure and the temperature computed by the two models at the control point 1 (rotating with the tool). The rotation speed of the tool is 40 RPM and the initial position of this point is defined in Fig. 5 Thermomechanical contact A perfect sticking thermomechanical contact is considered between the tool and the work-piece. It means that the temperature field and the displacement field are continuous through the interface between the tool and the work-piece. Like some authors [START_REF] Simar | Effect of rotational material flow on temperature distribution in friction stir welds[END_REF][START_REF] Gerlich | Strain rates and grain growth in Al 5754 and Al 6061 friction stir spot welds[END_REF][START_REF] Guerra | Flow patterns during friction stir welding[END_REF], we assume that the heat produced by the friction between the tool and the work-piece is negligible versus the heat generated by plastic deformations.

Comparison of numerical results

In this paper, the numerical results of the solid approach are compared with the already validated model based on the fluid approach (see [START_REF] De Saracibar | Computational modeling and sub-grid scale stabilization of incompressibility and convection in the numerical simulation of friction stir welding processes[END_REF][START_REF] Santiago | Modelado termo-mecánico del proceso de friction stir welding utilizando la geometría real de la herramienta[END_REF][START_REF] Dialami | An apropos kinematic framework for the numerical modelling of friction stir welding[END_REF][START_REF] Chiumenti | Numerical modeling of friction stir welding processes[END_REF]). In this example, the section of the pin is an equilateral triangle (Figs. 4,5). The dimensions of the tool are presented in Fig. 4. The width of the two plates is 50 mm, the thickness is 4.7 mm and the simulated length is 100 mm. The rotation axis of the tool is located at the centre of the simulated region (see Fig. 5). Fig. 9 Evolution of the pressure and the temperature computed by the two models at the control point 2 (rotating with the tool). The rotation speed of the tool is 40 RPM and the initial position of this point is defined in Fig. 5 The most important parameters of the considered FSW process are the rotation speed (40 RPM or 100 RPM) and the welding speed (400 mm min -1 ). The thermomechanical properties of the plates are the following:

• density: 2700 kg m -3 • bulk modulus: 69 GPa (used only with the solid approach) • thermomechanical Norton-Hoff law (presented in the page 4) with μ = 100 MPa, m = 0.12, • heat conductivity: 120 W m

-1 K -1 • thermal expansion coefficient: 1 × 10 -6 K -1 • heat capacity: 875 J kg -1 K -1
The thermomechanical properties of the tool are the following:

• density: 7800 kg m -3 ; • heat conductivity: 43 W m -1 K -1 ; Fig. 10 Evolution of the pressure and the temperature computed by the two models at the control point 3 (fixed in space). The rotation speed of the tool is 40 RPM and the initial position of this point is defined in Fig. 5 • heat capacity: 460 J kg -1 K -1 . During all the computation, the room temperature is considered constant at 20 • C. The thermal boundary conditions of this problem are the following (heat transfer is neglected along the external perimeter of the two plates):

• Conduction on the lower side of both plates (approximation of the thermal behaviour of the backing plate), exchange coefficient: 4500 W m -2 K -1 ; • Convection and radiation on the free upper side of both plates (except the part in contact with the tool), convection coefficient: 10 W m -2 K -1 , emissivity coefficient: 0.2.

With the rotation speed of 40 RPM, the total time of the simulation is 15 s which corresponds to 10 revolutions for the tool. With the rotation speed of 100 RPM, the solution is computed for 20 revolutions for the tool (i.e. the total time of the simulation is 12 s). Fig. 11 Final values of the pressure and the temperature computed by the two models along the control line (i.e. after 15 s). The rotation speed of the tool is 40 RPM and the considered line is defined in Fig. 5 Figures 6 and7 expose respectively the mesh of the solid model and the one of the fluid model. Figures 8,9, 10 and 11 show the evolution of the pressure and the evolution of the temperature computed by the two models with the rotation speed of 40 RPM at the control points and along the control line defined in Fig. 5. Figures 12,13, 14 and 15 present the same comparison with the rotation speed of 100 RPM. Points 1 and 2 have the same rotational speed as the tool (these points move according to the mesh). Point 3 is fixed in space.

After a transient phase which depends on the numerical strategy adopted for each approach the results of both models are very similar for the two values of the rotation speed of the tool (see Figs. [START_REF] Heurtier | Mechanical and thermal modeling of friction stir welding[END_REF][START_REF] He | A review of numerical analysis of friction stir welding[END_REF][START_REF] Dialami | An apropos kinematic framework for the numerical modelling of friction stir welding[END_REF][START_REF] Toumpis | Thermomechanical deformation behaviour of dh36 steel during friction stir welding by experimental validation and modelling[END_REF][START_REF] He | A multi-scale model for description of strain localization in friction stir welded magnesium alloy[END_REF][START_REF] Simoes | Material flow and thermo-mechanical conditions during friction stir welding of polymers: literature review, experimental results and empirical analysis[END_REF]. The difference of frequency between the pressure at point 3 and the pressure and the temperature at points 1 and 2 is explained by the fact that point 3 is fixed in space while points 1 and 2 have the same rotational velocity as the tool. On the one hand, the pressure at point 3 is affected by the three corners of the pin. On the other hand, the frequency of the pressure and the temperature at points 1 and Fig. 12 Evolution of the pressure and the temperature computed by the two models at the control point 1 (rotating with the tool). The rotation speed of the tool is 100 RPM and the initial position of this point is defined in Fig. 5 2 are controlled by the rotation speed of the tool. Consequently, the pressure frequency at point 3 is three times higher than the frequency of the pressure or the temperature at points 1 or 2.

In addition, the temperature field is a good indicator of the mechanical intermixing of the material. The small differences between the values of the temperature fields along the control line (Figs. 11,[START_REF] Chiumenti | Numerical modeling of friction stir welding processes[END_REF] show that the two models essentially deliver the same results. The good agreement between the temperature field is confirmed thanks to the value of the temperature field on the top of the plates (see Figs. [START_REF] Donea | Arbitrary Lagrangian Eulerian methods[END_REF], 17 for a rotation speed of 40 RPM and Figs. 18, 19 for a rotation speed of 100). The small gap between the temperature fields can be explained by the differences between the element types used by both models. Moreover, the value of the pressure field is similar around the tool (see Figs. 11,[START_REF] Chiumenti | Numerical modeling of friction stir welding processes[END_REF]. The difference between both formulations can explain an important part of the gap between the pressure field of both models. Indeed, the pressure is computed directly at the nodes Fig. [START_REF] He | A multi-scale model for description of strain localization in friction stir welded magnesium alloy[END_REF] Evolution of the pressure and the temperature computed by the two models at the control point 2 (rotating with the tool). The rotation speed of the tool is 100 RPM and the initial position of this point is defined in Fig. 5 with the fluid model and is evaluated at the location of the quadrature points with the solid model. Those values are then extrapolated to the nodes. These comparisons with two values of the rotation speed of the tool prove that the solid and the fluid models essentially deliver the same results as far as the temperature field and the pressure field prediction during welding are concerned. Nevertheless, from the point of view of CPU time, the fluid model is faster than the solid one (see Table 1). On the other hand, in the solid model the behaviour of the work-pieces can be modelled by a thermo-elasto-visco-plastic law. Then, the solid model is able to compute directly the residual stress. Figure 20 presents the value of the temperature and the von Mises stress over the control line of Fig. 21 with yield stress equal to 280 MPa and a rotation speed of 40 RPM at the end of the simulation (i.e. after 15 s). After a cooling period, the residual stresses are directly computed (see Fig. 22). The shape of the computed residual stresses corresponds to the classical shape observed after such a process [START_REF] Hattel | Modelling residual stresses in friction stir welding of al alloys-a review of possibilities and future trends[END_REF][START_REF] Sonne | The effect of hardening laws and thermal softening on modeling residual stresses in FSW of aluminum alloy 2024-T3[END_REF]. Fig. 14 Evolution of the pressure and the temperature computed by the two models at the control point 3 (fixed in space). The rotation speed of the tool is 100 RPM and the initial position of this point is defined in Fig. 5 The solutions with a rotation speed of 40 RPM are the extension in 3D of the results computed in 2D (see in [START_REF] Bussetta | Comparison of a fluid and a solid approach for the numerical simulation of friction stir welding with a non-cylindrical pin[END_REF]). The results computed by the 3D models are very different from the ones given by the 2D ones. The maximal temperature at points 1, 2 and 3 are respectively about 230 • C, 250 • C and 200 • C for the 3D models versus 140 • C, 110 • C and 105 • C in the case of the 2D models. These differences are explained by the approximations of the 2D models, like the effects of the shoulder and the thermal boundary conditions. Indeed, the shoulder that cannot be modelled in 2D has an important effect on the mechanical intermixing of the two materials, as well as on the value of the temperature field. The deformation rates in the neighbourhood of the shoulder are more important than the ones around the pin. Therefore, the main part of the heat generated during the FSW process is produced in the neighbourhood of the shoulder. Consequently, the heat generated by the shoulder explains that the values of the temperature fields computed [START_REF] Chiumenti | Numerical modeling of friction stir welding processes[END_REF] Final values of the pressure and the temperature computed by the two models along the control line (i.e. after 12 s). The rotation speed of the tool is 100 RPM and the considered line is defined in Fig. 5 with 3D models are more important than the ones of 2D models. So it is clear that a 3D model is more relevant than a 2D model.

Conclusion

The phenomena happening during the friction stir welding (FSW) process are at the interface between solid mechanics and fluid mechanics. In this paper, two different formulations are presented to simulate the FSW process numerically. One 3D model is based on a solid approach which computes the position and the temperature fields and another one is based on a fluid approach written in terms of velocity, pressure and temperature fields. Both models use advanced numerical techniques such as the Arbitrary Lagrangian Eulerian formalism or remeshing operations or an advanced stabilization algorithm. These advanced numerical techniques allow the simulation of the FSW process with a tool with no rotational symmetry. The aim of the paper is to compare two computational models based respectively on a solid and a fluid approach for the solution of FSW process. Based On the other hand, the model based on the solid approach has the advantage that it can be used with any thermo-elasto-visco-plastic constitutive model. Therefore, the solid model can be used to predict the FSW process and to also compute the residual stresses after the end of the process. [START_REF] De Saracibar | On the formulation of coupled thermoplastic problems with phasechange[END_REF] Final values after cooling of the temperature and the von Mises stress computed by the solid models along the control line. The behaviour of the work-pieces is modelled by a thermo-elasto-visco-plastic law and the considered line is defined in Fig. 21

Fig. 1

 1 Fig. 1 Scheme of the FSW process

Fig. 2

 2 Fig.2Scheme of the FSW model. A rotation is imposed to the pin (blue arrow) while the advancing velocity of the pin is replaced by a velocity imposed to the plates in the opposite direction (red arrows)

3 2 1 Fig. 3

 313 Fig.[START_REF] Feulvarch | A simple and robust moving mesh technique for the finite element simulation of friction stir welding[END_REF] Scheme of the different zones of the numerical model. The different zones of the model: ALE formulation is used on the red region (1), the transition zone corresponds to the white region (2), and the Eulerian formulation is applied on the grey region[START_REF] Feulvarch | A simple and robust moving mesh technique for the finite element simulation of friction stir welding[END_REF]

Fig. 4 Fig. 5

 45 Fig. 4 Scheme of the tool

Fig. 6 Fig. 7

 67 Fig. 6 Initial quadrangular mesh (21,980 elements) used in the solid model. Global view and zoom around the tool position

  Fig.[START_REF] Chiumenti | Numerical modeling of friction stir welding processes[END_REF] Final values of the pressure and the temperature computed by the two models along the control line (i.e. after 12 s). The rotation speed of the tool is 100 RPM and the considered line is defined in Fig.5

Fig. 16 Fig. 17

 1617 Fig.[START_REF] Donea | Arbitrary Lagrangian Eulerian methods[END_REF] Value of the temperature field computed by the solid model on the top of the plates at the end of the simulation (i.e. after 15 s). The rotation speed of the tool is 40 RPM

Fig. 18 Fig. 19

 1819 Fig. 18 Value of the temperature field computed by the solid model on the top of the plates at the end of the simulation (i.e. after 12 s). The rotation speed of the tool is 100 RPM

Fig. 20

 20 Fig.[START_REF] Cervera | Mixed linear/linear simplicial elements for incompressible elasticity and plasticity[END_REF] Final values of the temperature and the von Mises stress computed by the solid models along the control line (i.e. after 15 s). The behaviour of the work-pieces is modelled by a thermo-elasto-visco-plastic law and the considered line is defined in Fig.21

Fig. 21

 21 Fig. 21 Mesh of the solid model with the position of the control line used to plot the residual stress

Table 1 CPU time in hours versus the model type

 1

	Rotation speed	Solid model	Fluid model
	40 RPM	59	37
	100 RPM	131	79

Acknowledgements

The Belgian authors wish to acknowledge the Walloon Region for its financial support to the STIRHETAL project (WINNOMAT program, convention number 0716690) and to the FSW-PME project (Programme de recherche collective 2013 -convention number 1217826) in the context of which this work was performed.

Authors' contributions

The Belgian authors worked on the model based on the solid approach and the Spanish authors worked on the one based on the fluid approach. All authors have read and approved the final manuscript

Competing interests

The authors declare that they have no competing interests.