Vyshakh Palli-Thazha

David Filliat

Javier Ibañez-Guzmán

Applying map-masks to Trajectory Prediction for Interacting Traffic-Agents

Autonomous vehicles perceive their surroundings and foresee the future behaviour of all other relevant interacting traffic-agents in order to navigate safely and to operate in the public road networks. Trajectory prediction is difficult due to the stochastic manner these different types of trafficagents interact with each other and the changing navigation context. In this work, we propose a recurrent artificial neural network LSTM encoder-decoder based architecture to predict the movement of traffic agents. It uses map-masks of the area surrounding the ego vehicle and previous trajectory information to predict the trajectory of interacting traffic agents. This paper compares the proposed approach with LSTM baselines,using the NuScenes dataset which includes LiDAR point-cloud ground-truth data for traffic agents plus map information. Experimental results show that the proposed method outperforms the baselines based on the prediction accuracy.

I. INTRODUCTION

Sequence prediction problems have been an important area of research with multiple applications like Natural Language Processing, Weather Forecast, Economics and Intelligent Systems. Autonomous vehicle systems come across this problem when trying to understand the behaviour of all other traffic-agents that it interacts with. This could be a pedestrian, car, bike or any object that appear in its navigation space. Predicting the trajectory of the involved agents allow the autonomous vehicle to prepare itself for future decisions and manoeuvres, after assessing the risk involved in the current scenario. The goal is to try and make the vehicle understand the environment like humans do. The high level of interactions, occlusions and changing context makes it hard to solve this prediction problem with high accuracy. Therefore multiple inputs, patterns and behaviours must be studied to come to an operational and acceptable accurate prediction. Researchers have addressed the problem focusing on different aspects. Sequence predictions for single-agent behaviours have been solved using classical techniques but in real life this is unusually the case. There are multiple agents and their behaviour changes according to the class of the agent -whether it is a car, a pedestrian or a bike. Our research takes the class of the agent into account in the problem statement and tries to predict their different vyshakh.palli-thazha@ensta-paris.fr, vyshakh.v.palli-thazha@renault.com David Filliat, Director, U2IS Lab, ENSTA Paris, Institut Polytechnique de Paris david.filliat@ensta-paris.fr Javier Ibañez-Guzmán, Groupe Renault javier.ibanez-guzman@renault.com behaviours. Another factor that plays an important role in real life traffic scenarios are the interactions between different classes of agents. Our current research also gives importance to this and has incorporated an interactionaware model. Thirdly, our model takes into account, the surroundings of the ego-motion car, by using a map-mask that provides information on the driveable and non-driveable areas for different classes of agents. These map-masks also provide information on the path that could be taken by different agents. In the recent years, the boom in high-computational-power processors and GPUs has lead to faster research and development cycles. More and more data-sets are being released for different applications -KITTI, Apolloscape, NuScenes and NGSIM. This has lead to a focus on deep learning algorithms which solve classical problems with extremely high rate of success and accuracy. This success of Learning Algorithms have attracted research on several methods like Convolutional Neural Networks(CNN), Recurrent Neural Networks(RNN), Long Short-Term Memory(LSTM) and Generative Adversarial Networks(GANs). RNNs and LSTMs have shown promising results with solving sequence prediction problems. We use one such LSTM encoder-decoder architecture (see Figure 3 for an overview) as the baseline for our study and use the NuScenes dataset [START_REF] Caesar | nuscenes: A multimodal dataset for autonomous driving[END_REF] which provides 3D bounding-boxes in LiDAR point-clouds for different traffic agents, object/agent IDs and do a preliminary study for improving trajectory prediction using map-masks (binary maps) for the environment which we use in our experiments.

LiDAR Input

In this work, the use of map-mask patches to improve the prediction of trajectories for different classes of interacting traffic-agents is proposed. Specifically, our contributions include:

• A new LSTM encoder-decoder architecture that uses Map-Mask patches to make trajectory predictions for different classes of traffic agents in drivable and nondrivable areas • We evaluate the proposed model in comparison with LSTM baselines and is found to be superior in performance both in single-agent and multi-agent interaction scenarios. The remainder of the paper is organised as follows: Section II mentions relevant research papers and results that we have taken inspiration from or has guided our research. Research related to single-agent predictions, multiagent predictions, interaction-aware models, context-aware models are discussed in this section. Section III describes our map-mask based LSTM encoder-decoder architecture and the working of our complete trajectory prediction approach. Section IV discusses the results and the performance of our model and finally section V concludes the paper and mentions the future work.

II. RELATED WORKS A. Classical Methods

Trajectory prediction for pedestrians, vehicles and other traffic agents has been studied using classical techniques such as Constant Velocity, Constant Acceleration, Linear Regressions, Kalman Filters [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF], Monte Carlo Simulation [START_REF] Danielsson | Monte carlo based threat assessment: Analysis and improvements[END_REF], Time-series methods and Hidden Markov Models [START_REF] Firl | Predictive maneuver evaluation for enhancement of car-to-x mobility data[END_REF]. Most of these predictions are limited to short term predictions and even if these are long-term, their accuracy is affected because of the lack of context information and from ignoring the interaction between different classes of traffic-agents. A survey on motion prediction can be found in [START_REF] Lefèvre | A survey on motion prediction and risk assessment for intelligent vehicles[END_REF] where Lefèfre et. al. have analysed different trajectory prediction methods based on model completeness and risk assessment.

B. Learning based Models

Trajectory prediction methods have been applied to different sensor outputs, some solely with vision based detection outputs, others with 3D bounding-boxes obtained from Li-DAR point-clouds. We base our work with the detections on 3D point-cloud data. Much research has gone into developing deep learning algorithms to infer such predictions. Altché et. al. [START_REF] Altché | An LSTM network for highway trajectory prediction[END_REF] uses LSTMs to predict the future trajectory of a single target vehicle. It predicts the longitudinal and lateral trajectories of vehicles on a highway. In [START_REF] Milan | Online multi-target tracking using recurrent neural networks[END_REF] A. Milan et. al. proposed an end-to-end learning approach for online multi-target tracking which uses RNNs to predict the state of each target step by step and uses LSTMs to achieve data association. They work on vision data obtained from the 3D MOTChallenge dataset. Nikhil et. al. [START_REF] Nikhil | Convolutional neural network for trajectory prediction[END_REF] uses a Convolutional Neural Network based approach to predict trajectories of surrounding vehicles. This model utilises highly parallelisable convolutional layers to handle temporal dependencies. The trajectory histories are embedded to a fixed size tensor and stacked convolutional layers are used to ensure temporal consistency. Jawed et. al. [START_REF] Jawed | Datadriven vehicle trajectory forecasting[END_REF] also uses a similar architecture based on convolutional layers.

C. Interaction and/or Context Aware Models

There are several examples in the literature that tries to solve the prediction problem which involves interactions between traffic agents. Some of them focus on pedestrians alone or vehicles alone while others try to solve cases with interaction between multiple classes of agents.

Several models study interaction between vehicles: In [10] Hu et. al. uses a generative model to jointly predict the sequential motions of each pair of interacting vehicles and it uses as input past trajectories and environment information. Kim et. al. [START_REF] Kim | Probabilistic vehicle trajectory prediction over occupancy grid map via recurrent neural network[END_REF] uses occupancy grid and RNNs to predict the trajectory of vehicles on highways. In [START_REF] Deo | Convolutional social pooling for vehicle trajectory prediction[END_REF] Deo et. al. uses an an LSTM encoder-decoder architecture that uses convolutional social pooling for learning inter-dependencies between vehicles. Here the spatial configuration of the agents in a scene is embedded into a spatial grid around the egomotion vehicle and this grid is passed into convolutional and pooling layers to obtain the social context encoding. This along with an LSTM encoding of the agents' trajectory are concatenated and passed through a decoder to obtain the final predicted trajectory. A similar architecture is used by Park et. al. [START_REF] Park | Sequence-to-sequence prediction of vehicle trajectory via lstm encoder-decoder architecture[END_REF] with the decoder producing K likely trajectories over an occupancy grid using the beam search technique. Some models study pedestrian interactions: Kooij et. al. [START_REF] Kooij | Contextbased pedestrian path prediction[END_REF] proposes a Dynamic Bayesian Network which uses the situational awareness and spatial layout to predict pedestrian paths. It uses pedestrian head orientation, distance between vehicles and pedestrians at expected point of closest approach and distance of pedestrian to curbside as latent states on a Switching Linear Dynamical System (SLDS) to predict changes in pedestrian dynamics. In [START_REF] Alahi | Social lstm: Human trajectory prediction in crowded spaces[END_REF] Alahi et. al. tries to learn general human movement and predicts their future trajectories. They model pedestrian motion by a pooling based LSTM architecture and jointly predicts the trajectories of all the people in a scene.

Models also use various context information to improve trajectory prediction: Lee et. al. [START_REF] Lee | Desire: Distant future prediction in dynamic scenes with interacting agents[END_REF] applies a Conditional Variational AutoEncoder (CVAE) based RNN encoderdecoder to make prediction for interacting agents. They use past trajectories as input along with a feature map, generated from scene elements like roads and sidewalks, that provides context information to the prediction model. Habibi et. al. [START_REF] Habibi | Context-aware pedestrian motion prediction in urban intersections[END_REF] uses relative distance to curbside and status of pedestrian traffic lights as additional information to provide context to predict the pedestrian path. Ridel et. al. [START_REF] Ridel | Understanding pedestrian-vehicle interactions with vehicle mounted vision: An lstm model and empirical analysis[END_REF] implicitly models pedestrian interactions with vehicles to predict pedestrian behaviour. They use pedestrian head orientation as an input along with past trajectories. Another interesting study is done by Ma et. al. [START_REF] Ma | Trafficpredict: Trajectory prediction for heterogeneous traffic-agents[END_REF] which uses 4D graphs to model interactions and classes of the interacting agents.

Our approach addresses the problem by using map-masks to facilitate the understanding of the scene environment. This study is a first step towards using semantic information in the form of maps for improving trajectory prediction. The trend is to add more information to maps. Works like OpenStreetMaps and NuScenes are releasing maps with information other than road geometry. Details such as crosswalk, sidewalk, traffic lights, stop lines, lanes are being added to maps. These semantics could be used to improve the perception of the environment.

III. PROBLEM FORMULATION AND METHOD

In a driving scenario, several objects might be present, sharing the same ego-vehicle work-space it navigates. These include pedestrians, other vehicles or cyclists. State-of-the-art perception algorithms provide bounding boxes/segmentation for each class of objects. There are algorithms which can perceive the said objects from camera images or from LiDAR point-cloud data. In our case, we use LiDAR point-cloud data and off-the-shelf detection algorithms to extract the features of all the traffic-agents within the field of view of the sensor used. For experimental purposes, the PointPillars object-detection network by Lang et. al. [START_REF] Lang | Pointpillars: Fast encoders for object detection from point clouds[END_REF] is adopted. This network is chosen as it shows superior performances over other relevant frameworks and they have released their code on GitHub. The observed trajectory is taken from the output of the object detection framework. For training the prediction network, the ground truth trajectory available for the situation from the NuScenes dataset is used.

The feature set of each traffic agent k at time t can be defined by

f k t = (p k t , c k t , M k t ) where p k t = [x k t , y k t ],
x and y are the spatial coordinates of the detected traffic-agent and c is the class of the traffic agent, c i ∈ (1, 2, 3) where 1 is for pedestrians, 2 is for cars and 3 is for cyclists and M k t is a 128*128 pixels binary map patch of the area surrounding each detected object. The task is to observe the features of all traffic-agents in the interval [1 : T obs ] and predict their positions in [T obs+1 : T pred ].

A. Trajectory Prediction using LSTMs

The prediction framework developed by T. Fernando et. al. [START_REF] Fernando | Soft + hardwired attention: An lstm framework for human trajectory prediction and abnormal event detection[END_REF] that applies LSTMs for predicting trajectories is applied, along with the training methods adopted by Bahdanau et. al. in their Neural Machine Translation work [START_REF] Bahdanau | Neural machine translation by jointly learning to align and translate[END_REF]. For each time-step, an object-pool is created, associated with a trajectory prediction module. It regresses short term and long-term trajectories for each traffic-agent in the objectpool. This is based on the historic trajectory of that agent and neighbouring trajectories. Therefore, each traffic-agent that resides in the pool is associated with a probable short term and long-term trajectory. This approach is advantageous because raw detection outputs can be noisy due to occlusions, false alarms, inaccurate bounding boxes, and missing detec-tions. Studies show that occlusions last around 3-10 frames [START_REF] Sadeghian | Tracking the untrackable: Learning to track multiple cues with long-term dependencies[END_REF]. Let the historical trajectory of a traffic-agent k, from timestep 1 to time-step T obs be given by,

p k = [p 1 , ..., p T obs ],
where,

p i = [x i , y i ]
for each time-step i. p i composes of points in a Cartesian grid obtained from the top-down view of the LiDAR pointcloud detections. These historical trajectories are passed through the LSTM encoder of each respective traffic-agent to generate its historical embeddings as follows,

h k t = LST M (p k t , h k t-1
), generating a sequence of historical embeddings.

A historical context vector C H,k t is defined to encode the trajectory information from the traffic-agent of interest (k), which can be computed as a weighted sum of hidden states from t = [1 :

T obs ], C H,k t = T obs j=1 α tj h k j
and the weight α tj as shown in [START_REF] Bahdanau | Neural machine translation by jointly learning to align and translate[END_REF] can be computed by,

α tj = exp(e tj ) T l=1 exp(e tl ) e tj = a(h k t-1 , h k j )
where the function a is a feed forward neural network which is trained jointly with the whole network.

The spatial context vector C S,k is used for embedding the neighbouring trajectories. The spatial weights, denoted by w n j , can be computed as,

w n j = 1 dist(n, j)
where dist(n, j) is the Euclidean distance between the n th neighbour and the traffic-agent of interest at the j th time-step, and w n j is the generated spatial attention weight. When there are N neighbouring trajectories in the local neighbourhood, and h n j is the encoded hidden state of the n th neighbour at the j th time-step, then the context vector for the spatial model is defined as,

C S,k = N n=1 T obs j=1 w n j h n j
The merged context vector, C * ,k t , computed by,

C * ,k = tanh([C H,k ; C S,k ])
is then passed through the LSTM decoder to predict the future trajectory for the traffic-agent of interest. This is explained in Section III-C. It is by modifying the context vector that the information from the binary map-masks is input to the prediction model. This is explained in the next sub-section. Fig. 3: Passing the map-patches through a CNN AutoEncoder generates the latent vector that stores information about the drivable area and this vector is passed on to the LSTM encoder architecture For understanding the behaviour of different classes of traffic-agents in drivable and non-drivable spaces, the input of map information in the form of binary patches is studied. We take square patches of the size 128*128 pixels (see Figure 2) around each detected traffic-agent in the perception space (see Figure 2a) of the ego-motion vehicle and for each timestep in the observed time interval. Each of these patches are passed through a Convolutional AutoEncoder producing a latent vector C C t which holds the information about the drivable space associated with the respective traffic-agent (see Figure 3). This vector is then concatenated with the merged context vector C * t from the LSTM Model to obtain the new latent vector C M,k t .

B. Map Input for Trajectory Prediction

C. LSTM Decoder

Now, the final predicted trajectory q t is obtained by passing the context vector C M,k t in the case of map input or C * ,k t in the case of only LSTMs, through the LSTM decoder and two fully connected layers. The overall system architecture is illustrated in Figure 1.

q t = LST M c (h k t-1 , q t-1 , C M/ * ,k t
) q k t is composed of points in a Cartesian grid. The decoder used for each traffic-agent class is different. Hence, based on c, we choose the LSTM to be used to decode. We let T obs = 3 for short term predictions and T obs = 10 for long term trajectory predictions.

q t = [x Obs+1 , y Obs+1 , ..., x Obs+N , y Obs+N ]
where N is the prediction horizon.

The CNN AutoEncoder is tested both with and without pre-training and no significant difference has been found. Therefore, this model is attached to the LSTM encoderdecoder and trained end-to-end with Adam optimiser and a learning rate of 1 × e -4 for 100 epochs and fine-tuned with a learning rate of 1 × e-5.

IV. RESULTS AND EVALUATION

The prediction accuracy is measured based on two criteria as described in [START_REF] Gupta | Social gan: Socially acceptable trajectories with generative adversarial networks[END_REF], [START_REF] Zhao | Multi-agent tensor fusion for contextual trajectory prediction[END_REF]: Average Displacement Error (ADE) and Final Displacement Error (FDE) in pixels with respect to each time step t within the prediction horizon:

1) Average Displacement Error (ADE): Average distance between ground truth and our prediction over all predicted time-steps for the traffic agent.

ADE(i) = 1 T j=1,2...T (x j i -x j i ) 2 -(ŷ j i -y j i ) 2 ADE = 1 n i=1,2,...n ADE(i)
2) Final Displacement Error (FDE): The distance between the predicted final destination and the true final position of the agent.

F DE(i) = 1 n (x T i -x T i ) 2 -(ŷ T i -y T i ) 2 F DE = 1 n i=1,2,...n F DE(i)
where n is the total number of interacting agents in the test set, x ij and y ij denote the coordinates of the i th agent in the predicted time-step j and T denotes the final predicted time-step.

As training data, the NuScenes dataset is used. From this we obtain 3200 different trajectories with enough length (at least 30 time frames) for the experiments. These trajectories are obtained from 850 scenes of 20s recordings from the NuScenes data-collection car setup. These are split into train, validation and test sets of size 2100, 550 and 550 trajectories respectively. For single-agent trajectory prediction, these sequences are considered separately, and for the multiagent trajectory predictions, the information of all interacting agents are taken into account. Each trajectory sequence is associated with the relevant interacting agent information also.

In order to illustrate the interest of using drivable/nondrivable map information, four variants of our approach are implemented. The performance of the model in different scenarios is shown in Figure 4. Figure 4a shows an example of straight trajectories and the prediction giving results in the nondrivable area when only LSTMs are used and better predictions confined to the drivable area when map-mask is added as input. Also the fitting of the trajectory prediction to road curves is better predicted when map-masks are introduced as can be seen in 4b. Another example of prediction being made in the drivable area is seen in Figure 4c.

The comparison of the ADE/FDE values for the LSTM baselines and the models with map-mask input is presented in Table I.

As can be seen in each case, single agent and multiagent predictions, the addition of the map-mask patches The green path is the Ground truth and the blue path is the predicted trajectory help the system make more accurate predictions. This gives us confidence to study the use of a map with much more information in the future to improve the predictions even more.

V. CONCLUSION AND FUTURE WORK

The proposed approach discusses and experimentally verifies the idea of using binary map-masks to improve the prediction of trajectories for different traffic agents. It compares single-agent and multi-agent trajectory prediction models and asserts the overall improvement of the prediction when mapmasks are given as input. This study serves as a preliminary step in the direction of using more environmental context information for improving prediction accuracy. As a future work, we plan to study OpenStreetMaps and use the multitude of information that can be obtained from their wide network of mapped areas.

Fig. 1 :

 1 Fig. 1: System Architecture: The prediction module has two inputs: sequential historical trajectories of all traffic-agents and the associated map-mask patches of each agent for the observed time instances

Fig. 2 :

 2 Fig. 2: Binary map-mask and map patches: A trajectory and the map-patches surrounding the agent are marked by blue squares, centred around the traffic agent marked in green

• Method 1 :

 1 Single Agent LSTM Model for all traffic agents: Separate LSTM models are used to predict the motion of single traffic agents. No interaction is involved. • Method 2: Single Agent with LSTM Model and Map-Mask input: Map-mask is taken as input to predict the trajectories of single agents, not considering interactions. • Method 3: Multi-Agent interactions and Model with LSTM only: Interactions are taken into consideration and only LSTM encoder-decoder architecture is used. This is similar to the method in [26]. • Method 4: Multi-Agent interactions and Model with LSTM and Map-Mask input: Multi-agent interactions are taken into account in the LSTM encoder-decoder architecture and map-masks are used to improve the predictions of traffic-agents involved in such interactions.

Fig. 4 :

 4 Fig. 4: Behaviour of Trajectory Prediction Models without(Left) and with Map-mask(Right) input. The green arrows in the images show the direction of prediction. The green path is the Ground truth and the blue path is the predicted trajectory

TABLE I :

 I Comparison of Prediction Accuracy for different models


domain: spi phys chim

text:
Véronique Kehr-Candille
email: veronique.kehr-candille@onera.fr

Modelling the damping at the junction between two structures by non-linear models: improving the model and the resolution

We are interested in the modelling of the damping at the junction between two substructures. We model the connection by a meta-model which takes into account both dissipative and non-linear aspects of the connection. In this work we use a Iwan-Jenkins model, inserted in a finite element system, and we solve the time-domain system with an algorithm using the notion of sub-differential. This algorithm is compared with other more usual time-domain solvers.

In structural dynamics, vibratory levels depend directly on damping. So it is necessary to have at its disposal, tolls and models allowing to correctly represent the damping from the design phase.

The origin of energy dissipation in structures is double: on one hand, materials intrinsic damping, and on the other hand, dissipation generated by the friction phenomena at the interfaces between the sub-structures. We are interested here in the modelling of this second source of dissipation. For the metallic structures, the dissipation generated at the interfaces is the main damping source.

Modelling friction and interface contact has been the subject of numerous studies. Most finite element softwares contains contact modules (Nastran, Abaqus, Aster). However, these approaches, used for the calculation of dynamic answers, lead to prohibitive calculation times. Our objective is to develop a relatively simplified model (called meta-model) of the junction between two sub-structures, which allows to correctly represent dissipative and non-linear aspects, without leading to excessively long time in the calculations of dynamic answer.

Previous work 1.2

In previous studies [START_REF] Kehr-Candille | Modelling the Damping at the Junction between Two Substructures[END_REF], we have developed such a meta-model: an equivalent Bouc-Wen model or Dahl model is inserted between two degrees of freedom, in a linear finite element system. The final system is nonlinear and is solved with classical Runge-Kutta algorithms, with the help of Craig-Bampton reduction methods. In parallel, methods for identifying Bouc-Wen or Dahl model parameters from the hysteresis curve have been developed [START_REF] Kehr-Candille | Identification of hysteretic systems for damping modelling[END_REF].

We have shown that this approach allowed to correctly representing both the non-linear and dissipative behaviors of the junction, without leading to excessively long computational times.

Current work 1.3

The results obtained in the previous studies were not completely satisfactory. Indeed, in many experimental investigations (cf. figure 1, extracted from [START_REF] Kehr-Candille | Modelling the Damping at the Junction between Two Substructures[END_REF]), we observe on the frequency response curves the following non-linear trend : when the amplitude of the excitation increases, the resonance frequency decreases (the curved is shifted on the left) and damping increases (the resonance peak is lower). With the Bouc-Wen or Dahl models, the second trend (damping increases with amplitude) could have been reproduced, but the first trend (resonance frequency decreases with amplitude) could not be simulated (in numerical simulations, on the opposite, sometimes resonance frequency increases with amplitude). Moreover, it has been shown [START_REF] Lamarque | Non-smooth Deterministic or Stochastic Discrete Dynamical Systems[END_REF][START_REF] Bastien | An analysis of the modified Dahl and Masing models : Application to a belt tensioner[END_REF] that the equations obtained with the Iwan-Jenkins model can be solved with an original and promising algorithm, developed for non-smooth model.

Resolution of non-smooth models 2

Differential inclusion 2.1

The concept of differential inclusion (maximal monotone operator), the definition of the resolvent function and the corresponding computational algorithm are extensively described in [START_REF] Lamarque | Non-smooth Deterministic or Stochastic Discrete Dynamical Systems[END_REF], with rigorous mathematical demonstrations. We present here a brief survey of these notions, and a description of the algorithm.

A differential inclusion is a multi-valued operator whose graph is maximal monotone. For instance, the operator 𝜎 defined by equation ( 1) is a differential inclusion (cf. figure 2 left).

(1)

An other example is the 𝛽 operator, defined by equation ( 2), which can be considered as the inverse of operator 𝜎 (cf. figure 2 right).

(2)

Figure 2 -Multi valued operators

Classical algorithms used for solving non linear models (Runge-Kutta, Newmark, …) rely on the concept of Taylor expansion, and therefore must be used only for at least 𝐶 1 operators. Strictly speaking, for nonsmooth operators like 𝜎 and 𝛽, these numerical schemes cannot be used.

Reference [START_REF] Lamarque | Non-smooth Deterministic or Stochastic Discrete Dynamical Systems[END_REF] presents a method for solving these non-smooth operators, if they can be expressed as maximal monotone operators.

For a differential inclusion 𝐴, the function 𝐽 𝜆 = (𝐼𝑑 + 𝜆 𝐴) -1 is called the resolvent function. Now, we consider the differential system:

𝑈 ̇= 𝑓(𝑡, 𝑈) where 𝑓(𝑡, 𝑈) corresponds to a differential inclusion. This system can be solved in the time domain by the Euler implicit scheme: [START_REF] Lamarque | Non-smooth Deterministic or Stochastic Discrete Dynamical Systems[END_REF] where 𝐽 ℎ is the resolvent function and ℎ is the time interval.

It is shown in [START_REF] Lamarque | Non-smooth Deterministic or Stochastic Discrete Dynamical Systems[END_REF] that the numerical solution obtained with this algorithm exists and is unique.

Iwan-Jenkins model 2.2

The Jenkins model, also called Masing model or Prandlt model, is the association in sequence of a spring (stiffness 𝑘) and a dry friction element (threshold 𝐴). The Iwan model is the association in parallel of several Jenkins models. With the addition of a spring (stiffness 𝑘 0 ) in parallel, we obtain the generalized Iwan model (cf. figure 3). For the elementary Iwan model, the differential system is written :

(4)
denoting 𝑓 the restoring force (other quantities are defined in figure 3). For the dry friction element, the Coulomb friction is written using the differential inclusion 𝜎 defined by equation (1). This system can be set under the form 𝑈 ̇= 𝑓(𝑡, 𝑈) with [START_REF] Lacayo | Nonlinear modeling of structures with bolted joints: A comparison of two approaches based on a time-domain and frequency-domain solver[END_REF] In the same way, the differential system for the generalized Iwan model is: [START_REF] Festjens | A numerical toll for the design of assembled structures under dynamic loads[END_REF] This model has been widely studied in literature. We show in figure 4 the hysteresis cycles (nonlinear force versus displacement) obtained for a Jenkins model and for an Iwan model with two elementary models. With an increasing number of elementary models in parallel, a lot of hysteresis cycles can be reproduced with the generalized Iwan model. Reference [START_REF] Bastien | An analysis of the modified Dahl and Masing models : Application to a belt tensioner[END_REF] presents a method for the identification of the generalized Masing model from experimental hysteresis curve. Our purpose is to represent the junction by a simple non linear model (meta model) inserted between two degrees of freedom of an initial linear finite elements system (see figure 5). The method is the same as the method used in [START_REF] Kehr-Candille | Modelling the Damping at the Junction between Two Substructures[END_REF].

(7)

In [START_REF] Kehr-Candille | Modelling the Damping at the Junction between Two Substructures[END_REF], we used a Bouc-Wen model, and here we use an Iwan model instead. So, the differential system is written as follows:

(8)

The Iwan model introduces as many additional degrees of freedom as elementary Jenkins models. Algorithms 2.4

Differential inclusion

The operator corresponding to the elementary Iwan model is maximal monotone, and therefore the numerical described in 2.1 can be used.

Moreover, in this case the resolvent function is explicit (resolvent function associated to operator 𝜷).

According to [START_REF] Lamarque | Non-smooth Deterministic or Stochastic Discrete Dynamical Systems[END_REF], we have:

(9)
So, for a one degree of freedom system, we have:

(10)
And for a multiple degree of freedom system (finite element system as in 2.3), we have:

(11)
In reference [START_REF] Lamarque | Non-smooth Deterministic or Stochastic Discrete Dynamical Systems[END_REF], it is shown that the solution to the problem with differential inclusion exists and is unique for a one degree of freedom system. We have shown (demonstration not detailed here) that this result of existence and unicity could be generalized to our case.

Runge-Kutta

The multi values operator 𝜎 can be regularized using the following approximation:

(12) with coefficient 𝐾 having a sufficient "high" value (cf. figure 6). And for a multiple degrees of freedom system, the system (7) became: (14) These systems involve now continuous functions, and can be solved using classical non-linear algorithms like Runge-Kutta schemes.

Applications 3

One degree of freedom 3.1

The algorithm using differential inclusions has been applied first to a one degree of freedom system of mass 𝑚 = 1 kg.

First we use an elementary Jenkins model with stiffness 𝑘 = 10 N.m, and different values for threshold 𝐴. The excitation force is a swept sine of amplitude 𝐹 0 and speed 0.2 octave/minute, between 0.3 Hz and 0.6 Hz. Figure 7 on the left shows the time response curve for different values of threshold 𝐴 : a decreasing value of 𝐴 leads to an increasing damping. Figure 7 on the right shows normalized frequency responses for 𝐴 = 50 and this highlights the non-linear behaviour of the Jenkins model (damping increases with amplitude).

Hysteresis curves are plotted for a sinusoidal excitation (frequency closed to the eigen frequency). For comparison, we also used a more classical algorithm to solve the regularized system with the help of equation ( 12). We use a Runge-Kutta-Fehlberg method (order 4, error estimator of order 5) referenced in the following RKF45. The two algorithms (resolvent method with differential inclusion and RKF45) are compared in terms of results in figure 9 and in terms of computation times in table 1, for an elementary Jenkins model (𝑘 = 10 𝑁. 𝑚, 𝐴 = 2) with a swept sine excitation. The time step is 10 -3 s. Concerning computation times, it is clear that the resolvent method is much faster than the RKF45 algorithm. Concerning the accuracy of the results obtained by the RKF45 method, the choice of the value of the regularization parameter 𝐾 is crucial, but this choice is not at all obvious.

However, the explicit Euler scheme used with the resolvent method needs small or even very small time steps (this is a well known drawback). On the other hand, the Runge-Kutta method doesn't require so small time step and moreover, the time step can be adjusted automatically. So, in some cases, the advantage of the resolvent method in terms of computation time is not so obvious. But its advantage in terms of accuracy of the results remains (the unicity has been proved mathematically).

Multiple degrees of freedom 3.2

The method has been applied to a multiple degrees of freedom system : the assembly of two tension bars by an Iwan-Jenkins model (see figure 10). The system is equivalent to the assembly of springs (stiffness 𝑘 = 4. 10 4 N. m) and masses ( 𝑚 = 10 kg). The normalized frequency response (FRF) at the end point is plotted figure 11. The non-linear behavior of the Iwan model is highlighted: with increasing excitation amplitude, the resonance frequency decreases, and the damping increases. This simulation result is consistent with usual experimental results.

The results plotted figure 11 has been obtained for a small value of time step: 5. 10 -6 s, needed for convergence of Euler explicit scheme. We tried to solve with the RKF45 algorithm with an adaptive time step for various values of coefficient 𝐾 (cf. equation ( 12)), but we failed to reproduce the results obtained with the resolvent method (cf. figure 12). We have improved the previous works concerning the modelling of the damping at the junction by a meta-model. We used an Iwan-Jenkins model which allows reproducing non-linear behaviours observed in experimental tests. The numerical computation has been improved by an algorithm using differential inclusions.

This method has been applied to an academic example. The results are encouraging to apply it to an industrial structure. However, the Euler scheme used by the computation method needs small time steps, and this is an important drawback of the method. This point must be investigated in the future.

Figure 1 -

 1 Figure 1 -Experimental results

Figure 3 -

 3 Figure 3 -Jenkins model (left) and generalized Iwan model (right)

Figure 4 -

 4 Figure 4 -Theoretical hysteresis curve for the Iwan-Jenkins model

Figure 5 -

 5 Figure 5 -Meta-model inserted in a finite element system

Figure 6 -

 6 Figure 6 -Regularization of multi valued operator 𝜎

  Figure 8 shows two curves: one for an elementary Jenkins model referenced 𝑁 = 1 (𝑘 = 10 𝑁. 𝑚, 𝐴 = 5), and the other one for an Iwan model with five elementary models referenced 𝑁 = 5 (𝑘 1 = ⋯ = 𝑘 5 = 𝑘 5 ⁄ = 2 𝑁. 𝑚 and [𝐴 𝑖 ] = [1,2,3,5,7]).

Figure 7 -

 7 Figure 7 -Time response (left) and frequency response (right)

Figure 9 -

 9 Figure 9 -Comparison of algorithms

Figure 10 -

 10 Figure 10 -Finite elements system

Figure 11 -

 11 Figure 11 -FRF at the end point

Table 1 :

 1 Computation times

	computation	resolvent algo.	RKF45
	time (s)	0.06	K=20 K=200	0.22 0.42


text:
Keywords: Timoshenko beam, non-selfadjoint, Riesz basis

. The outline of the proof is as follows. First we write a variational formulation of the problem which is transformed into a first order evolution equation. The unbounded operator of this evolution equation is the generator of a contraction semigroup so that the problem is well-posed. This operator has a compact resolvent but it is non-self-adjoint. The spectrum of this operator is studied by making a change of variable in order to reduce the problem to a first-order system asymptotically linear in the spectral parameter. By applying some results of [2], we obtain the asymptotics of the spectrum of this operator and we show that its eigenvalues are asymptotically algebraically simple. Then we show that the system of root vectors of the operator is complete and, applying a theorem from [3], that they form a Riesz basis. As a consequence of this result, one obtains an asymptotic expansion of the solution of the initial problem.

Spectral analysis of a nonhomogeneous rotating Timoshenko beam: Riesz basis property

Jean-Luc Akian 1, * 1 DMAS, ONERA, Université Paris Saclay F-92322 Châtillon -France

Introduction

In what follows, l is the length of the beam, A its cross-sectional area, E its Young modulus, ρ its mass density, I the moment of inertia of its cross-section, G its shear modulus, α s its shear correction factor, w its infinitesimal transverse displacement, ϕ the infinitesimal rotation of its cross-section, k w , k ϕ , c w and c ϕ the corresponding spring and damping con-stants at the junction. All the coefficients and their inverses are assumed to be in W 1,∞ (0, l). If u 0 , v 0 are given in subspaces of [H 1 (0, l)] 2 and [L 2 (0, l)] 2 , the corresponding variational problem reads as follows: find

u = t (w, ϕ) ∈ C 1 ([0, +∞[; [H 1 (0, l)] 2 ) ∩ C 2 ([0, +∞[; [L 2 (0, l)] 2 ) such that ∀δu ∈ [H 1 (0, l)] 2 , m(ü, δu) + b( u, δu) + a(u, δu) = 0 and such that u(0) = u 0 ; u(0) = v 0 (Problem (P )), where ∀ u, δu ∈ [H 1 (0, l)] 2 , a(u, δu) = l 0 [T w δw + EIϕ δϕ ]dx + l 0 α s GA× [(w -ϕ)(δw -δϕ)]dx + k w w(0)δw(0)+ k ϕ ϕ(0)δϕ(0) -Ω 2 l 0 ρIϕδϕdx, b(u, δu) = c w w(0)δw(0) + c ϕ ϕ(0)δϕ(0), and ∀ u, δu ∈ [L 2 (0, l)] 2 , m(u, δu) = l 0 ρIϕδϕdx + l 0 ρAwδwdx. Setting v = t (z, ψ), U = u v , U 0 = u 0 v 0 -Ωu 0
, after replacing u(t) by e Ωt u(t), (P ) boils down to

U = AU, U(0) = U 0 , where A is an unbounded operator of H = [H 1 (0, l)] 2 × [L 2 (0, l)] 2 with domain D(A). A * is an unbounded operator of H with a do- main D(A * ) = D(A) so that A is non-self- adjoint.
On the other hand 0 ∈ ρ(A) and A -1 is a compact operator of H therefore the spectrum of A consists of isolated eigenvalues of finite algebraic multiplicity. A is maximaldissipative hence A is the generator of a contraction semigroup on H. Returning to the initial problem (P ) it follows that if U 0 ∈ D(A), there exists a unique solution to

(P ) such that u ∈ C 0 ([0, +∞); [H 2 (0, l)] 2 ) ∩ C 1 ([0, +∞); [H 1 (0, l)] 2 ) ∩ C 2 ([0, +∞); [L 2 (0, l)] 2 ) and such that U ∈ C 0 ([0, +∞); D(A)).
Suggested members of the Scientific Committee: Vincent Pagneux, Cédric Bellis

2 Spectral theory of A Setting α 1 = ρA T +αsGA , α 2 = ρ E , Λ1 = diag (α 1 , α 2 , -α 1 , -α 2 ) there exist matrices T 0 , Λ0 , Λ-1 , B1 , B0 , C1 , C0 such that for λ ∈ C * , the conditions: U ∈ D(A), U = 0 and (A -λId)U = 0 are equivalent to the condi- tions: Z ∈ [H 1 (0, l)] 4 , Z = 0 and T (λ)Z = 0 where Z = T -1 0 u u /λ , T (λ) = T D (λ) T R (λ) , T D (λ)Z = Z -(λ Λ1 + Λ0 + 1 λ Λ-1 )Z, T R (λ)Z = ( B1 + 1 λ B0 )T 0 Z(0)+( C1 + 1 λ C0 )T 0 Z(l).
Assume that there exists a constant

C > 0 such that |α 1 -α 2 | ≥ C > 0 on [0, l]. Then from [2],
Theorem 2.8.2, there exists an asymptotic fundamental matrix function of T D (λ)Z = 0 denoted Ẽ(., λ) such that for λ ∈ C * , the conditions: 0) and c ϕ = (α 2 EI)(0) the eigenvalues of A are given asymptotically by the formulas If Ω is an open nonempty subset of C, E and F are Banach spaces, µ ∈ Ω, S ∈ H(Ω, L(E, F )) and x ∈ H(Ω, E) are holomorphic families, x is called a root function of S at µ if x(µ) = 0 and (Sx)(µ) = 0. The order of the zero of Sx at µ is called the multiplicity of x with respect to S at µ. One can also define the geometric and algebraic multiplicities of S at µ [2], Chapter I. By comparing the multiplicities of the root functions and the geometric multiplicities of the three holomorphic families

Z ∈ [H 1 (0, l)] 4 , T D (λ)Z = 0 are equivalent to the condition: ∃ c(λ) ∈ C 4 with Z(., λ) = Ẽ(., λ)c(λ). It follows that λ is an eigenvalue of A iff detM (λ) = 0 where M (λ) = T R (λ) Ẽ(., λ). If c w = (α 1 (T + α s GA))(
λ 1 n = µ 1 -Ω+ν 1 n , λ 2 n = µ 2 -Ω + ν 2 n where µ 1 , µ 2 ≤ 0, ν 1 n , ν 2 n ∈ iR, |ν 1 n |,
λ ∈ C * → T (λ) ∈ L([H 1 (0, l)] 4 , [L 2 (0, l)] 4 ×C 4 ), λ ∈ C * → A(λ) = A -λId ∈ L(D(A), H) and λ ∈ C * → M (λ) ∈ L(C 4 , C 4
) at an eigenvalue of A, it can be shown that if µ 1 = µ 2 , the eigenvalues of A are asymptotically algebraically simple (in the "classical" sense). By estimating the resolvent of A first on the negative real axis (comparing A with the self-adjoint operator obtained by taking c w = c ϕ = 0) and then on the right half-plane (using the Hille-Yosida Theorem) and by applying the Phragmen-Lindelöf Theorem to a certain function of exponential type, it is shown that the system of root vectors of A is complete in H.

Recall that a sequence (h n ) n∈N of a Hilbert space V is a Riesz basis of V iff there exists an isomorphism U of V and a Hilbert basis (e n ) n∈N of V such that ∀n ∈ N, h n = U e n . Applying Theorem 1.1 of [START_REF] Xu | The expansion of a semigroup and a Riesz basis criterion[END_REF], it follows that there exists a sequence of root vectors of A which forms a Riesz basis of H. The following theorem gives an expansion of the solution of (P ).

Theorem 1 Assume that U 0 ∈ D(A). Let {λ k , k ∈ N} be the set of eigenvalues of A. For all k ∈ N, let r k be the geometric multiplicity of λ k , U k,i,j , i = 1, . . . r k , j = 0, . . . , m k,i -1 a familly of root vectors corresponding to λ k such that (A -λ k Id)U k,i,j = U k,i,j-1 , i = 1, . . . r k , j = 0, . . . , m k,i -1 (with U k,i,-1 = 0) and set U k,i,j = u k,i,j v k,i,j .

Then there exists a familly of polynomial functions P k,i,j of degree ≤ m k,i -1 -j, k ∈ N, i = 1, . . . r k , j = 0, . . . , m k,i -1 such that if for all n ∈ N and t ∈ R we set where a k,i,j (t) = e (Ω+λ k )t P k,i,j (t), then for all T > 0, u n → u in C 0 ([0, T ]; [H 2 (0, l)] 2 ) ∩ C 1 ([0, T ]; [H 1 (0, l)] 2 ) ∩ C 2 ([0, T ]; [L 2 (0, l)] 2 ).

|ν 2 n

 2 | → +∞ when |n| → +∞. If µ 1 = µ 2 the zeros of det M (λ) are asymptotically simple and inf n =m |λ n -λ m | > 0.

  i,j (t)u k,i,j


Ph.D Joseph Holung

Analysis of the magnetic field and torque in radially magnetized permanent magnet Couplings: Exact 2D formulation

A two-dimensional formulation for the magnetic field and torque in radially magnetized, permanent magnet couplings is given. Expressions for the field quantities are developed from the magnetic scalar potential in the air gap between the shaft and cylinder. The cogging torque in such multi-pole couplings is derived and compared with experimental results.

INTRODUCTION

Magnetic couplings are interesting because they can transmit torque without mechanical contact. A radial magnetic coupling consists of a shaft, fitted with a circular array of permanent magnet arc segments, which is separated by an air gap from a similar array of permanent magnet arc segments attached to the bore of a cylinder. The magnetization of each arc segment is assumed to be in a purely radial direction, either positive or negative, with the number and arrangement determining the number of poles in the coupling.

In the case of axial magnetic couplings, closed-form expressions for the force and torque have been derived from the magnetic vector potential obtained by solving Laplace's and Poisson's equation for a twodimensional (2D) analytical model [START_REF] Lubin | Simple Analytical Expressions for the Force and Torque of Axial Magnetic Couplings[END_REF]. The approach is well-suited to parametric studies often carried out in early engineering investigations to size geometry and lead to simpler expressions than other approaches such as a Coulombian formulation [START_REF] Ravaud | Torque in PM Couplings: Comparison of Uniform and Radial Magnetization[END_REF][START_REF] Mack | Electricity and Magnetism[END_REF][START_REF] Mack | Mathematics for the advanced Course in Engineering[END_REF].

In this paper, the magnetic field and torque in radially magnetized couplings is derived from a twodimensional analytical model. Exact closed-form expressions are obtained by solving Poisson's equation for the magnetic scalar potential in the permanent magnet regions and Laplace's equation in the air gap region.

PROBLEM DESRIPTION AND ASSUMPTIONS

A cross-section of an ideal, radially magnetized coupling is shown in Fig. 1. A circular array of permanent magnet arc segments extends radially to radius R2, from a shaft of radius R1. An air gap separates a corresponding array of permanent magnet arc segments attached to the bore of a cylinder with bore radius R4. The air gap radial clearance is R3-R2. The pole arc to pole pitch ratio is denoted by α. The alternating Fig 1 . Cross-section of a radially magnetized coupling. The pole arc to pole pitch ratio is α, the number of poles is p (p=6). The shaft and housing have infinite permeability (µ=∞, at r= R1 and R4). δ is an angular offset (torque angle) used as a parameter in the torque expression. colors denote the alternating positive and negative radial direction of the magnetization in each of the adjacent arc segments. With integer p denoting the number of pole pairs, Fig. 1 shows a p=6 pole-pair, radial, magnetic coupling. The cylinder is assumed to be sufficiently long and losses due to end effects and fringing negligible.

In order to simplify the two-dimensional analysis, the following additional assumptions are made:

1. The iron shaft and cylinder have infinite magnetic permeability, i.e., µ=∞ 2. The radially magnetized permanent magnets have relative recoil permeability µr =1. 1The analysis domain will be separated into three regions. Region 1 is the region occupied by permanent magnets attached to the shaft, R2 ≤ r ≤ R1, region 2 is the air gap region, R3 ≤ r ≤ R2, and region 3 is occupied by the outer constellation of permanent magnets, R4 ≤ r ≤R3. The magnets of region 3 are shown rotated by angle δ (torque angle) relative to the magnets in region 1. It is clear from the geometry that the magnetic field distribution is periodic with period T= 2π/p. The constitutive law relating the magnetic flux density vector, B, to the magnetic field intensity vector, H, and the magnetization vector, M, in the permanent magnet regions is

𝑩 = 𝜇( 𝑯 + 𝑴) ( 1
)
where µ =4πx10 -7 henrys/meter, is the permeability of space. As explained in [START_REF] Woodson | Electromechanical Dynamics[END_REF], this material model, Fig. 2. (Reproduced from [START_REF] Woodson | Electromechanical Dynamics[END_REF]). Hard magnetic material magnetized and then demagnetized by the current I.

shown in Fig. 2., represents "hard" magnetic materials with linear demagnetizing curve. The curve results when these materials are subjected to a field intensity H by means of an extremely large current as indicated by the inset. They explain that when the current (i.e., H) is removed, there is a residual flux density Br. Then, if the sample is removed from the circuit and subjected to other magnetic fields, the B-H curve settles down to operate at some point, such as "A" shown in the figure. The straight line indicated by the dotted line is the representative model of the behavior.

ANALYSIS

The magnetic field problem will be formulated in terms of the scalar magnetic potential. The field intensity is defined as the negative gradient of a magnetic scalar potential, ψ

𝑯 = -𝛁ψ (2)
Gauss's law of magnetism states that the magnetic field has divergence equal to zero,

𝛁 ⋅ 𝐁 = 𝟎 (3)
Substituting ( 2) into (1) and then imposing condition (3) results in the governing equations in each region. Thus, in the permanent magnet regions, 1 and 3, the governing equation is Poisson's equation, which, expressed in the two cylindrical coordinates (r, θ), is

𝛁 ⋅ 𝛁𝜓(𝑟, 𝜃) = 𝜕 2 𝜓 𝜕𝑟 2 + 1 𝑟 𝜕𝜓 𝜕𝑟 + 𝜕 2 𝜓 𝜕𝜃 2 = 𝛁 ⋅ 𝑴 (4)
In the air gap region, region 2, M=0, and there is no current source, so the governing equation is Laplace's equation,

𝜕 2 𝜓 𝜕𝑟 2 + 1 𝑟 𝜕𝜓 𝜕𝑟 + 𝜕 2 𝜓 𝜕𝜃 2 = 0 (5)
The magnetization vector to be considered has only a radial component, (M = Mr er), where er is a unit vector in the radial direction, and is distributed ideally as shown in Fig. 3, for region 3. The square wave Mr Fig. 3. Radial component of magnetization vector in region 3.

can be represented as a Fourier series and substituted into the right side of Eq. ( 4); the result is

𝜕 2 𝜓 𝜕𝑟 2 + 1 𝑟 𝜕𝜓 𝜕𝑟 + 𝜕 2 𝜓 𝜕𝜃 2 = 1 𝑟 𝑀 𝑟 (6)
𝑀 𝑟 = ∑ 𝑀 𝑘 * sin 𝑘𝑝 ( 𝜃 -𝛿) ∞ 𝑘=1,3,5,… (7)
where

𝑀 𝑘 = 4 𝐵 𝑟 𝑘 𝜋 µ * cos k 𝜋 2 (1 -𝛼) (8)
Mr consists of an infinite sum of sine waves of odd harmonics. The governing equation in region 1 is similar and is obtained by setting δ equal to zero. To complete the analytical model, certain boundary conditions must be satisfied. The tangential component of H must be zero at the soft iron boundaries (since µ= ∞), and the tangential component of H and the normal component of B must be continuous at the interfaces.

Using superscripts and subscripts to identify the regions and components, these boundary conditions are summarized in Table 1.

=0 r=R1 14

General solution of the field equations.

The general solution for the potential will be obtained by superposition. The analytic model will be split into two parts which, when solved and summed together, provide the solution. The geometry of the first excludes the magnets of region 1, and the second excludes the magnets of region 3. The two cases are shown in Fig. 4. In each case, a preliminary solution is first obtained by considering a single term in the driving function. The complete solution will be obtained by summation of additional terms as indicated by Eq. ( 7) and Eq. (8).

Solution for case (a) and kp≠1

Referring to Fig. 4(a), and considering a single term in the driving function, the governing equation in the permanent magnet outer ring is Poisson's equation,

𝜕 2 𝜓 𝜕𝑟 2 + 1 𝑟 𝜕𝜓 𝜕𝑟 + 𝜕 2 𝜓 𝜕𝜃 2 = 1 𝑟 𝑀 𝑘 𝑠𝑖𝑛 𝑘𝑝(𝜃 -𝛿) ( 15
)
The solution is

𝜓 3 (𝑟, 𝜃) = [𝑎 ( 𝑟 𝑅 4 ) 𝑘𝑝 + 𝑏 ( 𝑅 3 𝑟 ) 𝑘𝑝 + 𝑀 𝑘 * 𝑟 1 -(𝑘𝑝) 2 ] sin 𝑘𝑝(𝜃 -𝛿) ( 16
)
where a superscript on the scalar potential is used to denote the region, a and b are two arbitrary constants in the complementary solution, and the last term is the particular solution. This particular solution is valid for kp ≠1. The solution for kp=1 will be given later as a special case.

In the air region, (R3 ≤ r ≤ R1), the governing equation is Laplace's equation, Eq. ( 5), whose general solution is well-known. The solution must also be chosen to satisfy the boundary conditions, which for case (a) is given by Eqs. (9-11), and Eq. ( 14), in Table 1, adjusting the superscripts for the interface boundaries. Accordingly, the solution is of the form

𝜓 (𝑎) (𝑟, 𝜃) = [𝐴 ( 𝑟 𝑅 1 ) 𝑘𝑝 + 𝐵 ( 𝑅 3 𝑟 ) 𝑘𝑝 ] [𝐶 sin(𝑘𝑝 * 𝜃) + 𝐷 cos(𝑘𝑝 * 𝜃)] ( 17
)
where a superscript (a) is used in this case to represent the scalar potential in the air region for case (a), and A, B, C, and D, are arbitrary constants. All the constants in Eq. ( 16) and Eq. ( 17) can be determined by applying the boundary conditions. The details are provided in the Appendix. The result for the scalar potential in the permanent magnet region is

𝜓 3 (𝑟, 𝜃) = 𝐽(𝑟) * sin 𝑘𝑝(𝜃 -𝛿) ( 18
)
where

𝐽(𝑟) = 𝑏 [ ( 𝑅 3 𝑟 ) 𝑘𝑝 -( 𝑅 3 𝑅 4 ) 𝑘𝑝 ( 𝑟 𝑅 4 ) 𝑘𝑝 ] -𝑅 4 ( 𝑟 𝑅 4 ) 𝑘𝑝 * 𝑀 𝑘 1-(𝑘𝑝) 2 + 𝑀 𝑘 * 𝑟 1-(𝑘𝑝) 2 (19) and 𝑏 = 𝑀 𝑘 * 𝑅 3 2 [1-(𝑘𝑝) 2 ][ ( 𝑅 3 𝑅 1 ) 2𝑘𝑝 -( 𝑅 3 𝑅 4 ) 2𝑘𝑝 ] * [ 2 ( 𝑅 3 𝑅 4 ) 𝑘𝑝-1 -(1 + ( 𝑅 3 𝑅 1 ) 2𝑘𝑝 ) -(1 -( 𝑅 3 𝑅 1 ) 2𝑘𝑝 ) (𝑘𝑝)] (20)
The result for the scalar potential in the air gap is

𝜓 (𝑎) (𝑟, 𝜃) = 𝐽(𝑟 = 𝑅 3 ) (1 -( 𝑅 3 𝑅 1 ) 2𝑘𝑝 ) [ ( 𝑅 3 𝑟 ) 𝑘𝑝 -( 𝑅 3 𝑅 1 ) 𝑘𝑝 ( 𝑟 𝑅 1 ) 𝑘𝑝 ] 𝑠𝑖𝑛 𝑘𝑝(𝜃 -𝛿) (21)
The two components of the magnetic flux density vector in the air region can be computed from Eq. ( 21). The radial component is

𝐵 𝑟 (𝑎) (𝑟, 𝜃) = -µ 𝜕𝜓 (𝑎) 𝜕𝑟 = 𝜇 𝐽(𝑟 = 𝑅 3 ) * (𝑘𝑝) (1 -( 𝑅 3 𝑅 1 ) 2𝑘𝑝 ) * 𝑅 3 [ ( 𝑅 3 𝑟 ) 𝑘𝑝+1 + ( 𝑅 3 𝑅 1 ) 𝑘𝑝+1 ( 𝑟 𝑅 1 ) 𝑘𝑝-1 ] sin 𝑘𝑝(𝜃 -𝛿) (22)
and the tangential component is

𝐵 𝜃 (𝑎) (𝑟, 𝜃) = -µ 𝜕𝜓 (𝑎) 𝑟 𝜕𝜃 = -𝜇 𝐽(𝑟 = 𝑅 3 ) * (𝑘𝑝) (1 -( 𝑅 3 𝑅 1 ) 2𝑘𝑝 ) * 𝑟 [ ( 𝑅 3 𝑟 ) 𝑘𝑝 -( 𝑅 3 𝑅 1 ) 𝑘𝑝 ( 𝑟 𝑅 1 ) 𝑘𝑝 ] cos 𝑘𝑝(𝜃 -𝛿) (23)

Solution for case (b) and kp≠1

For case (b), set δ=0 in Eq. ( 15) to obtain the governing equation in the permanent magnet inner ring of region 1,

𝜕 2 𝜓 𝜕𝑟 2 + 1 𝑟 𝜕𝜓 𝜕𝑟 + 𝜕 2 𝜓 𝜕𝜃 2 = 1 𝑟 𝑀 𝑘 𝑠𝑖𝑛 (𝑘𝑝𝜃) ( 24
)
The solution is

𝜓 1 (𝑟, 𝜃) = [𝑎 ( 𝑟 𝑅 1 ) 𝑘𝑝 + 𝑏 ( 𝑅 2 𝑟 ) 𝑘𝑝 + 𝑀 𝑘 * 𝑟 1 -(𝑘𝑝) 2 ] sin(𝑘𝑝𝜃) (25)
where a and b are two arbitrary constants (different from case (a)). In the air region, (R2 ≤ r ≤ R4), the solution to Laplace's equation must be chosen to satisfy the boundary conditions, which for case (b) is given by Eqs. (12-14), and Eq. ( 9) in Table 1. The solution is

𝜓 (𝑏) (𝑟, 𝜃) = [𝐴 ( 𝑟 𝑅 4 ) 𝑘𝑝 + 𝐵 ( 𝑅 2 𝑟 ) 𝑘𝑝 ] [𝐶 sin(𝑘𝑝 * 𝜃) + 𝐷 cos(𝑘𝑝 * 𝜃)] ( 26
)
where the superscript (b) is used to represent the scalar potential in the air region for case (b), and A, B, C, and D, are arbitrary constants. All the constants in Eq. (25) and Eq. ( 26) can be determined by applying the boundary conditions as before. Considerable work can be avoided by noting similarities in the boundary conditions and equations compared with case (a). For example, we can replace R4 and R3 in the magnet region for case (a) with R1 and R2, to obtain the corresponding solution for case (b). Similarly, by replacing R3 and R1 in the air region for case (a) with R2 and R4, respectively, we can obtain the corresponding solution for case (b). With this in mind, we can write the solution for the scalar potential in the permanent magnet region

𝜓 1 (𝑟, 𝜃) = 𝐿(𝑟) * sin(𝑘𝑝𝜃) (27)
where

𝐿(𝑟) = 𝑏 [ ( 𝑅 2 𝑟 ) 𝑘𝑝 -( 𝑅 2 𝑅 1 ) 𝑘𝑝 ( 𝑟 𝑅 1 ) 𝑘𝑝 ] -𝑅 1 ( 𝑟 𝑅 1 ) 𝑘𝑝 * 𝑀 𝑘 1-(𝑘𝑝) 2 + 𝑀 𝑘 * 𝑟 1-(𝑘𝑝) 2 (28) and 𝑏 = 𝑀 𝑘 * 𝑅 2 2 [1-(𝑘𝑝) 2 ][ ( 𝑅 2 𝑅 4 ) 2𝑘𝑝 -( 𝑅 2 𝑅 1 ) 2𝑘𝑝 ] * [ 2 ( 𝑅 2 𝑅 1 ) 𝑘𝑝-1 -(1 + ( 𝑅 2 𝑅 4 ) 2𝑘𝑝 ) -(1 -( 𝑅 2 𝑅 4 ) 2𝑘𝑝 ) (𝑘𝑝)]
(29)

The scalar potential in the air gap is

𝜓 (𝑏) (𝑟, 𝜃) = 𝐿(𝑟 = 𝑅 2 ) (1 -( 𝑅 2 𝑅 4 ) 2𝑘𝑝 ) [ ( 𝑅 2 𝑟 ) 𝑘𝑝 -( 𝑅 2 𝑅 4 ) 𝑘𝑝 ( 𝑟 𝑅 4 ) 𝑘𝑝 ] 𝑠𝑖𝑛 (𝑘𝑝𝜃) (30)
For the magnetic flux vector in air, we have for the radial component

𝐵 𝑟 (𝑏) (𝑟, 𝜃) = -µ 𝜕𝜓 (𝑏) 𝜕𝑟 = 𝜇 𝐿(𝑟 = 𝑅 2 ) * (𝑘𝑝) (1 -

Scalar potential and torque, kp≠1

The solution for all field quantities is obtained by superposition. For example, the magnetic scalar potential in the air gap region of the coupling is

𝜓(𝑟, 𝜃) = ∑ (𝜓 (𝑎) + (𝜓 (𝑏) ) (33) ∞ 𝑘=1,3,5…
where Eqs. ( 21) and (30) can be substituted in the above, and as many harmonic terms of the forcing function can be included as desired, in accordance with Eqs. ( 7) and (8).

Of special interest is the torque in the coupling; this can be computed from the Maxwell stress tensor, which, in electromagnetism with only a magnetic field, takes the form

𝜎 𝑖𝑗 = 1 µ 𝐵 𝑖 𝐵 𝑗 - 1 2µ 𝐵 2 𝛿 𝑖𝑗 ( 34
)
where Bi, Bj are components of the flux density vector, and δ ij is the Kronecker delta. The magnetic torque is obtained by considering a circular integration path of any radius, R, within the air gap around the shaft and the tangential component of stress σ rθ from Eq. ( 34). The differential torque, dT, acting on an element of area Rdθdz is dT= R *σ rθ * (Rdθdz). Thus, the total magnetic torque per unit length acting on the cylindrical coupling is

𝑇 𝐿 ⁄ = 𝑅 2 µ ∫ 𝐵 𝑟 (𝑅, 𝜃) * 𝐵 𝜃 (𝑅, 𝜃)𝑑𝜃 (35) 𝜋 -𝜋
where L is the axial length of the magnet segments. By adding all the contributions to the field in the same manner as was done for the scalar potential, we now express the torque per unit length as An expression for the torque is obtained by first substituting for each term using Eqs. ( 22), ( 23), (31), and (32), and carrying out the integration. The result is the first and fourth integral is zero and only the two middle terms contribute to the torque; therefore, it is the interaction of the two fields that give rise to the torque. The final resulting expression for the torque in the coupling can be written compactly as

𝑇 𝐿 ⁄ = 𝑅 2 µ ∑ ∫ [
𝑇 𝐿 ⁄ = 𝜋𝑅 2 µ ∑ ( 𝐵 𝑟0 (𝑏) * 𝐵 𝜃0 (𝑎) -𝐵 𝑟0 (𝑎) * 𝐵 𝜃0 (𝑏) ) sin 𝛿 ∞ 𝑘 ( 38
)
where the extra subscript 0 indicates the amplitude of each respective sinusoid term, which we evaluated before; namely, 𝐵 𝑟0 (𝑎) = 𝜇

Special case: kp=1, case (a)

The results obtained in the preceding sections must be modified when kp=1. Poisson's equation for case (a) now becomes

𝜕 2 𝜓 𝜕𝑟 2 + 1 𝑟 𝜕𝜓 𝜕𝑟 + 𝜕 2 𝜓 𝜕𝜃 2 = 1 𝑟 𝑀 𝑘 𝑠𝑖𝑛(𝜃 -𝛿) ( 43
)
The solution is

𝜓 3 (𝑟, 𝜃) = [𝑎 ( 𝑟 𝑅 4 ) + 𝑏 ( 𝑅 3 𝑟 ) + 𝑀 𝑘 * ln 𝑟 * 𝑟 2 ] sin(𝜃 -𝛿) (44)
Thus, it is seen that the new solution is obtained from the previous by setting kp=1, and replacing the term Mk/(1-(kp) 2 ) in the particular solution with Mk/2* ln r.

Similarly, in the air region the solution to Laplace's equation is now written as ))] (48)

The result for the scalar potential in the air gap is

𝜓 (𝑎) (𝑟, 𝜃) = 𝐽(𝑟 = 𝑅 3 ) (1 -( 𝑅 3 𝑅 1 ) 2 ) [ ( 𝑅 3 𝑟 ) -( 𝑅 3 𝑅 1 ) ( 𝑟 𝑅 1 ) ] 𝑠𝑖𝑛 (𝜃 -𝛿) (49)
The two components of the magnetic flux density vector in the air region can be computed from Eq. ( 49). The radial component is

𝐵 𝑟 (𝑎) (𝑟, 𝜃) = -µ 𝜕𝜓 (𝑎) 𝜕𝑟 = 𝜇 𝐽(𝑟 = 𝑅 3 ) (1 -( 𝑅 3 𝑅 1 ) 2 ) * 𝑅 3 [ ( 𝑅 3 𝑟 ) 2 + ( 𝑅 3 𝑅 1 ) 2 ] sin(𝜃 -𝛿) (50)
and the tangential component is

𝐵 𝜃 (𝑎) (𝑟, 𝜃) = -µ 𝜕𝜓 (𝑎) 𝑟 𝜕𝜃 = -𝜇 𝐽(𝑟 = 𝑅 3 ) (1 -( 𝑅 3 𝑅 1 ) 2 ) * 𝑟 [ ( 𝑅 3 𝑟 ) -( 𝑅 3 𝑅 1 ) ( 𝑟 𝑅 1 ) ] cos(𝜃 -𝛿) (51)
Special case: kp=1, case (b)

Poisson's equation for case (b) is

𝜕 2 𝜓 𝜕𝑟 2 + 1 𝑟 𝜕𝜓 𝜕𝑟 + 𝜕 2 𝜓 𝜕𝜃 2 = 1 𝑟 𝑀 𝑘 𝑠𝑖𝑛𝜃 (52)
Obtained by letting δ=0. The solution is

𝜓 1 (𝑟, 𝜃) = [𝑎 ( 𝑟 𝑅 1 ) + 𝑏 ( 𝑅 2 𝑟 ) + 𝑀 𝑘 * ln 𝑟 * 𝑟 2 ] sin 𝜃 (53)
Similarly, in the air region the solution to Laplace's equation is now written as

𝜓 (𝑏) (𝑟, 𝜃) = [𝐴 ( 𝑟 𝑅 4 ) + 𝐵 ( 𝑅 2 𝑟 ) ] [𝐶 sin 𝜃 + 𝐷 cos 𝜃] (54)
All the constants are determined by applying the boundary conditions as before. The result for the scalar potential in the permanent magnet region is

𝜓 1 (𝑟, 𝜃) = 𝐿(𝑟) * sin 𝜃 (55)
where

𝐿(𝑟) = 𝑏 [ ( 𝑅 2 𝑟 ) -( 𝑅 2 𝑅 1 ) ( 𝑟 𝑅 1 ) ] -𝑅 1 ( 𝑟 𝑅 1 ) * 𝑀 𝑘 2 * ln 𝑅 1 + 𝑀 𝑘 * 𝑟 2 ln 𝑟 (56)
and

𝑏 = 𝑀 𝑘 * 𝑅 2 2 [2][ ( 𝑅 2 𝑅 4 ) 2 -( 𝑅 2 𝑅 1 ) 2 ] * [(1 -( 𝑅 2 𝑅 4 ) 2 -2 * ln( 𝑅 2 𝑅 1 ))] (57)
The result for the scalar potential in the air gap is

𝜓 (𝑏) (𝑟, 𝜃) = 𝐿(𝑟 = 𝑅 2 ) (1 -( 𝑅 2 𝑅 4 ) 2 ) [ ( 𝑅 2 𝑟 ) -( 𝑅 2 𝑅 4 ) ( 𝑟 𝑅 4 ) ] 𝑠𝑖𝑛 𝜃 (58)
The two components of the magnetic flux density vector in the air region can be computed from Eq. ( 58). The radial component is

𝐵 𝑟 (𝑏) (𝑟, 𝜃) = -µ 𝜕𝜓 (𝑎) 𝜕𝑟 = 𝜇 𝐿(𝑟 = 𝑅 2 ) (1 -( 𝑅 2 𝑅 4 ) 2 ) * 𝑅 2 [ ( 𝑅 2 𝑟 ) 2 + ( 𝑅 2 𝑅 4 ) 2 ] sin 𝜃 ( 59
)
and the tangential component is

𝐵 𝜃 (𝑎) (𝑟, 𝜃) = -µ 𝜕𝜓 (𝑎) 𝑟 𝜕𝜃 = -𝜇 𝐿(𝑟 = 𝑅 2 ) (1 -( 𝑅 2 𝑅 4 ) 2 ) * 𝑟 [ ( 𝑅 2 𝑟 ) -( 𝑅 2 𝑅 4 ) ( 𝑟 𝑅 4 ) ] cos 𝜃 (60)

RESULTS

To illustrate the utility of the derived closed-form expressions, computations were done on a radially magnetized, multi-pole coupling, with parameters given in Table 2.

The torque calculations were carried out by taking a contour in the air gap with radius equal to the mean radius, (i.e., R= (R2+R3)/2). Table 3 shows the amplitude of the cogging torque per unit length, computed from the torque expressions, for different number of pole-pairs as well as the contribution of the first four terms of k included in the summation. The sum converges rapidly and even a single term (k=1) approximation would yield a sufficiently accurate estimate for most preliminary sizing study. Within the range examined for this coupling, the maximum cogging torque that could be achieved is 29.37 N-m per meter of axial length of the magnet segments. This corresponds to a 3 pole-pair design, as shown in Fig. 5.

Torque measurements were taken on several couplings prepared by BJA Magnetics, [7]. The magnet segments were constructed from N52 (NdFeB) magnets and all samples correspond to the 6 pole-pair geometry of Table 2. The steel shaft and housing were made from SUS403. The magnet arc segments were magnetized diametrically (i.e., in a uniform direction) rather than in a truly radial direction. As others have noted [START_REF] Ravaud | Torque in PM Couplings: Comparison of Uniform and Radial Magnetization[END_REF], this is due to difficulty in producing a radial magnetization field. As the number of arc segments increases, the field will better approximate a true radial field.

The torque measurements are shown in Fig. 6. Measurements were made in both clockwise and counterclockwise directions and cover one full revolution of the shaft. The average amplitude was approximately 70 Kgf-mm. From Table 3, the predicted amplitude of the cogging torque for an ideal, radially magnetized, 6 pole-pair coupling with this geometry would be 95.9 Kgf-mm. It is interesting to note that the prediction could be adjusted and brought into better correlation with the experimental measurements by replacing the value of Br with a reduced value B0, (Refer to Fig. 2 and the discussion of it). 2.

Fig. 6. Torque measurements on several samples of the 6 pole-pair, radially magnetized coupling of Table 2.

CONCLUSION

The scalar potential for multi-pole, radially magnetized couplings have been derived from the solution of Poisson's equation in the permanent magnet regions, and Laplace's equation in the air-gap region. All field quantities of interest can be obtained from the scalar potentials. In particular, the closed-form expression for the cogging torque was developed, which shows the sinusoidal behavior. The theoretical prediction for a 6 pole-pair coupling was carried out, compared with measurements taken on several samples, and the utility of the closed-form expression illustrated with a simple parametric study. Finally, it was explained how the model can be correlated to experimental findings in practice. Substitute Eq. (A2) into Eq. ( 17)

𝜓 (𝑎) (𝑟, 𝜃) = [ ( 𝑅 3 𝑟 ) 𝑘𝑝 -( 𝑅 3 𝑅 1 ) 𝑘𝑝 ( 𝑟 𝑅

Fig. 4 .

 4 Fig. 4. The analytic model split into two parts: (a) the model geometry without the inner magnets of region 1, and (b), the model geometry without the outer magnets of region 3. The iron shaft has radius r= R1, and the outer boundary is at r= R4 for both.

Fig. 5 .

 5 Fig. 5. Torque amplitude computed from closed-form expressions for the coupling of Table2.

(

  sin(𝑘𝑝 * 𝜃) + 𝐵𝐷 cos(𝑘𝑝 * 𝜃)] (𝐴3) Apply boundary condition, Eq. (10), cos(𝑘𝑝 * 𝜃) -𝐵𝐷 sin(𝑘𝑝 * 𝜃)] (𝐴4) Equating coefficients of the sine and cosine terms, we obtain the two relations 𝐽(𝑟 = 𝑅 3 ) cos 𝑘𝑝𝛿 = [We can eliminate δ by squaring both sides of Eqs. (A5) and (A6) and adding to get [𝐽(𝑟 = 𝑅 3 Apply boundary condition, Eq. (11), 𝐵 𝑟 3 =𝐵 𝑟 𝑎 or (-𝜕𝜓3(𝑅3) 𝜕𝑟 + 𝑀 𝑟 = -sin(𝑘𝑝 * 𝜃) + 𝐵𝐷 cos(𝑘𝑝 * 𝜃)] (𝐴8) Equating coefficients of the sine and cosine terms, we obtain the two relations We can now eliminate [BC 2 + BD 2 ] using Eqs. (A7) and (A11) to get 𝐽(𝑟 = 𝑅 3 ) By substituting Eq. (19) into Eq. (A12), the equation can be solved to determine the constant b. The result is Eq. (20).With b determined, we can compute J(r=R3) from Eq. (19). Then, for every value of δ, we can compute BC and BD from Eqs. (A5) and (A6). Then we have ψ, Br and Bθ. For example, substituting for BC and BD from Eqs. (A5) and A(6) into Eq. (A3), we arrive at Eq. (21):

Table 1 .

 1 Boundary Conditions at Interface

	Domain1	Domain2	Relation		Interface	Equation
	Iron	Magnet	𝐻 𝜃 3 =-1 𝑟	𝜕𝜓3 𝜕𝜃	=0	r=R4	9
	Magnet	Air gap	𝐻 𝜃 3 =𝐻 𝜃 2			r=R3	10
			(-1 𝑟	𝜕𝜓3 𝜕𝜃	= -1 𝑟	𝜕𝜓2 𝜕𝜃	)
	Magnet	Air gap	𝐵 𝑟 3 =𝐵 𝑟 2			r=R3	11
	Air gap	Magnet	𝐻 𝜃 2 =𝐻 𝜃 1			r=R2	12
			(-1 𝑟	𝜕𝜓2 𝜕𝜃	= -1 𝑟	𝜕𝜓1 𝜕𝜃	)
	Air gap	Magnet	𝐵 𝑟 2 =𝐵 𝑟 1			r=R2	13
	Magnet	Iron	𝐻 𝜃 1 =-1 𝑟	𝜕𝜓1 𝜕𝜃

  All the constants are determined by applying the boundary conditions as before. The result for the scalar potential in the permanent magnet region is

	𝜓 (𝑎) (𝑟, 𝜃) = [𝐴 ( 𝑅 1 𝑟	) + 𝐵 ( 𝑅 3 𝑟	) ] [𝐶 sin 𝜃 + 𝐷 cos 𝜃]	(45)
	𝜓 3 (𝑟, 𝜃) = 𝐽(𝑟) * sin(𝜃 -𝛿)			(46)
	where
	𝐽(𝑟) = 𝑏 [ ( 𝑅 3 𝑟 ) -( 𝑅 3 𝑅 4	) ( 𝑅 4 𝑟	) ] -𝑅 4 ( 𝑅 4 𝑟	) *	𝑀 𝑘 2 * ln 𝑅 4 +	𝑀 𝑘 * 𝑟 2 ln 𝑟	(47)
	and
	𝑏 =	𝑀 𝑘 * 𝑅 3 2 [2][ ( 𝑅 3 𝑅 1 ) 2 -(	𝑅 3 𝑅 4	) 2	]	* [(1 -( 𝑅 3 𝑅 1	) 2	-2 * ln(	𝑅 3 𝑅 4

Table 2 .

 2 Parameters of radial coupling studied

	Item	Description	Value
	R1	Inner radius of shaft magnets	2.5 x 10 -3 m
	R2	Outer radius of shaft magnets	3.0 x10 -3 m
	R3	Inner radius of cylinder magnets	3.25 x10 -3 m
	R4	Outer radius of cylinder magnets	5 x10 -3 m
	L	Axial length of magnets	37 x 10 -3 m
	Br	Residual flux density, N52 magnet	1.44 Tesla
	α	Pole arc to pole pitch ratio	0.9
	p	Number of pole-pairs	variable

Table 3 .

 3 Computed torque amplitude per unit length, N-m/m, on coupling of Table2.

				Number of pole-pairs, p
	Index, k 1	2	3	4	5	6
	1	15.95	24.60	27.33	27.27	26.08	24.48
	3	1.76	2.21	1.76	1.39	1.09	0.86
	5	0.11	0.37	0.25	0.17	0.11	0.07
	7	0.07	0.06	0.03	0.02	0.01	0.01
	∑(T/L) = 17.89	27.24	29.37	28.85	27.29	25.42

This assumption is reasonable for the neodymium-iron-boron (NdFeB) magnets which have µr =1.05.

ACKNOWLEDGEMENT

The author gratefully acknowledges the funding support from Kenneth Seethaler, Manager, Strategic Technology and Innovation Center, Lenovo Corporation, Morrisville, NC, USA. Thanks also to Richard Helle and Daniel Bliss, BJA Magnetics, for many hours spent on this project, and for preparing all the samples for evaluation, and to Dr. Ryan Murray, Department of Mathematics, NC State University, for his helpful advice.

APPENDIX

All the constants can be determined from the boundary conditions and will be illustrated for case (a), kp ≠ 1. First, apply boundary condition, Eq. ( 9), (


J Bischoff

C Delafoy

N Chaari

C Vauglin

K Buchanan

P Barberis

F Schuster

J-C Brachet

K Nimishakavi

CR-COATED CLADDING DEVELOPMENT AT FRAMATOME

Framatome is actively developing an evolutionary enhanced Accident Tolerant Fuel (E-ATF) design showing enhanced performance in both nominal and accidental conditions: Cr-coated M5 cladding combined with Cr2O3-doped UO2 fuel. The Cr-coated M5 cladding consists of a thin, dense and extremely adherent chromium layer deposited on the external surface of the M5 cladding. The main advantage of this solution is that it significantly reduces the high temperature steam oxidation, and thus limits the heat and hydrogen production usually observed in the oxidation reaction of steam with zirconium. Furthermore, recent out-of-pile results have shown improved behavior of the Cr-coated cladding with respect to debris fretting, which suggests that this solution will enhance the reliability of fuel. Additionally, Framatome has developed a large scale Cr-coating deposition to demonstrate the feasibility of PVD deposition on full length tubes and to use this facility to produce the Cr-coated cladding tubes for Lead Fuel Rod (LFR) irradiations such as the one planned for insertion in 2019.

Introduction

The Fukushima-Daiichi accident triggered a worldwide research and development effort in the nuclear fuel industry to increase fuel margins in severe accidents leading to the development of Enhanced Accident Tolerant Fuels (E-ATF). In this context Framatome is actively developing several concepts ranging from short-term evolutionary (Cr-coated zirconium alloy cladding and Cr2O3-doped UO2 fuel) to long term revolutionary (SiC/SiC composite cladding) solutions, relying on its worldwide teams and partnerships, with programs and irradiations planned both in Europe and the United States. The short-term solution is defined as evolutionary since it very similar to the current Zr-based cladding/UO2 pellet fuel rods, and thus corresponds to an improvement or evolution from that basis. SiC/SiC composite cladding on the other hand corresponds to a significant shift since we would move from a metallic to a ceramic cladding and all the current regulations may have to change. Concerning the evolutionary solution, the advantage of combining the Cr-coated cladding with the Cr2O3-doped UO2 fuel is that both solutions provide complementary benefits, such as for example reduced fission gas release from the Cr2O3-doped UO2 pellet in transients and accidents, thus reducing internal rod pressure, which when combined to the reduced ballooning of the Cr-coated cladding provides significant improvements in maintaining the coolable geometry of the fuel rod in accident conditions. The benefits of this short-term evolutionary solution are described in more detail in another paper [START_REF] Delafoy | Benefits of Framatome's E-ATF evolutionary solution: Cr-coated cladding with Cr2O3-doped UO2 fuel[END_REF].

This article focuses on the developments of the Cr-coated M5 cladding, which is the most mature cladding solution developed by Framatome. It consists of a relatively thin Cr-coating (10-20µm) deposited on the surface of the M5 cladding using a Physical Vapor Deposition (PVD) technique.

The coating thus deposited is very dense and adherent therefore making it very protective. Additionally, the PVD technique does not modify the underlying substrate microstructure since no heat treatment is applied on the tubes during deposition and the increase in temperature due to the incident Cr atoms is relatively small. Figure 1

2.

Motivation: improved accident and normal operation behavior

High temperature steam oxidation resistance

The main goal of enhanced accident tolerant fuel cladding is to improve the behavior of the fuel in accidental conditions, which has always been an important part of fuel development, with in this case a particular focus on reducing the high temperature (HT) steam oxidation of the current zirconium alloy cladding. This focus significantly increased since the Fukushima accident where the rapid zirconium alloy oxidation led to the production of H2 which then resulted in the reactor building explosion. The Cr-coated cladding was developed as a solution to reduce the high temperature steam oxidation of the cladding and shows significant benefits in HT steam oxidation resistance as demonstrated by current results obtained for long exposure times at the Framatome Research Center in Paimboeuf. Figure 2 shows the compared oxidation kinetics at 1100°C of Crcoated and uncoated cladding up to 18000s and the visual aspect of the samples after quench. cladding samples and the visual aspect after quench In this figure the uncoated sample fails after quench after exposure to 1100°C steam for 12000s while the Cr-coated sample maintains very low oxidation kinetic (reduced by at least 15x compared to uncoated samples) up to at least 18000s with no failure after quench. The main reason for the difference in behavior is that the Cr-coating protects the cladding from oxidation and from the diffusion of oxygen within the underlying zirconium substrate which leads to cladding embrittlement and therefore rupture after quench. By reducing the HT steam oxidation significantly, the Cr-coated cladding also maintains post-quench ductility for longer times thus providing significant benefits in accidental conditions. This test confirms previous results obtained at the CEA at other temperatures and times [START_REF] Brachet | On-going Studies ar CEA on Chromium Coated Zirconium Based Nuclear Fuel Claddings for Enhanced Accident Tolerant LWR Fuel[END_REF][START_REF] Brachet | Behavior under LOCA Conditions of Enhanced Accident Tolerant Chromium Coated Zircaloy-4 Claddings[END_REF][START_REF] Brachet | Behavior of Chromium Coated M5 Claddings under LOCA Conditions[END_REF].

Debris wear resistance in normal operation

Although the Cr-coated cladding was developed initially to improve performance in accident conditions, it also brings significant benefits in normal operating conditions [START_REF] Bischoff | AREVA NP's Enhanced Accident Tolerant Fuel Developments: Focus on Cr-coated M5 Cladding[END_REF]. For example, since Cr is harder than Zr, it will protect the cladding for debris wear. To demonstrate this a parametric test was developed at the Framatome Technical Center of Le Creusot to simulate an extremely harsh debris wear in order to compare the behavior of Cr-coated and uncoated cladding. The debris used was a stainless steel wire that was fixed and a sliding displacement of +/-200µm was imposed on the cladding with a contact force on the wire of 1.5N. The set-up was placed in an autoclave at 300°C with a loop controlling the water chemistry to be representative of PWR chemistry ([B]=1000ppm and [Li]=2.2ppm). The goal of this test was not to be fully representative of conditions in reactor but was meant to be extremely harsh so as to highlight the difference in behavior of the two types of cladding. Similar tests to simulate more precisely grid-to-rod fretting had already been performed and are presented other articles [START_REF] Delafoy | Benefits of Framatome's E-ATF evolutionary solution: Cr-coated cladding with Cr2O3-doped UO2 fuel[END_REF][START_REF] Bischoff | AREVA NP's Enhanced Accident Tolerant Fuel Developments: Focus on Cr-coated M5 Cladding[END_REF].

For the uncoated sample, the test had to be stopped after only 23h due to significant wear while the Cr-coated sample was tested for 100h. After 23h, the maximum wear depth on cladding was measured around 500µm corresponding to almost the entire cladding thickness, while for the Crcoated sample after 100h the maximum wear depth was around 280µm so almost half of the uncoated sample. Figure 3 shows a schematic of the test set-up and the visual aspect of the Crcoated and uncoated tubes after the tests. In the visual aspects after the tests it can be clearly seen that the notch on the Cr-coated sample is reduced in size and depth compared to the uncoated sample. It is also interesting to notice that even though the coating has been penetrated fully the wear rate is still low, which may be explained by the fact that the steel wire still has contact with the Cr-coating on the edges of the worn area.

The main cause of fuel rod failures in reactors is fretting wear, whether it be debris or grid-to-rod fretting wear. This result combined with previous results on grid-to-rod fretting wear resistance [START_REF] Bischoff | Development of Cr-coated Zirconium Cladding for Enhanced Accident Tolerance[END_REF] demonstrates that the Cr-coated can reduce cladding wear significantly and therefore its implementation will most likely reduce the number of failed fuel rods. This is an important benefit for the utilities since a Cr-coated fuel rod will improve the reliability of the fuel.

3.

Manufacturing developments

All the out-of-pile characterizations were performed using samples fabricated in a PVD machine that limited the size of the samples to 50cm in length. Foreseeing the need of insertion of LFRs in a very short time frame to obtain in-pile results as soon as possible, Framatome launched the fabrication of a full length prototype in 2016. The fabrication of the prototype to coat full-length (4m long) tubes was performed with two specific goals in mind:

-Demonstrate the feasibility of full length Cr-coated cladding using PVD -Support the insertion of full-length Cr-coated lead fuel rods (LFRs) in commercial PWRs.

The prototype machine was finished being built in summer of 2017 with the first Cr-coating tubes produced in the fall 2017. Several months were spent testing the machine after start-up and performing an optimization of the parameters to obtain the same quality of coatings as on the smaller length, and to make sure the coating was uniform in thickness and quality over the full length. Figure 4 shows an example of full-length Cr-coated tubes fabricated using the prototype. The tubes have a uniform visual aspect along the full length suggesting a coating with a homogeneous quality over the full length. This visual uniformity is also confirmed through coating thickness measurements using the calotest method (often used for coatings) during trial tests given in Figure 5, which show similar thicknesses within the uncertainty of the measurement (+/-1.5µm) over the whole length of the tube. The calotest method corresponds to the grinding of the coating surface up to the substrate with a ball of known diameter, and the radius of the various "circles" obtained by the grinding, corresponding to the edges of the coating, allow the operator to calculate back the thickness of the coating.

Additionally, HT air oxidation tests at 1100°C for 850s were performed on segments of the tubes at various axial positions to make sure the behavior of the coating was the same along the length. Figure 6 shows the visual aspect of the samples oxidized in these conditions and immediately quenched in room temperature water. Both samples have a similar visual appearance and the weight gains were close with 4.97 and 4.52 mg/cm 2 for extremity and center positions, respectively. The oxidation was performed on open tube samples, therefore double-sided oxidation, with the Cr-coating only on the external surface, so the majority of the weight gain measured comes from the oxidation of the unprotected internal zirconium alloy surface. These weight gains are also similar to ones obtained for Cr-coated samples obtained previously using a smaller PVD chamber, thus confirming the reproducibility on the full length scale of the coating quality previously obtained on small length. Furthermore, the immediate quench in room temperature water is very harsh for the samples and highlights the excellent adherence of the coating and the similarity in behavior along the length.

Consequently, the fabrication of Cr-coated cladding on full length is feasible and shows similar properties as for samples coated previously in a smaller PVD chamber, thus confirming Framatome's ability to reproduce the same coating quality with the change of scale. Framatome is therefore currently preparing for the qualification of the full length Cr-coated tubes before the fabrication of the Cr-coated tubes that will be used for the LFRs in 2019. Appearance of Cr-coated tube samples after double-sided HT oxidation in air at 1100°C for 850s followed by a quench at room temperature (4.97 (extremity) and 4.52 (center) mg/cm 2 in weight gain).

4.

Obtaining in-pile data for future licensing

With the goal of accelerating the implementation of Cr-coated cladding in commercial reactors, Framatome has started irradiations of this solution early-on to obtain the necessary in-pile data to support future licensing needs. Framatome was therefore the first to irradiate E-ATF solutions in a commercial reactor with the insertion of Material Test Rods in the Gösgen reactor in 2016 as part of the IMAGO irradiation (Irradiation of Materials for Accident tolerant fuels in GOesgen). After 1 cycle the samples were inspected visually and showed good behavior with no delamination or deterioration of the Cr-coating, as presented more in detail in another article [START_REF] Girardin | Inspection Capabilities and In-pile Experience with Innovative and Enhanced Accident Tolerant Fuel Materials at KKG[END_REF]. Some samples were extracted after 1 cycle to be analyzed in hot cells but most of the samples were reinserted to continue the irradiation.

Additional irradiations of Cr-coated cladding have started or are planned in research reactors to evaluate the behavior of a fuel rod and the interaction of fuel pellets with the Cr-coated cladding:

-As part of the French Joint Research program with CEA and EDF, the irradiation of UO2-Cr-coated cladding rodlets (containing both Zy4 and M5 substrates) in the HALDEN reactor started in summer 2017. Several pellet diameters are tested to investigate the impact of pellet-cladding interaction on the Cr-coated cladding behavior. -As part of phase 2 of the DOE ATF program, an irradiation of Framatome's short term E-ATF concept (Cr2O3-doped UO2 fuel with Cr-coated M5 cladding) is planned to start in 2018 in the Advanced Test Reactor (ATR) at Idaho National Laboratory in the US as part of the ATF-2 irradiation. This will be the first irradiation of the Framatome's complete E-ATF fuel rod concept with chromia doped fuel and Cr-coated cladding.

These irradiations will provide complementary information to the IMAGO irradiation for the justification of future LFR irradiations with data more representative of the fuel rod by investigating the overall pellet-cladding interaction. Furthermore, since Framatome had launched several irradiations in various reactors, and especially commercial reactors, the closure of HALDEN will not impact significantly the Framatome Cr-coating irradiation program and the results that would have been obtained in HALDEN can be obtained through other irradiations. Thus the closure of HALDEN will have minimal impact of Framatome's irradiation feedback and will not delay the implementation of its evolutionary solution. The rodlets irradiated in HALDEN will nonetheless be investigated and will be useful to obtain irradiated data on Cr-coated cladding at an intermediate dose.

5.

Conclusion

Framatome is strongly engaged in the E-ATF developments relying on its worldwide teams with programs both in Europe and the USA. Its objective is to accelerate the implementation of its near-term solution consisting of Cr2O3-doped UO2 fuel combined with Cr-coated M5 cladding, which provides significant benefits in both normal and accidental conditions. The main focus of the current developments is on the Cr-coated cladding since the doped fuel has already been industrialized. For example, the Cr-coating increases the wear resistance of the cladding suggesting improved reliability of the fuel in normal operation, as well as reduced steam oxidation resistance at high temperature increasing safety margins in design basis accidents or providing additional coping time in severe accidents.

Due to the good behavior of the Cr-coated cladding in out-of-pile tests, Framatome has fabricated a PVD prototype machine in partnership with the CEA to be able to coat full length cladding tubes.

The goal was on one hand to demonstrate the feasibility of PVD coatings on the full length and to develop a facility to be able to fabricate full length Cr-coated tubes. The construction was finalized in summer 2017 and Framatome has since been performing tests and optimization of the prototype to obtain the same quality of coatings as previously observed on short length with good axial uniformity over the full length. Framatome is now gearing up to qualify and fabricate the full length Cr-coated cladding tubes that will serve for the LFRs planned in 2019 in the USA.

Additionally irradiations were started or are going to be launched in commercial and material test reactors to obtain the necessary in-pile data to support future licensing of the fuel rod. Framatome was for example the first to irradiate E-ATF concepts in a commercial reactor with the irradiation of test samples in the Gösgen reactor (IMAGO). These data will be essential to implement Framatome's near-term solution in the upcoming years.

Fig 1 .

 1 Fig 1. Description of the Cr-coated M5 cladding

Fig 2 .

 2 Fig 2. Comparison of the HT steam oxidation kinetics at 1100°C of Cr-coated and uncoated M5cladding samples and the visual aspect after quench In this figure the uncoated sample fails after quench after exposure to 1100°C steam for 12000s while the Cr-coated sample maintains very low oxidation kinetic (reduced by at least 15x compared to uncoated samples) up to at least 18000s with no failure after quench. The main reason for the difference in behavior is that the Cr-coating protects the cladding from oxidation and from the diffusion of oxygen within the underlying zirconium substrate which leads to cladding embrittlement and therefore rupture after quench. By reducing the HT steam oxidation significantly, the Cr-coated cladding also maintains post-quench ductility for longer times thus providing significant benefits in accidental conditions. This test confirms previous results obtained at the CEA at other temperatures and times[START_REF] Brachet | On-going Studies ar CEA on Chromium Coated Zirconium Based Nuclear Fuel Claddings for Enhanced Accident Tolerant LWR Fuel[END_REF][START_REF] Brachet | Behavior under LOCA Conditions of Enhanced Accident Tolerant Chromium Coated Zircaloy-4 Claddings[END_REF][START_REF] Brachet | Behavior of Chromium Coated M5 Claddings under LOCA Conditions[END_REF].

Fig 3 .

 3 Fig 3. Schematic of the debris wear set-up and visual aspect of Cr-coated and uncoated cladding after the test.

Fig 4 .

 4 Fig 4. Example of Cr-coated full length tubes fabricated using the prototype showing uniformity of the visual aspect over the full length.

Fig 5 .

 5 Fig 5. Cr-coating thicknesses along the length of tubes showing thickness uniformity.

Acknowledgements

The authors would especially like to thank Eric Monsifrot for all his help and expertise on the fabrication of coatings and also Kernkraftwerk Gösgen-Däniken AG for their collaboration and financial contribution in the irradiation of E-ATF samples in the Gösgen reactor, and their support for the licensing of Material Test Rods and future examination of the samples. Some of the data presented in this article is based upon work supported by the U.S. Department of Energy under Award Number DE-NE0008220 with Framatome Inc.

M5 is a trademark or a registered trademark of Framatome or its affiliates in the USA or other countries.