Frédéric Mallet
email: fmallet@sophia.inria.fr

Robert De Simone

Laurent Rioux
email: laurent.rioux@thalesgroup.com

Event-triggered vs. Time-Triggered communications with UML MARTE

In the real-time and embedded domain, systems tend to combine periodic and aperiodic computations. This leads to mixing event-triggered with timetriggered communications with their pros and cons. Then, modeling standards of the domain must provide mechanisms to support both kinds whereas historically they pertain to different communities: asynchronous and synchronous designers. In this paper, we compare the expressiveness of two standards of the domain (AADL and MARTE) to model these two kinds of communications. Specifically we focus on the Time facilities of MARTE and on AADL models amenable to end-to-end flow latency analyses.

Introduction

Embedded applications often combine aperiodic (or sporadic) and periodic computations. In the automotive industry, this has lead to mixing eventtriggered communications (for aperiodic computations) with time-triggered communications (for periodic computations) in bus standards like FlexRay (http://www.flexray.org) or TT-CAN [START_REF] Albert | Comparison of Event-Triggered and Time-Triggered Concepts with Regard to Distributed Control Systems[END_REF]. In the avionic industry, generally an application mixes aperiodic events (e.g., interactions with the pilot generate aperiodic events, plane modes (air or ground)) and periodic events coming with system updating (e.g., fuel quantity, update system data…). Time-triggered approaches enhance predictability by reducing latency jitters and provide higher dependability by making it easier to detect missed messages or illegal accesses to the bus. However, event-triggered systems are more flexible to support configuration changes without a complete redesign and adapt faster to asynchronous events. In Electronic Design Automation (EDA) eventdriven simulators (like those for VHDL or Verilog) provide a large flexibility and support the design of both synchronous and asynchronous architectures. Though, cycle-based simulators have better performances provided that architectures are mainly synchronous.

In EDA, Avionic and automotive industries, we need models able to describe these two communication models. Additionally, considering the large number of actors in the design of the very large systems (or even systems of systems) we need standard-based approaches to provide interoperability between models and to cover the whole design flow, from systems requirements to code generation. These models must be precise enough to support various analyses at different refinement levels. In this paper, we focus on two particular standards, AADL (Architecture Analysis & Design Language) [START_REF] Albert | Comparison of Event-Triggered and Time-Triggered Concepts with Regard to Distributed Control Systems[END_REF] standardized by the Society of Automotive Engineers (SAE) and the UML (Unified Modeling Language) profile for MARTE (Modeling and Analysis of Real-Time and Embedded systems) [START_REF] Omg | UML profile for Modeling and Analysis of Real-Time and Embedded systems (MARTE)[END_REF], recently adopted by OMG (Object Management Group). Both standards focus on the modeling and analysis of embedded systems. Both offer constructs to model the application, the execution platform, and to allocate the former on to the latter.

In this paper, we compare their expressiveness to combine periodic computations together with aperiodic ones and to describe the induced prominent communication schemes: event-driven/time-triggered. We particularly emphasize on MARTE Time Model [START_REF] André | Modeling Time(s)[END_REF] that was specifically devised to specify in a formal way timed domains of computation and communication. This is the continuation of some of our previous work ( [START_REF] André | Modeling of Immediate vs. Delayed Data Communications: from AADL to UML MARTE[END_REF], [START_REF] Lee | Dealing with AADL End-to-end Flow Latency with UML MARTE[END_REF]) to compare both formalisms. We illustrate this comparison by using examples devised to show how to perform end-to-end flow latency analysis on AADL models [START_REF] Feiler | Flow Latency Analysis with the Architecture Analysis and Design Language (AADL)[END_REF]. In this example, several threads, periodic or not are connected through event, data or event-data ports. We advocate that AADL introduces avoidable redundancies that, at the very least, obscure the model and may even lead to inconsistencies. We also claim that a subset of UML/MARTE can be combined with AADL to cover a larger scope than the one currently covered by AADL, thus benefiting to AADL users. Such a combination would also benefit to UML/MARTE users because some of their models could then be analyzed by existing AADL tools.

MARTE Time Model

Time and time-related concepts of the UML profile for MARTE are further described in another paper [START_REF] André | Modeling Time(s)[END_REF]. This section recalls the time structure definition and focuses on time relations required to represent the two different kinds of communications: event-triggered or time-triggered.

Definitions

In MARTE, Time can be physical, and considered as continuous or discretized, but it can also be logical, and related to user-defined clocks. Time may even be multiform, allowing different times to progress in a non-uniform fashion, and possibly independently to any (direct) reference to physical time. The time structure is defined by a set of clocks and relations on these clocks. Here clock is not a device used to measure the progress of physical time. It is rather a mathematical object lending itself to formal processing. A clock that refers to physical time is called a chronometric clock. A distinguished chronometric clock called idealClk is provided as part of the MARTE time library. This clock represents the "ideal" physical time used, for instance, in physical and mechanics laws. At the design level most of the clocks are logical ones. For instance, we consider the processor cycle or the bus cycle as been logical clocks. For each clock, we consider an ordered set of instants (, ) where  is an order relation on . Clocks are independent of each others unless some instant relations are imposed. Three kinds of instant relations have been defined: coincidence (), (weak) precedence () and strict precedence (). We have also defined clock relations that are a convenient way to impose many-often infinitely many-instant relations at once. A Time Structure is a set of clocks and the partial ordering relation induced by the instant relations on the clocks.

Event-triggered communications

A task T starts at time Ts and finishes at time Tf.

For each task, we consider the clock ^Ts as the set of instants at which the task T starts and the clock ^Tf as the set of instants at which the task T finishes. A task cannot end before having started and every time a task starts it must end, in one way or another (normal ending, abortion, interrupted).

We use the clock relation alternatesWith (denoted by ) to represent this causality relation between the start and the end (Eq. 1) of tasks. Eq. 1 denotes that

( i* ) ( ^Ts[i]  ^Tf[i]  ^Ts[i+1]
), i.e., every i th instant of ^Ts strictly precedes every i th instant of ^Tf which in turns (weakly) precedes every (i+1) th instant of ^Ts. This relation is not symmetrical and does not assume the task T is periodic.

^Ts alternatesWith ^Tf (^Ts  ^Tf) ( 1
)
This relation is very general and can also represent an event-triggered communication from a task T1 to a task T2. The same relation holds between ^T1f and ^T2s, Task T2 is executed when and as soon as the task T1 completes. In that example, the termination of T1 triggers asynchronously the start of T2. Note that we only have partial orders and i.e., no instant relation is induced between the start or end of T2 and the next start of T1, even though the figure may seem to imply one.

Time-triggered communications

With time-triggered communications, the data is sampled from a buffer according to a triggering condition. We use the clock relation sampledOn to represent this kind of sampling and the triggering condition is given by instants of clocks.

Following our previous example, we replace Eq. 2 by Eq. 3. clk is the sampling condition, i.e., the triggering clock. This clock relation is equivalent to the instant relations of Eq. 4.

^T2s  ^T1f sampledOn clk

(3)

( i* ) ( j,k* ) ( ^T2s[i]  clk[j] && clk[j-1]  ^T1f[k]  clk[j]) (4)
Figure 2 illustrates the use of the clock relation sampledOn. It does not show the start of T1 since it is not relevant here. The start of task T2 is precisely given by the sampling clock clk, however, some events may be missed if the sampling clock is not fast enough.

Periodic tasks and physical time

Logical clocks are infinite sets of instants but we do not assume any periodicity, i.e., the distance between successive instants is not known. The relation discretizedBy is used to discretize idealClk, a dense chronometric (related to physical time) perfect (with no jitter or any other flaw) clock. Eq. 5 shows how to use the clock relation discretizedBy to create a 100Hz clock.

c 100  idealClk discretizedBy 0.01 [START_REF] André | Modeling of Immediate vs. Delayed Data Communications: from AADL to UML MARTE[END_REF] Eq. 5 states that the distance (duration) between two successive instants of clock c 100 is 0.01s. The unit second (s) is implied by the use of idealClk.

A brief AADL overview

AADL supports the modeling of application software components (thread, subprogram, process), execution platform components (bus, memory, processor, device) and the binding of software onto execution platform. Each model element (software or execution platform) must be defined by a type and comes with at least one implementation.

Threads are executed within the context of a process, therefore the process implementations must specify the number of threads it executes and their interconnections.

Type and implementation declarations also provide a set of properties to characterize model elements. For threads, the AADL standard properties include the dispatch protocol (periodic, aperiodic, sporadic, background), the period (if the dispatch protocol is periodic or sporadic), the deadline, the minimum and maximum execution times, along with many others.

AADL end-to-end flows explicitly identify a datastream from sensors to the external environment (actuators). Figure 3 illustrates the example under consideration that derives from [START_REF] Feiler | Flow Latency Analysis with the Architecture Analysis and Design Language (AADL)[END_REF].

This flow starts from a sensor (an aperiodic device instance) and sinks in an actuator (also aperiodic) through three process instances. The first process executes the first two threads and the last thread is executed by the second process. The two devices are part of the execution platform and communicate via a bus (db1) with two processors (cpu1 and cpu2), which host the three processes with several possible bindings. All processes are executed by either the same processor, or any other combination. One possible binding is represented by the arrows on this figure. The component declarations and implementations are not presented here. The full AADL code is available in [START_REF] Feiler | Flow Latency Analysis with the Architecture Analysis and Design Language (AADL)[END_REF].

Figure 3. The example in AADL.

There are three kinds of ports: data, event and event-data. Data ports are for data transmissions without queuing. Connections between data ports are either immediate or delayed. Event ports are for communications of events that may be queued. The size of the queue may induce transfer delays that must be taken into account when performing latency analysis. Event data ports are for message transmission with queuing, here again the queue size may induce transfer delays. On our example, all components have data ports represented as a solid triangle. We have © 2008 IEEE omitted the ports of the process since they are required to be of the same type than the connected port declared within the thread declaration and are therefore redundant.

MARTE for AADL

MARTE is expected to be the basis for UML representation of AADL models [START_REF] Faugère | MARTE: Also an UML Profile for Modeling AADL Applications[END_REF]. The adopted MARTE OMG specification provides guidelines in this direction. The main goal of this paper is to further investigate how specific AADL concepts required for end-to-end flow latency analysis can be represented in MARTE. As such, this work may be integrated in the official 1.0 standard annex.

In this section, we recall the principles we presented previously [START_REF] Lee | Dealing with AADL End-to-end Flow Latency with UML MARTE[END_REF] to build a model library for AADL with MARTE and that should be used as a black block by end-users. The following section illustrates the use of this library on two selected examples. For brevity, we only present model elements required for dealing with our example.

AADL application software components with MARTE.

First, we need to create classifiers to represent AADL threads. In this example, we only need PeriodicThread and AperiodicThread. We use the stereotype SwSchedulableResource from the Software Resource Modeling sub-profile together with its properties deadlineElements and periodElements that help model transformation tools to extract the right property. Periodic threads have a property called period. The MARTE equivalent to the AADL type Time is NFP_Duration, defined in the MARTE::BasicNFP_Types (Non Functional Property Types) model library. An NFP_Duration value is defined as a tuple containing a real value and a time unit, among others.

AADL ports with MARTE.

UML component diagrams provide ports and connectors to connect components. The queuing policy should rather be represented on the algorithm itself, i.e., on a UML activity diagram. Activities are composed of actions. Ordering in which the actions are executed are given by a control flow. Data communications between the actions are represented with object flows. By default, an object flow has a queue, the size of which can be parameterized with its property upperBound. So object flows can be used to represent AADL communications using either event or event-data ports. UML allows the specification of a customized selection policy to select which one among the tokens stored in the object node is read. Unfortunately, the selection behavior must select only one token making it impossible to represent the AADL dequeue protocol AllItems. This protocol dequeues all items from the port every time the port is read. Thus, only the dequeue protocol OneItem is supported.

To model data ports, UML provides DataStore nodes. On these nodes, the tokens are never consumed thus allowing for multiple readings of the same token. Using a data store node with an upper bound equal to one is a good way to represent communications through data ports.

The difference between immediate and delayed communications is addressed in the next sub section, since it is not really a structural matter but rather a temporal aspect.

AADL MoCC with MARTE.

Aside the model elements, the time semantics of these elements must be defined. On one hand, the model of computation, i.e., when the processing starts, finishes or is aborted. On the other hand, the model of communications, i.e., what kind of communication is used. The MARTE Time subprofile, inspired from the theory of tag systems [START_REF] Sangiovanni-Vincentelli | A framework for comparing models of computation[END_REF], provides a set of general mechanisms to define MoCC. These modeling aspects should be hidden to end-users and we show here how to use MARTE, as a model architect, to build a partial MoCC suitable for AADL. Time constraints are specified in MARTE using the stereotype ClockConstraint together with a specification language called Clock Constraint Specification Language (CCSL). The clock relations presented in Section 2 are part of CCSL.

We consider only two kinds of communications, the ones that are possible in AADL. Event-triggered communications and sampled communications. Note that the nature (event, event-data, or data) of the ports involved in the communication is not enough to determine its kind.

For instance, event-triggered communications exist in chains of aperiodic tasks (devices or threads) connected by event or event-data ports. They also exist with periodic tasks connected by data-ports through an immediate connection. In that latter case, the consuming task becomes aperiodic and its execution is triggered by the completion of the producing task. The CCSL clock relation alternatesWith models data-driven communications.

Sampled communications occurs in various cases with data ports and periodic threads with delayed

Two examples

The MARTE representation (fully asynchronous case)

We start by describing the model algorithm with an UML activity diagram (see Figure 4, upper-most part). All communications are through event-data ports with infinite queues. Two actions (acquire and release) have been added as the behavior of devices. AADL software (Figure 4, middle part) components are modeled with MARTE composite structure diagrams using the classifiers defined in Section 4. The bottom layer represents the execution platform (processors and bus). This layer-oriented approach significantly differs from AADL two-layer models and gives flexibility to change one layer independently of the others. AADL models does not consider the pure applicative part and merge this information either within the second or the third level (compare with Figure 3).

The AADL binding mechanism finds its equivalent in the MARTE allocation package. Actions and object nodes are allocated (dashed arrows on Figure 4) to software components.

All threads are aperiodic, therefore all communications are asynchronous and we only use the clock relation alternatesWith (Eq. 6-9). ^Ds alternatesWith ^T1 [START_REF] Lee | Dealing with AADL End-to-end Flow Latency with UML MARTE[END_REF] ^T1 alternatesWith ^T2 [START_REF] Feiler | Flow Latency Analysis with the Architecture Analysis and Design Language (AADL)[END_REF] ^T2 alternatesWith ^T3

^T3 alternatesWith ^Da [START_REF] Sangiovanni-Vincentelli | A framework for comparing models of computation[END_REF] All these annotations (stereotypes) can be extracted using model-driven engineering techniques and fed into time analysis tools, including AADL latency analysis tool. Then, we go a bit further than AADL, by bringing back the latency analysis results into UML and MARTE in the form of timing diagrams (Figure 5). The timing diagram represents a family of possible schedules for a given execution flow and a given pair application/execution platform. Computation execution times (thick horizontal lines) are equal to the latency for devices and range between the MinimumExecutionTime and the Deadline for threads. Oblique lines linking two computation lines represent the communications and the sampling delays. For sampled communications, this amounts to wait for the next tick of the receiver clock. The maximal sampling delay is when the communication waits for the full sampling period because the previous tick has just been missed. It is not normative in UML timing diagrams to have these "oblique" lines, but it is a convenient notation to represent intermediate communication states between two steady processing states (e.g., between Ds and t1). Assuming, as in [START_REF] Feiler | Flow Latency Analysis with the Architecture Analysis and Design Language (AADL)[END_REF], that the sampling delays are always maximal, we get the same formulas (reproduced below) as the AADL latency analysis tool.

The MARTE representation (Mixed Event-data flow case)

We study here a second possible configuration extracted from [START_REF] Feiler | Flow Latency Analysis with the Architecture Analysis and Design Language (AADL)[END_REF] that only differs by making periodic the thread t2 (Figure 6). Few other cases involving data ports are studied in [START_REF] André | Modeling of Immediate vs. Delayed Data Communications: from AADL to UML MARTE[END_REF]. We also get a different timing diagram and different flow latency formulas (not shown here).

Conclusion

AADL offers lots of features very important to model and analyze computations and communications of embedded systems. However, combining all of these features without a guideline (not part of the standard) can lead to model completely meaningless and impossible to analyze. We have shown how the MARTE Time model could be used to have the same expressiveness with less modeling concepts. More generally, MARTE and its time model could be used to model various timed models of computation and communication.

We also think that it is important to have specifications free, as much as possible, of implementation choices (platform independent models). To achieve this goal, we need model elements of a higher level of abstraction than AADL threads. AADL two-level models assume that part of the application has already been allocated to a software execution platform made of threads. Our approach makes that allocation explicit when required. We propose to use for that purpose UML activities. Making a link to the software execution platform (runtime executive) is not a refinement but rather an allocation. The former implies models of the same nature, whereas the latter make links between models of different natures. If system level models are needed, we can use SysML for that purpose. UML activities integrate very well with SysML models and some experimentations [START_REF] Demathieu | Experimentation of MARTE in the Industry[END_REF] have shown that MARTE can also be used at system level together with SysML.

Figure 1

 1 illustrates graphically the clock relation alternatesWith. Horizontal lines represent the clocks and their instants. Vertical lines are coincidence relations. Dashed arrows with a filled triangle as an arrowhead are strict precedence relations whereas arrows with a hollow triangle as an arrowhead are (weak) precedence relations. The precedence relations are directly induced by the three clock relations alternatesWith given on the right-hand side of the figure.

Figure 1 .

 1 Figure 1. The clock relation alternatesWith.

Figure 2 .

 2 Figure 2. ^T2s  ^T1f sampledOn clk.

Figure 4 .

 4 Figure 4. MARTE model, fully aperiodic case.

Figure 5 .

 5 Figure 5. Timing diagrams, all aperiodic case.

Forum

  on specification, verification & Design Languages (FDL'08), 23-25 September 2008, pp. 154-159 The original publication is available in IEEE explore (http://dx.doi.org/10.1109/FDL.2008.4641438) © 2008 IEEE

Figure 6 .

 6 Figure 6. MARTE model, mixed case.Induced clock relations and a graphical representation of the implied partial ordering is illustrated in Figure7. The communication from step1 to step2 becomes a sampled communication and is shown as a plain (green) arrow.clk  idealClk discretizedBy P[START_REF] Demathieu | Experimentation of MARTE in the Industry[END_REF] Eq. 10 declares the periodic clock that triggers t2.

Figure 7 .

 7 Figure 7. CCSL constraints, mixed case.

  The combination of Forum on specification, verification & Design Languages (FDL'08), 23-25 September 2008, pp. 154-159 The original publication is available in IEEE explore (http://dx.doi.org/10.1109/FDL.2008.4641438)

	various parameters induces either asynchronous or
	sampled communications.
	© 2008 IEEE

  Forum on specification, verification & Design Languages (FDL'08), 23-25 September 2008, pp. 154-159 The original publication is available in IEEE explore (http://dx.doi.org/10.1109/FDL.2008.4641438)

© 2008 IEEE

  Forum on specification, verification & Design Languages (FDL'08), 23-25 September 2008, pp. 154-159 The original publication is available in IEEE explore (http://dx.doi.org/10.1109/FDL.2008.4641438)

  Forum on specification, verification & Design Languages (FDL'08), 23-25 September 2008, pp. 154-159 The original publication is available in IEEE explore (http://dx.doi.org/10.1109/FDL.2008.4641438)

	communications or with immediate communications in
	case of oversampling.
	© 2008 IEEE


Nadia Curie Dempowo

François Valdivieso

Julien Bruchon

Jacques Lechelle

A ONE-DIMENSIONAL MODELING OF U/Pu INTERDIFFUSION IN MIXED OXIDES (U,Pu)O2 AND ANALYZIS OF HOMOGENIZATION WITHIN A MICROSTRUCTURE USING A FINITE DIFFERENCE METHOD

Keywords: Modeling, Interdiffusion, MOX, homogenization, microstructure, finite difference

One-dimensional modeling of bulk cationic interdiffusion

in mixed oxides (U,Pu)O_2 and analysis of homogenization within a microstructure using a finite difference method

Nadia Dempowo, François Valdivieso, Julien Bruchon,

Introduction

Sintering is a process frequently encountered in the manufacturing process of ceramics materials. The sintering of ceramics constituted by only one chemical element is a non-reactive sintering. However, very often, ceramics are constituted by more than one chemical element and are therefore heterogeneous from a chemical point of view. Reactive sintering has long been successfully experimentally used in the field of traditional ceramics and has been developed for technical ceramics since the last quarter of the last century. Mullite-zirconia is often used as a model system [START_REF] Slamovich | Densification Behavior of Single-Crystal and Polycrystalline Spherical Particles of Zirconia[END_REF][START_REF] Descamps | High-Temperature Characterization of Reaction-Sintered Mullite-Zirconia Composites[END_REF][START_REF] Prochazaka | Microstructure of Sintered Mullite-Zirconia Composites[END_REF][START_REF] Raghdi | Mullite-zirconia composites prepared from halloysite reaction sintered with boehmite and zirconia[END_REF][START_REF] Duran | Templated grain growth of textured mullite/zirconia composites[END_REF]. Sintering is a key step in the manufacturing process of mixed oxides used as fuels in nuclear reactors. Since Pu has been reprocessed as MOX fuel, the sintering stage has become a reactive sintering in which a chemical reaction between the initial constituents of the green compact (made of UO2 and PuO2 powders) and its densification occurs during the same thermal treatment. So far, in order to construct a sintering model, all numerical studies have supposed that the ceramics are homogeneous [START_REF] Rojek | Modeling of a Sintering Process at Various Scales[END_REF][START_REF] Pino-Muñoz | Sintering at Particle Scale: An Eulerian Computing Framework to Deal with Strong Topological and Material Discontinuities[END_REF][START_REF] Wakai | Sintering mechanics of ceramics: a short review[END_REF][START_REF] Zhang | Sintering of two particles by surface and grain boundary diffusiona threedimensional model and a numerical study[END_REF][START_REF] Zhang | The sintering of two particles by surface and grain boundary diffusion-a twodimensional numerical study[END_REF]. Unlike these latter, the sintering of mixed oxides is reactive, due to the number of chemical elements in presence, precisely uranium, plutonium and oxygen. A sub-granular model for solid state free sintering of mixed oxides has been developed in order to follow the evolution of two grains in references [START_REF] Ajdour | Développement d'un code de calcul pour la simulation du frittage en phase solide[END_REF][START_REF] Léchelle | A sub-granular scale model for solid state free sintering: results on the evolution of two grains[END_REF][START_REF] Léchelle | A mechanistic approach of the sintering of nuclear fuel ceramics[END_REF]. This model, named SALAMMBO, describes the behavior of two grains by solving Fick's laws with fluxes arising from different diffusion mechanisms (grain boundary, surface and volume). In addition to the chemical origin of diffusion, a mechanical part leading to shrinkage is considered. In fact, the flux is written 𝑗 ⃗ = -𝐿𝛻 ⃗⃗ µ 𝑇 . The chemical potential in this equation consists of a mechanical and a chemical part. This means µ= µmecha + µchem with µmecha = -Ω 𝟐 (𝜎 ̅ ̅ : ℇ ̅ ̅ ) where 𝜎 ̅ ̅ stands for the stress field whereas ℇ ̅ ̅ stands for the strain field and µchem is constant (taken as zero) when the material is homogeneous (current SALAMMBO consideration) and varying when the material is heterogeneous (MOX case). The Navier-Lame equations are solved at each time step in the bulk of each grain, providing the 3D-displacement field and hence the density of the elastic mechanical energy. The gradient of energy is the driving force of diffusion phenomena occurring during densification. As any other sintering model, this model was first developed for the simplified case of a chemically homogeneous material. As (U,Pu)O2 ceramics are initially heterogeneous and homogeneity of nuclear is required to avoid localized high burnup areas during irradiation, it is necessary to take into account the bulk diffusion of each of the three chemical elements in order to construct a reactive sintering model of (U,Pu)O2. The understanding of this reactive sintering can help to improve cationic homogenization experimental conditions. To do so, the interdiffusion mechanism which is one of the physical phenomena involved and which conditions the formation of the (U,Pu)O2 solid solution requires attention. This reactive sintering supposes to solve Fick's laws not only for each diffusion mechanism, but also for each chemical element. Some studies in the literature deal with the resolution of Fick's laws in the binary or ternary systems depending on their aims. In a binary system, only one diffusion coefficient comes into play. In contrast, in a ternary system in which two or more cations are being transported, new complications arise : a single interdiffusion coefficient is no longer appropriate to describe the diffusion process [START_REF] Darken | Application of the Gibbs-Duhem Equation to Ternary and Multicomponent Systems[END_REF][START_REF] Cussler | Multicomponent diffusion -Chemical Engineering Monographs[END_REF]. This makes the use of a diffusion matrix necessary [START_REF] Philibert | Diffusion et transport de matière dans les solides[END_REF]. The diffusion matrix contains four diffusion coefficients namely the "main-term" diffusion coefficients (for diagonal ones) and "cross-term" diffusion coefficients (for offdiagonal ones). In general, the cross-terms are not symmetric. Some works reported in the literature neglect the cross-terms [START_REF] Kajihara | Analysis of dissolution of alpha phase in gamma/alpha/gamma Diffusion couples of the Fe-Cr-Ni system using analytical solutions for semi-infinite diffusion couples[END_REF][START_REF] Vitek | Numerical modeling of diffusion-controlled phase transformations in ternary systems and application to the ferrite/austenite transformation in the Fe-Cr-Ni system[END_REF]. Such an approximation cannot be made for U-Pu interdiffusion because these cations diffuse on the same cationic sites of the fluorine structure with a trace amount of cationic vacancies. A non-zero value of a cross-term gives a measure of a flux of one component generated by the concentration gradient of a second component. That means the multicomponent diffusion is a non-linear problem and the interactions between components are significant.

The The

Methodology 2.1. Multicomponent flux equations: the ternary system (U,Pu,O)

Fick's laws expressing the fluxes of species account for diffusion phenomena. In the context of interdiffusion, the process consists in bringing two single-crystals into contact. Each single-crystal has its own initial chemical composition.

𝝏𝒄 𝒌 𝝏𝒕 = -𝜵. 𝑱 𝒌 (2)
where Jk represents the diffusional flux of the chemical element k; n stands for the number of chemical element in the system; ck is the concentration of the chemical element k; and 𝑫 𝒌𝒋 𝒏 is the diffusivity that can be expressed as [START_REF] Andersson | Models for numerical treatment of multicomponent diffusion in simple phases[END_REF]:

𝐷 𝑘𝑗 𝑛 = 𝐷 𝑘𝑗 -𝐷 𝑘𝑛 ( 𝑉 𝑗 𝑉 𝑛 ) (3)
where Vj represents the partial volume of chemical element j and Vn is the total volume of the pattern (i.e one mole of U1-xPuxO2).

In addition,

𝐷 𝑘𝑗 = ∑ (𝛿 𝑖𝑘 -𝑥 𝑘 )𝑥 𝑖 𝑀 𝑖 𝜕𝞵 𝑖 𝜕𝑥 𝑗 𝑛 𝑖=1 (4)
Here 𝑥 𝑘 is the mole fraction of chemical element k, 𝛿 𝑖𝑘 represents the Kronecker symbol, equal to 1 if i=k, and 0 to otherwise.

𝜕𝞵 𝑖 𝜕𝑥 𝑗

is the partial derivative of the chemical potential of the species i with respect to the mole fraction xi and corresponds to the thermodynamic factor of Darken. This quantity can therefore be evaluated from the thermodynamic description of the system. The expressions for the chemical potentials must be given in the form µ 𝑘 (𝑥 1 , 𝑥 2 , … , 𝑥 𝑛 ). Mi represents the mobility of element i. In this framework, there is no question of determining diffusion coefficients from mobilities and the thermodynamic model of the solid solution; these values will be obtained by interpolating experimental data.

𝑱𝒌= -∑ 𝑫 𝒌𝒋 𝒏 𝜵𝒄 𝒋 𝒏-𝟏 𝒋=𝟏

, which gives the expression of flux, has a limit: one of the elements is not taken into account into the summation. The choice of the last element that is excluded in the computation of interdiffusion flux depends on the type and the characteristics of the system. Generally, in the systems build up of solute and solvent, the flux of solvent is neglected. So, we will determine in the n-component system, (n-1) interdiffusion fluxes and (n-1) interdiffusion coefficients. In addition, in the volume center frame, the fluxes are related to each other and we can numerically deduce the flux of the last element (ncomponent) as ∑ 𝑽 𝒊 𝑱 𝒊 𝟑 𝒊=𝟏 =𝟎 [START_REF] Philibert | Diffusion et transport de matière dans les solides[END_REF]. Therefore, as explained in the introduction, Fick's laws will apply on U and Pu and the flux of O could be deduced if necessary from this relation.

In order to simulate 1D interdiffusion, the two Fick's laws are combined to determine the concentrations, leading to: The values of interdiffusion coefficients are obtained by interpolating Matzke graph [START_REF] Matzke | Atomic Transport Properties in U 0 2and Mixed Oxides (U,Pu)O2[END_REF] that gives the representation of interdiffusion coefficient at T=1600°C with an oxygen potential of ~90 kcal/mol (~ -375kJ/mol). The magnitude order is around 10 -4 cm²/s. Thus, for the modeling of the interdiffusion between the couple UO2/PuO2, the following diffusion coefficients will be used: 𝐷 UU O = 6*10 -14 𝑐𝑚 2 /𝑠, 𝐷 PuPu O = 12*10 -14 𝑐𝑚 2 /𝑠, 𝐷 UPu O =10 -14 𝑐𝑚 2 /𝑠, 𝐷 PuU O = 2*10 -14 𝑐𝑚 2 /𝑠.

𝜕𝑐 𝑘 𝜕𝑡 = ∑ 𝐷 𝑘𝑗 𝑛 𝜕 2 𝑐 𝑗 𝜕𝑥 2 𝑛-1 𝑗=1 (
Furthermore, for the modeling of the interdiffusion between the couple of two pellets UO2/U0.55Pu0.45O2, the magnitude order provided by the plot of Noyau [START_REF] Noyau | Etude des phénomènes d'autodiffusion et d'interdiffusion du plutonium dans les céramiques de type U1-yPuyO2±x[END_REF], for a value of oxygen potential of -395 kJ.mol -1 , will be used. The values that will be considered to reproduce the concentration profile of UO2 are: 2 /𝑠 whereas those that will reproduce the concentration profile of U0.55Pu0.45O2 are:

D UU O = 10 -15 𝑐𝑚 2 /𝑠, D UPu O =3*10 -16 𝑐𝑚 2 /𝑠, D PuU O =5*10 -16 𝑐𝑚 2 /𝑠, D PuPu O =10 -16 𝑐𝑚
D UU O = 10 -14 𝑐𝑚 2 /𝑠, D UPu O =10 -15 𝑐𝑚 2 /𝑠, D PuU O =8*10 -15 𝑐𝑚 2 /𝑠, D PuPu O =12 * 10 -15 𝑐𝑚 2 .

Initial Conditions

Like every time-dependent problem, the first step is the definition of initial conditions, which define the composition of each single-crystal (each diffusion half-couple).

 Case of UO2/PuO2

The first half-couple is uranium dioxide. In the fluorine structure (case of UO2), the expression of concentration is:

c= 𝑛 𝑉 = 𝑧 𝑁 𝑎 * 𝑎 3
c is the concentration of U in the fluorine structure, z is the number of formula units in the conventional cell and its value is 4 for the UO2 system, Na is the Avogadro number. The lattice parameter a is 5.47Å at sintering temperature (around 1600 or 1700°C) [START_REF] Andersson | Models for numerical treatment of multicomponent diffusion in simple phases[END_REF]. The corresponding concentration of U is 0.04 mol/cm 3 . By analogy, the lattice parameter of PuO2 is 5.396Å and the concentration of Pu is 0.042 mol/cm 3 . Those initial conditions are applied to solve the previous equations.

 Case of UO2/U0.55Pu0.45O2

The initial conditions of this couple are also required because they will be used at , where y is the Pu content, Vm the molar partial volume, Vm = 𝑎 3 (1 + 𝛼 𝑇 ) 3 𝑁 𝑎 . The lattice parameter a is given by a = 5.470 -0.074y +(0.301 + 0.11y)x Å [START_REF]Thermodynamic and transport properties of Uranium dioxide and related phases[END_REF], with x the deviation from oxygen stoichiometry.

Taking into account a 45at%Pu content, the lattice parameter becomes 5.436Å. The thermal expansion coefficient αT, for this non-stoichiometric component is: α T = 1.1833 × 10 -5 -5.013 × 10 -9 𝑇 + 3.756 × 10 -12 𝑇 2 -6.125×10 -17 𝑇 3 for 923≤T≤3120 K, leading to 𝑐 𝑈 = 0.02273 mol/cm 3 𝑐 𝑃𝑢 = 0.0186 mol/cm 3 .

Boundary conditions

Boundary conditions are required to close previous system of equations. The choice of the type of boundary conditions depends on the type of problem to be solved.

 Dirichlet boundary conditions

For the sintering of large single-crystals, because interdiffusion zone is smaller than the sample, Dirichlet boundary conditions will be first used. Note that the finite difference scheme that will be used later does not define the concentration values at the two external nodes, corresponding to k=1 and k=n. The Dirichlet boundary conditions considered here assume that the concentration in U and Pu does not change at the boundaries of spatial domain.

 Neumann boundary conditions

The Neumann boundary conditions consists in imposing the normal flux of the species at the boundary of the domain. Here, a flux equal to zero is considered, which means, in a 1D-situation, that no matter moves through the ends of the sample.

Numerical resolution of the problem: the finite difference discretization

This section describes the numerical method used for solving multicomponent diffusion partial differential equations. Diffpack library [START_REF] Langtanten | [END_REF] is used in order to solve the resulting linear system of equations. The purpose is to analyze the evolution of concentration profiles of U and Pu chemical elements.

A finite difference scheme is applied to Equation 8, in order to approximate the unknown c' and then rebuild c. More precisely, the time domain [0, T] is discretized by a set of ordered points, t0=0, t1,…,tn, equally spaced. The time step is then defined by 𝛥𝑡 = 𝑡 𝑙 -𝑡 𝑙-1 for any l> 0. Similarly, the space domain is discretized by a set of equidistributed points x1 = 0, x2

(𝑃𝑢) = 𝜆 𝑃𝑢 × 𝛥𝑡 𝛥𝑥 2 𝑐 𝑘-1,𝑙 ′ (𝑃𝑢) + (1 -2𝜆 𝑃𝑢 × 𝛥𝑡 𝛥𝑥 2 )𝑐 𝑘,𝑙 ′ (𝑃𝑢) + 𝜆 𝑃𝑢 × 𝛥𝑡 𝛥𝑥 2 𝑐 𝑘+1,𝑙 ′ (𝑃𝑢) 2≤k≤n-1 (8) { 𝑐 𝑈 = 𝑏 𝑈𝑈 𝑐 𝑈 ′ + 𝑏 𝑈𝑃𝑢 𝑐 𝑃𝑢 ′ 𝑐 𝑃𝑢 = 𝑏 𝑃𝑢𝑈 𝑐 𝑈 ′ + 𝑏 𝑃𝑢𝑃𝑢 𝑐 𝑃𝑢 ′ (9)
One of the most important steps of the resolution of time-dependent problem is the choice of the time step. In the explicit schemes, the Courant-Friedrichs-Lewy (CFL) condition requires that the number of Courant be less than 0.5. This number is defined for a one-dimensional system by: C = 𝜆 𝑈 × 𝛥𝑡 𝛥𝑥 2 or C=𝜆 𝑃𝑢 × 𝛥𝑡 𝛥𝑥 2 . Thus, the value assigned to C is 0.3. Knowing the value of 𝜆 𝑈 and 𝜆 𝑃𝑢 , the value of time step is deduced, which is the minimum between . The value of space step is 𝛥x = L/(n-1), with L the length throughout the two half-couples and n is the total number of discretization points. This finite difference scheme is associated with the Dirichlet boundary conditions, which can be written in terms of the finite difference scheme as follows:

𝑐 1,𝑙+1 ′ (U)= 𝑐 1,𝑙 ′ (U) 𝑐 1,𝑙+1 ′ (Pu)= 𝑐 1,𝑙 ′ (Pu) 𝑐 1,𝑙+1 (U)= 𝑏 𝑈𝑈 ×𝑐 1,𝑙+1 ′ (U)+ 𝑏 𝑈𝑃𝑢 ×𝑐 1,𝑙+1 ′ (Pu) 𝑐 1,𝑙+1 (Pu) =𝑏 𝑃𝑢𝑈 ×𝑐 1,𝑙+1 ′ (U)+𝑏 𝑃𝑢𝑃𝑢 ×𝑐 1,𝑙+1 ′ (Pu) 𝑐 𝑛,𝑙+1 ′ (U)= 𝑐 𝑛,𝑙 ′ (U) 𝑐 𝑛,𝑙+1 ′ (Pu)= 𝑐 1,𝑙 ′ (Pu) 𝑐 𝑛,𝑙+1 (U)= 𝑏 𝑈𝑈 ×𝑐 𝑛,𝑙+1 ′ (U)+ 𝑏 𝑈𝑃𝑢 ×𝑐 𝑛,𝑙+1 ′ (Pu) 𝑐 1,𝑙+1 (Pu) =𝑏 𝑃𝑢𝑈 ×𝑐 𝑛,𝑙+1 ′ (U)+𝑏 𝑃𝑢𝑃𝑢 ×𝑐 𝑛,𝑙+1 ′ (Pu)
Therefore, the system is solved numerically and the unknown concentrations are then deduced.

Concerning the numerical implementation based on Neumann boundary conditions, the mathematical equations can be written:

𝛁𝑐 = 0 at k=1 and k=n . This leads to 𝛁𝑐′ = 0 at k=1 and k=n. Assuming that at the end point at k = 1, the incoming flux is equal to the outgoing flux:

𝑐 1,𝑙+1 ′ (U) = 𝜆 𝑈 × 𝛥𝑡 𝛥𝑥 2 (𝑐 0,𝑙 ′ (U) -2×c'1,l(U) + c'2,l(U) ) + 𝜆 𝑈 × 𝛥𝑡 𝛥𝑥 2 𝑐 1,𝑙 ′ (𝑈).
The condition of flux null is applied between the point k=1 and a point out of the domain k=0. Assuming that the difference between the flux at k=1 and at the other outer points on the left is zero at all times. This means that c'1,l(U) -c'0,l(U) = 0. The concentration 𝑐 0,𝑙 ′ is the concentration at a point out of the domain. By replacing c'0,l(U) = c'1,l(U) in the previous relation, the deduction appears

𝑐 1,𝑙+1 ′ (U) = (1 -𝜆 𝑈 × 𝛥𝑡 𝛥𝑥 2 ) 𝑐 1,𝑙 ′ (U) + 𝜆 𝑈 × 𝛥𝑡 𝛥𝑥 2 𝑐 2,𝑙 ′ (U).
The process is the same for point k = n.

Thus, at the first point corresponding to k=1, 𝑐 1,𝑙+1

′ (U) = (1 -𝜆 𝑈 × 𝛥𝑡 𝛥𝑥 2 ) 𝑐 1,𝑙 ′ (U) + 𝜆 𝑈 × 𝛥𝑡 𝛥𝑥 2 𝑐 2,𝑙 ′ (U)
At the last one corresponding to k=n, 𝑐 𝑛,𝑙+1

′ (U) = 𝜆 𝑈 × 𝛥𝑡 𝛥𝑥 2 𝑐 𝑛-1,𝑙 ′ (U) + (1 -𝜆 𝑈 × 𝛥𝑡 𝛥𝑥 2 )𝑐 𝑛,𝑙 ′ (U) By analogy, 𝑐 1,𝑙+1 ′ (Pu) = (1 -𝜆 𝑃𝑢 × 𝛥𝑡 𝛥𝑥 2 ) 𝑐 1,𝑙 ′ (Pu) + 𝜆 𝑃𝑢 × 𝛥𝑡 𝛥𝑥 2 𝑐 2,𝑙 ′ (Pu) 𝑐 𝑛,𝑙+1 ′ (Pu) = 𝜆 𝑃𝑢 × 𝛥𝑡 𝛥𝑥 2 𝑐 𝑛-1,𝑙 ′ (Pu) + (1 -𝜆 𝑃𝑢 × 𝛥𝑡 𝛥𝑥 2 )𝑐 𝑛,𝑙 ′ (Pu)
These resulting equations associated with the finite difference scheme lead to a numerical system that is solved using the functionalities of the Diffpack library for the finite difference method.

Results and discussion

simulation of two single-crystals of UO2 and PuO2

So far, the physical model and the numerical method necessary to solve the interdiffusion phenomenon between two singlecrystals in contact were presented. The purpose is to find the unknown concentrations of U and Pu as a function of time. As mentioned at the beginning of this paper, we aim at analyzing the change in concentration of each chemical element during interdiffusion in a system constituted by two half-couples of single-crystals. One made of UO2 and the other PuO2. Based on the fact that the size of these crystallites as shown by scanning electron microscope or deduced from X-Rays peak width in diffraction patterns, lies between 100 and 300nm, we have simulated the interdiffusion using the size of 100nm. Thus, the entire length of the two half-couples is 200nm. The initial interface separating the two half-couples is located in the center of the domain.

: Concentration profiles of U and Pu as a function of distance between the two half-couples of interdiffusion using Neumann boundary conditions

The same remarks as previously mentionned apply: the diffusion effect appears clearly, while the mole number of each species is well preserved. Since the space derivative of both concentrations are set to zero at the boundaries of the domain, the steady state is a state of constant concentration. The Neumann boundary conditions give a good description of the physical reality of interdiffusion at the crystallite scale during sintering when the interdiffusion zone is significant in comparison to the total interdiffusion couple length.

Validation of the model based on Dirichlet boundary conditions

The Dirichlet boundary conditions used to simulate the interdiffusion between two single-crystals do not represent the reality of phenomena because the size of single-crystals and the interdiffusion zone are of the same order of magnitude. In contrast, we can use this type of boundary conditions when the interdiffusion zone is smaller than the entire length constituted by the two half-couples, for example the case of the single-crystals with bigger size. Therefore, in order to validate our model, we will try to reproduce the results acquired during an experiment based on two interdiffusion half-couples, one consisting of a pure UO2 half pellet, the other one consisting of a half-pellet of the mixed oxide U0.55Pu0.45O2. In this aim, the initial conditions as described at Section 2.3 are used. The required interdiffusion coefficients were presented with the physical model (see Section 2.2).

So far, the only work in the literature dealing with the concentration profiles of UO2 and U0.55Pu0.45O2 is that of Noyau [START_REF] Noyau | Etude des phénomènes d'autodiffusion et d'interdiffusion du plutonium dans les céramiques de type U1-yPuyO2±x[END_REF]. He plotted the evolution of concentration as a function of distance along the interdiffusion couple. To cope with the very limited interdiffusion zone which extends from 2 up to 5 µm, the author selected an acquisition line for concentration measurements by means of a Castaing EPMA (Electron Probe Micro Analyzer) making a very small 𝛼 angle with the bonding interface of the two half-pellets. He managed in this way to obtain a large number of points with varying concentrations within the interdiffusion zone. The total distance between extreme acquisition points measured perpendicular to the bonding interface remains very small compared to the total length of the half-couples. Therefore, Dirichlet boundary conditions can be applied to these two half-couples to simulate U-Pu interdiffusion.

For a temperature of 1700°C and 50 hours of interdiffusion, the inclination angle is 𝛼=2,048 × 10 -2 rad. In the second case where the duration of interdiffusion is 143 hours, the inclination angle is 3,175 × 10 -2 rad. Figure 3 shows the two half-couples as well as the position of the acquisition line in relation to the bonding interface separating the two half-couples as described by Noyau [START_REF] Noyau | Etude des phénomènes d'autodiffusion et d'interdiffusion du plutonium dans les céramiques de type U1-yPuyO2±x[END_REF]. The interdiffusion between UO2 and U0.55Pu0.45O2 has been simulated for two durations of interdiffusion: 50 and 143 hours. The interdiffusion profiles are then compared to the experimental results of Noyau as shown on Figure 4.

Figure 4 : Concentration profiles of elements U and Pu at 1700°C after 50 hours and after 143 hours

In each figure, a good agreement can be observed between experiments and simulation, although some peaks can be observed on experimental results, probably due to grains boundary diffusion in the polycrystalline materials used. Such an agreement with experimental results validates the numerical model based on Dirichlet boundary conditions used for pellets or very large single-crystals in which interdiffusion zone is smaller than the total sample size.

Numerical study of homogenization within two microstructures

The second part of the results obtained in this paper is the evolution and the achievement of the homogenization in two different microstructures. The realistic simulation of sintering can be extended to more than two single-crystals or grains; it is a question of analyzing the homogenization phenomenon upon a set of grains consisting initially of UO2 and PuO2. The aim of this part is firstly, to better qualify this phenomenon with a given criterion of homogeneity and secondly, to compare the time of homogenization in different types of microstructures. Because of small size of single-crystals and in order to solve correctly the numerical problem, the Neumann boundary conditions will be used. Neumann conditions result in the nullity of flux in the center of each ending single-crystals, which are symmetry planes. Two types of microstructures will be studied: a unimodal distribution microstructure and a non-unimodal distribution microstructure of UO2 and PuO2.

 A unimodal distribution of UO2 and unimodal distribution of PuO2

During the step of grinding (step prior to sintering), the mechanical properties playing a role for UO2 crystallites are different from those of PuO2. As a consequence, UO2 crystallites are generally larger in size than those of PuO2. The first microstructure under consideration that appears is a set of single-crystals with various sizes: 300nm for UO2 single-crystals and 100nm for PuO2. Each type of assembly consists of four single-crystals. The size of the single-crystals at the ends is half of that of the corresponding crystal in order to represent a periodic infinite microstructure. For this microstructure, the concentration profiles are shown in : concentration profiles at different times After 5 minutes, the concentration, which was initially discontinuous at the grain boundary and constant in each grain, becomes smooth and remains periodic. This behavior highlights the effect of the diffusion phenomenon. A distribution of the concentration of each element upon the entire sample length is observed over time. After 1 hour, the profiles present a repartition of each chemical element, which can indicate that homogenization takes place. After 4 hours, a "steady state" seems to be reached. However, the time of the exact achievment of a satisfactory homogeneous mixing cannot be stated with certainty. To reach this goal, a homogeneity criterion will be used.

The homogeneity criterion depends on the scale of observation. In the 1D-representation of this work, this scale, hk, can simply be specified with respect to the space step ∆x: hk = 2 k × ∆x, with 0 ≤ k ≤

𝑙𝑛(𝑛) 𝑙𝑛(2)
-1, n is the number of points of the spatial discretization. In the following, the scale of observation will be the scale given by the spatial discretization, i.e. h0 = ∆x. At this scale, the standard deviation , i.e. the deviation from the average of concentration, is proposed as criterion. It has the same unity as the concentration, and its meaning is easy to interpret. For a ∆x scale, the standard deviation is written as:

𝜎(∆𝑥) = √ 1 𝑛 ∑ (𝑐 𝑈 (ⅈ) -𝑐̅ 𝑈 ) 2 𝑛 𝑖=1 (10)
The plot of  on a logarithm scale versus time on a linear scale leads to the graphs presented over the first four hours in Figure 6: . The plateau tends towards a horizontal asymptote close to 0 and not 0 due to the local truncation error of the numerical scheme [START_REF] Langtangen | Finite Difference Computing with PDEs: A Modern Software Approach[END_REF]. This kind of error is frequently encountered in the numerical resolutions. If the spatial discretization step increases (2∆x for example), the value of  will increase. The curve of  quickly becomes flat and the plateau corresponds to a  value (4.31*10 -3 mol/cm 3 ) greater than the scaled value ∆x. The representations of sigma versus time for a discretization of 2∆x ( Figure 7) can be used to assert that there is no consistency error. It emerges from these analyses that for this specific size and unimodal distribution of UO2 as well as unimodal distribution of PuO2, a significant homogenization occurs and ends within the duration of sintering. Homogeneity is achieved later than in the previous case. It comes from the fact that the mean crystallite size has increased whereas the diffusion coefficient remains the same. The value of sigma at the plateau is 2.35*10 -5 mol/cm 3 , i.e. smaller than its value for a unimodal distribution of UO2 and unimodal distribution of PuO2. One of the interesting remarks is that this significant homogenization occurs during the duration of the sintering process. Therefore, the increase of the size of a type of single-crystal (PuO2 in this case) delays the time to reach a significant homogenization.

 A non-unimodal distribution of UO2 and PuO2

The second microstructure consists of single-crystals with various sizes. The study will be divided into two parts: one consist of non-unimodal distribution with 25%Pu content and the other with 45% Pu content. Here, the Pu content stands for the proportion of PuO2 within the entire length. The aim is to compare the time to reach a significant homogenization within these two microstructures.

That of 25% consists of six single-crystals of UO2 and three of PuO2. Two of UO2 measure 300nm in size and while the two others measure 600nm in size. The two single-crystals of UO2 with 600nm are represented at the ends. On the other hand, all PuO2 crystallites measure 100nm in size. The numerical method (finite difference scheme), the initial conditions (cU = 0.04mol/l and cPu = 0.042mol/l) and the boundary conditions (Neumann boundary conditions) are identical as in the previous section.

The plots are presented on Figure 10: Figure 12 shows the standard deviation of these two non-unimodal microstructures in order to compare them quantitatively.

Figure 12 : Standard deviation of non-unimodal distribution microstructure with 45%Pu and 25%Pu content

The behavior of the concentration standard deviation curve of the two microstructures are identical. However, the microstructure with 25%Pu content decreases more quickly compared to that with a 45% Pu content. The value of standard of microstructure with 25%Pu content is 7.68*10 -4 mol/cm 3 whereas the other one is 2,66 * 10 -3 mol/cm 3 . Note that, at these two values, a significant homogenization is not yet achieved. Using the relation x = √𝐷𝑡 which relates the characteristic distance and interdiffusion time, the time necessary to achieve interdiffusion can be calculated. In this study, x =√𝜆𝑡, where 𝜆 is the smallest value between 𝜆U and 𝜆Pu. Thus, for the size of 300nm, the time is: t = (300 × 10 -7 ) 2 / 5.68338 × 10 -14 ≈ 4.4 hours. Microstructures showing several single-crystals with large size increases the achieving homogenization time. From these analyses, the influence of the grain size can be one of the cause of this later achievement of a significant homogeneity.

Conclusion

The paper has dealt with the interdiffusion phenomena observed during mixed oxides sintering. The multicomponent diffusion principle and the numerical method used to simulate U and Pu diffusion in (U,Pu)O2 were presented. We first noted a good agreement between our work using Dirichlet boundary conditions and the experimental work by Noyau with half-pellets used as interdiffusion couples, which allows us to validate the simulation of interdiffusion for large size samples. From a mathematical point of view, we have shown that taking into account the zero-flux Neumann boundary conditions is a good representation of the reality of interdiffusion (mass conservation) at the scale of crystallites. The numerical resolution was carried out using a finite difference scheme. Then, the second part of this study of homogenization within different types of microstructures allowed to highlight and use a homogeneity criterion which is the standard deviation of one of the diffusing species. The standard deviation of concentration was useful to quantify the achievement of the homogeneity within a microstructure. Finally, in this entire framework, only volume diffusion coefficients were used, since the simulations were carried out in one dimension. In a future 2D and 3D-work, we will also take into account grain-boundary diffusion, which is not negligible at lower temperatures and contributes strongly to diffusion phenomena. This part will constitute the next goals in the development of our reactive sintering model. This scheme of simulation of interdiffusion will be implemented together with the mechanical model into the sintering model (SALAMMBO) developed in our laboratory in order to able to simulate the reactive sintering of MOX fuel.

  the step of validation of our model based on Dirichlet boundary conditions. The values of the initial concentrations of the first half-couple UO2 at the left-hand side are identical to those shown above: cU = 0.04mol/cm 3 and cPu = 0. Concerning the half-couple U0.55Pu0.45O2 on the right-hand side, the concentrations of U and Pu are respectively: 𝑐 𝑈 =

Figure 1 Figure 1 :

 11 Figure 1 : Concentration profiles of U and Pu as a function of distance between the two half-couples of interdiffusion using Dirichlet boundary conditionsThis Figure shows the interdiffusion phenomenon. Initially, UO2 occupies only the left half-couple, and PuO2 occupies the right hand half (Figure1.a). These sharp concentrations are quickly smoothed out, due to the diffusion regularization effect, or mathematically, the Laplacian term in Fick's law (Figure1.a). The mole number of each species, i.e. the area under each curve, is conserved over time. Since the concentration of Pu and U is imposed at both domain boundaries, the steady state is reached with the straight profiles shown in Figure1.b: the time derivative is equal to zero, and consequently the spatial second-order derivative of the concentration is also equal to zero. This computation allows to highlight the manifestations of diffusion phenomena at any points of the single-crystals, except at the ends points. The other observation is that the complete diffusion takes less than 4 hours: the interdiffusion totally occurs during the sintering. However, the values of concentration at the ends, which remain fixed during the interdiffusion, do not represent the physical reality of interdiffusion within crystallites. An alternative would be to use another type of boundary condition based on the zero flux of each element at the ends of the domain: the Neumann boundary conditions described above (paragraph 2.4).The numerical implementation of Neumann boundary conditions as described above allows us to plot the concentration profiles of U and Pu shown in Figure2 :

Figure 3 :

 3 Figure 3 : Position of the acquisition line in relation to the bonding interface separating the two half-couples UO2/U0.55Pu0.45O2

Figure 5 :

 5 a) During the first five minutes b)During the first hour c) During the first four hours

Figure 5

 5 Figure 5 : concentration profiles at different times

Figure 6 :

 6 Figure 6 : Evolution of standard deviation  on a ∆x observation scale as a function of time

Figure 7 :Figure 8 :

 78 Figure 7 : Evolution of standard deviation  on a 2∆x discretization scale as a function of time

Figure 9 :

 9 Figure 9 : Evolution of standard deviation  on a ∆x observation scale as a function of time

a)

  During the first ten minutes b) During the first four hours

Figure 10 :Figure 11 :

 1011 Figure 10 : Concentration profiles at different timeIn this case, achieving the homogenization seems to take longer than in previous cases due to the larger mean size. After four hours, the concentration profile is not yet constant as it was in the case of unimodal distribution. This representation can demonstrate a lot in the non-achievement of homogenization, but is not sufficient to analyze it quantitatively.On the other hand, the plot of concentration profile for the other of 45% Pu content is shown in Figure11. Its consists of four single-crystals of UO2 (300nm, 250nm, 200nm and 100nm) and four of Puo2 (250nm, 200nm, 150nm and 100nm).

5) 2.2. Physical model

	D
	e
	p
	a
	r
	t i n g	𝜕𝑐 𝑘 𝜕𝑡	= ∑ 𝑛-1 𝑗=1	𝐷 𝑘𝑗 𝑛 𝜕 2 𝑐 𝑗 𝜕𝑥 2	⇒ { 𝜕𝑐 𝑃𝑢 𝜕𝑐 𝑈 𝜕𝑡 𝜕𝑡	= 𝐷 𝑈𝑈 𝑂 𝜕𝑐 𝑈 𝜕𝑥 2 + 𝐷 𝑈𝑃𝑢 𝑂 = 𝐷 𝑃𝑢𝑈 𝑂 𝜕𝑐 𝑈 𝜕𝑥 2 + 𝐷 𝑃𝑢𝑃𝑢 𝜕𝑐 𝑃𝑢 𝜕𝑥 2 𝜕𝑥 2 𝑂 𝜕𝑐 𝑃𝑢	⇒	𝜕𝑐 𝜕𝑡	= 𝐷	𝜕 2 𝑐 𝜕𝑥 2	( 6)
	f r with D= (	𝐷 𝑈𝑈 𝑂 𝐷 𝑃𝑢𝑈 𝑂	𝐷 𝑈𝑃𝑢 𝑂 𝐷 𝑃𝑢𝑃𝑢 𝑂	)		and c = ( 𝑐 𝑃𝑢 𝑐 𝑈	) ; thus, the matrix D is diagonalized, using the eigenvector matrix B and the
	o eigenvalue diagonal matrix 𝜆. The relations between the matrices B, D, c, c' and 𝜆 are:
	m t				𝑐 ′ = 𝐵 -1 𝑐 𝐷𝐵 = 𝐵𝜆			⇒		𝜕𝑐 ′ 𝜕𝑡	= 𝜆	𝜕 2 𝑐 ′ 𝜕𝑥 2	(7)
	h e with 𝑐 ′ = ( 𝑐 𝑃𝑢 ′ 𝑐 𝑈 ′ ) . Let us express the components of matrices B and λ:
	e q B= (	𝑏 𝑈𝑈 𝑏 𝑈𝑃𝑢 𝑏 𝑃𝑢𝑈 𝑏 𝑃𝑢𝑃𝑢	) and λ = ( 𝜆 𝑈 0 𝜆 𝑃𝑢 0	)
	u a t with 𝑏 𝑈𝑈 = 1, 𝑏 𝑈𝑃𝑢 =	𝐷 𝑈𝑈 -𝜆 𝑈 𝐷 𝑃𝑢𝑈	, 𝑏 𝑃𝑢𝑈 =	𝐷 𝑃𝑢𝑈 𝐷 𝑈𝑈 -𝜆 𝑃𝑢	, bPuPu = 1,
	i o n r	𝜆 𝑈 = 𝜆 𝑃𝑢 = 1 2 1 (𝐷 𝑈𝑈 + 𝐷 𝑃𝑢𝑃𝑢 -√(𝐷 𝑈𝑈 -𝐷 𝑃𝑢𝑃𝑢 ) 2 + 4𝐷 𝑃𝑢𝑈 𝐷 𝑈𝑃𝑢 2 (𝐷 𝑈𝑈 + 𝐷 𝑃𝑢𝑃𝑢 + √(𝐷 𝑈𝑈 -𝐷 𝑃𝑢𝑃𝑢 ) 2 + 4𝐷 𝑃𝑢𝑈 𝐷 𝑈𝑃𝑢
	e
	s
	u
	l
	t
	i
	n
	g
	f
	r
	o
	m
	F
	i
	c
	k
	'
	s
	l
	a
	w
	s
	E
	q
	u
	a
	t
	i
	o
	n


G Hassan
email: ghina.hassan@lirmm.fr

A Chemori
email: ahmed.chemori@lirmm.fr

M Gouttefarde
email: marc.gouttefarde@lirmm.fr

M El

C Francis
email: cfrancis@ul.edu.lb

P E Hervé

D Sallé

A Novel Extended Desired Compensation Adaptive Law for High-Speed Pick-and-Throw with PKMs

Keywords: Robust DCAL, adaptive gains, parallel kinematic manipulator, pick-and-throw, numerical simulations

This paper focuses on the development of a new revised Desired Compensation Adaptive Law (DCAL). DCAL is a model-based adaptive control strategy consisting of three main parts: (i) an adaptive feedforward term, (ii) a linear PD feedback term, and (iii) a nonlinear compensation term. In order to deal with highly nonlinear dynamic systems characterized by their abundant uncertainties and parameters variations, we propose to revise the original DCAL control law by adopting adaptive feedback gains depending on the system state errors. Besides, DCAL controller is known for its robustness against measurement noise thanks to its desired compensation design, but a large amount of external disturbances are still not compensated by such a design. Therefore, the proposed DCAL with adaptive gains (DCAL-AG) is extended with a sliding-based term to further improve its robustness and the overall performance. A model-based robust adaptive feedback controller appropriate to the control of nonlinear systems in real-time applications is thereby obtained. To demonstrate the improvements brought by the proposed control strategy, numerical simulations have been conducted on a Delta-link parallel robot named T3KR in a "Pick-and-Throw" application task at different operating conditions.

INTRODUCTION

Over the last decades, parallel kinematic manipulators (PKMs) have provided a part of industrial and research resources. These manipulators are characterized by their rigidity, high precision, high dynamics and lightness [START_REF] Merlet | Parallel robots[END_REF]. Nevertheless, PKMs inherit the complexity of their closed kinematic chain structure where two or more kinematic chains connect the moving plate to the fixed base. In addition, they are characterized by high uncertainties, parameter variations, external disturbances [START_REF] Pi | Trajectory tracking control of a 6dof hydraulic parallel robot manipulator with uncertain load disturbances[END_REF] and nonlinear dynamics, especially in high-acceleration applications [START_REF] Natal | Dual-space control of extremely fast parallel manipulators: Payload changes and the 100g experiment[END_REF]. As a result, the design of sophisticated control strategies for these kind of manipulators is a challenging task. In the literature, several control schemes have been proposed aiming to accurately drive PKMs. On the one hand, the kinematic controllers, known by non-model-based controllers, can achieve acceptable performance as long as the operating conditions do not change [START_REF] Saied | A new time-varying feedback rise control for second-order nonlinear mimo systems: theory and experiments[END_REF]. Nevertheless, as mentioned above, PKMs are nonlinear dynamic systems, subject to uncertainties and time-varying parameters, thus, a kinematic controller may deteriorate the performance and lead to undesirable behavior. On the other hand, research works show that designing a controller that is partially or fully rich in knowledge about the system dynamics can improve tracking performance by compensating for the system nonlinearities [START_REF] Codourey | Dynamic modeling of parallel robots for computed-torque control implementation[END_REF][START_REF] Shang | Augmented nonlinear pd controller for a redundantly actuated parallel manipulator[END_REF]. However, these types of controllers require an accurate dynamic model of the robotic system which is a difficult task or even an impossible one. Therefore, the need for adaptive control schemes arises. Modelbased adaptive controllers can online adjust the dynamic parameters, leading to an adequate compensation of dynamic nonlinearities and possible parameter variations and uncertainties of PKMs [START_REF] Bennehar | A novel adaptive terminal sliding mode control for parallel manipulators: Design and real-time experiments[END_REF][START_REF] Bennehar | Extended model-based feedforward compensation in 1 adap-tive control for mechanical manipulators: Design and experiments[END_REF](Bennehar et al., , 2018)). The Desired Compensation Adaptive Law (DCAL), developed by [START_REF] Sadegh | Stability and robustness analysis of a class of adaptive controllers for robotic manipulators[END_REF][START_REF] Sadegh | Stability and robustness analysis of a class of adaptive controllers for robotic manipulators[END_REF], has been applied, to drive a six-degrees-of-freedom (6-DOF) PKM named Hexaglide, for the first time in [START_REF] Honegger | Adaptive control of the hexaglide, a 6 dof parallel manipulator[END_REF]. It has shown a good tracking performance, while estimating in realtime all the inertial and friction parameters of the manipulator. Both in the control and in the adaptation laws, this controller uses the desired trajectories instead of the measured ones, which can explain its effectiveness. An extended version of DCAL has been proposed in [START_REF] Bennehar | A new extension of desired compensation adaptive control and its real-time application to redundantly actuated pkms[END_REF] to enhance accuracy of PKMs. Thanks to their above advantages, a very wide range of applications benefit from PKMs. Recently, PKMs have been used as a robotic solution for selective waste sorting (BHS, 2018). Such an application is considered as a difficult task for PKMs, since the manipulator has to handle different types of objects with different physical parameters, that may often be unknown or uncertain. Therefore, model-based adaptive schemes, characterized by dynamic parameter identification in an online algorithm, are the most appropriate control solutions for such kind of applications. For instance, the aforementioned DCAL may be a good candidate, thanks to its simple structure easy to implement, its real-time estimation of the model parameters, and its robustness against measurement noise. Nevertheless, DCAL is characterized by constant linear feedback gains. As known, static feedback control algorithms can provide good performance only in nominal steady state and when no changes in operating conditions occur. However, it has been shown that control methods with adaptive dynamic feedback gains can counteract external disturbances and accommodate the variations in dynamic parameters [START_REF] Gholami | Adaptive cascade control of the kntu cdrpm: a cable driven redundant parallel manipulator[END_REF][START_REF] Tijjani | Robust adaptive tracking control of underwater vehicles: Design, stability analysis and experiments[END_REF][START_REF] Escorcia-Hernández | Contribution to nonlinear robust control of parallel kinematic manipulators: design and experiments[END_REF]. In this paper, we aim to deal with uncertain and unknown objects in a Pick-and-Throw application using a PKM. Accordingly, for this task, online estimation of dynamic parameters is required. Therefore, we propose to exploit the advantages of the real-time estimation of the model parameters provided by DCAL and the corrective action produced by an adequate adaptation law for the feedback gains. In addition, to better counteract the external disturbances, we propose to extend the resulting controller by a nonlinear sliding-based term computed from the signum of the system state errors. The addition of this robustness related term will accommodate for the lack of robustness and can improve the overall tracking performance. In a real-time implementation, the discontinuous signum function may be replaced by a continuous sigmoid function to avoid chattering. For the adaptation law of the feedback gains, we propose to use the algorithm developed by [START_REF] Plestan | New methodologies for adaptive sliding mode control[END_REF]. It is a continuous adaptation law that ensures a nonoverestimation of the gains with respect to the perturbations. Numerical simulations have been conducted in different scenarios of a Pick-and-Throw application in order to investigate the enhancement and the robustness brought by the proposed control scheme. The rest of the paper is organized as follows: The proposed control strategy is introduced in Section 2. The description and modeling of T3KR parallel robot are provided in Section 3. In Section 4, the obtained numerical simulation results in different scenarios are discussed, and finally, some concluding remarks are drawn in Section 5.

PROPOSED CONTRIBUTION: ROBUST DCAL WITH ADAPTIVE FEEDBACK GAINS

The dynamics of a m-DOF kinematic manipulator controlled by n actuators can be described in joint space as follows [START_REF] Siciliano | Robotics: modelling, planning and control[END_REF]:

M(q) q +C(q, q) q + G(q) + Γ d (t) = Γ(t) (1)
where M(q) ∈ R n×n is the total mass and inertia matrix of the robot, C(q, q) ∈ R n×n denotes the Coriolis and centrifugal forces matrix, G(q) ∈ R n is the gravitational forces vector. q, q, q ∈ R n are the joint position, velocity and acceleration vectors, respectively. The vector Γ d (t) ∈ R n gathers a large class of nonlinear disturbances (i.e. external disturbances, unknown friction effects, unmodeled phenomena, etc.) and Γ(t) ∈ R n is the control input vector.

According to [START_REF] Craig | Adaptive control of mechanical manipulators[END_REF], the manipulator dynamic model is characterized by its linearity with respect to dynamic parameters such as inertia and masses. All constant parameters in the dynamic model are considered as coefficients of known functions (linear and nonlinear) of q, q, q. The external disturbances Γ d (t) are excluded from the linear reformulation of the dynamics since they are not modeled and cannot be written in a linear form of the parameters. Therefore, (1) can be rewritten as follows:

M(q) q+C(q, q) q+G(q)+Γ d (t) = W(q, q, q)Φ(t)+Γ d (t) (2)

where W(q, q, q) ∈ R n×p is called the regression matrix and is formed by known nonlinear functions of q, q, q. The vector Φ ∈ R p gathers the geometrical and dynamic parameters of the robot. In the sequel, a background on the standard DCAL will be provided. Then, the proposed control approach will be detailed.

General overview of DCAL control strategy

DCAL is a model-based adaptive control scheme developed by [START_REF] Sadegh | Stability and robustness analysis of a class of adaptive controllers for robotic manipulators[END_REF][START_REF] Sadegh | Stability and robustness analysis of a class of adaptive controllers for robotic manipulators[END_REF]. Its control law can be split up into three main parts: (i) a modelbased adaptive feedforward part, (ii) a linear feedback part, and (iii) an additional nonlinear feedback function. The relevance of DCAL lies in both the control and adaptation laws which use the desired joint trajectories instead of the measured ones. This is of great importance since the computational time may be significantly reduced and the effect of measurement noises is eliminated. The additional nonlinear term aims to accommodate for the errors resulting from using the desired states instead of the measured ones. The joint-space control law of DCAL is then expressed as follows [START_REF] Sadegh | Stability and robustness analysis of a class of adaptive controllers for robotic manipulators[END_REF]:

Γ DCAL = W (q d , qd , qd ) Φ(t)+Λ p e(t)+Λ v e v (t)+σ ∥e(t)∥ 2 e v (t) (3
) where e(t) = q d (t)q(t) is the joint position tracking error, with q d (t) ∈ R n is the vector of desired joint positions and q(t) ∈ R n is the vector of measured ones. e v (t) = ė(t) + λ e(t) is the combined position-velocity tracking error, λ ∈ R + is a positive design gain. Λ P , Λ v ∈ R n×n are positive-definite gain matrices, usually chosen diagonal. W (q d , qd , qd ) ∈ R n×p is the regressor matrix function depending on desired joint positions, velocities and accelerations. Φ(t) ∈ R n×p is an online estimation of the unknown parameters vector Φ, and σ ∈ R + is a positive design control parameter. The time-evolution of the estimated parameters Φ(t) in ( 3) is expressed by the following adaptation law:

Φ(t) = KW T (q d , qd , qd )e v (t) (4)
where K ∈ R p×p is a diagonal positive-definite adaptation gain matrix. As it can be seen, the regressor W in the adaptation law (4) is also evaluated based on the desired trajectories instead of the measurements.

Proposed robust DCAL with adaptive gains

Despite the efficiency of the standard DCAL, it exhibits a lack of performance due to the static linear feedback gains and the potential presence of external disturbances not compensated by the control law. To significantly improve the overall performance of such a controller, we first propose to replace the linear feedback term with an adaptive one where the gains are adjusted online according to the system state errors. Second, to further improve its robustness against disturbances, a slidingbased term depending on the combined error is added. The resulting expression of the proposed control law can be written as follows:

Γ RDCAL-AG = W (q d , qd , qd ) Φ(t) + Λ p (t) e(t) + Λ v (t) e v (t) + σ ∥e(t)∥ 2 e v (t) + β sgn(e v (t)) (5)
where β ∈ R n×n is a constant positive-definite gain matrix. Λ p (t) and Λ v (t) ∈ R n×n are time-varying gains matrices. One

IFAC ALCOS 2022

June 29 -July 1, 2022. Casablanca, Morocco interesting choice of the time evolution of the adaptive feedback gain matrices can be formulated based on the mechanism proposed in [START_REF] Plestan | New methodologies for adaptive sliding mode control[END_REF] as follows:

Λ p (t) = Λp |η p | + Λ pm , ηp = tanh(e) -η p (6) Λ v (t) = Λv |η v | + Λ vm , ηv = tanh(e v ) -η v (11)
Where Λp and Λv ∈ R n×n are positive-definite constant matrices, chosen to be diagonal. While Λ pm and Λ vm ∈ R n×n are other positive-definite diagonal matrices denoting the minimum value for each adaptive gain. η p and η v ∈ R n are nonlinear functions depending on the tracking error e and the combined tracking error e v , respectively. For more details on this adaptation law and the stability analysis using it in a robust control law, the reader can refer to Escorcia-Hernández (2020). It is worth to note that when the tracking error increases, the adaptive gains of the proposed control law produce a corrective action to reduce this large tracking error. Once it decreases, the adopted strategy begins to reduce the control action and adjusts the gains to avoid oscillations and sufficiently counteract the current uncertainties and disturbances. When it comes to the estimation of unknown dynamic parameters, the same adaptation law ( 4) is adopted for the proposed controller. Therefore, this new control technique inherits the advantages of the original DCAL in terms of noise measurement reduction and low computational time.

T3KR ROBOT: DESCRIPTION, MODELLING AND CONTROL APPLICATION

3.1 Description and kinematics of T3KR robot T3KR, shown in Fig. 1, is a rigid-link parallel robot designed within the framework of a cooperation between SATT AxLR, Tecnalia and LIRMM. It features an economical footprint with five DOF. Three translational motions along x, y and z axes, and one rotational motion, ψ, of the moving platform around the vertical z axis, are generated by the four main actuators placed on the fixed base. These motors are connected to the moving platform by four kinematic chains. Each of the kinematic chains is formed by a main actuator, a movable reararm and a forearm composed of two parallel rods (cf. Fig. 1). In addition, a rotational movement, φ , of the robot end-effector around the z axis, is provided by a further actuator fixed at the mobile platform.

It should be noted that the ψ rotation of the platform is kept at zero for all the proposed scenarios. It is worth to emphasize the innovative point of this robot: if the last coordinate is changed, the tool control point (TCP) does not move; in fact the rotation of the platform is a parallelogram mechanism movement, and the TCP is on the neutral axis of the mechanism. In our study, we are concerned only with the control of the parallel Delta-like positioning structure. Consider the vector X = [x, y, z, ψ] T as the Cartesian position and orientation of the end-effector, and the vector q = [q 1 , q 2 , q 3 , q 4 ] T as the actuated joint positions. The differential kinematic relationship between the Cartesian and joint velocities can be expressed, using the Jacobian matrix J, as follows: Ẋ = J q, where Ẋ and q are the Cartesian and joint velocities, respectively.

Dynamics of T3KR parallel robot

Based on the virtual work principle, the dynamics of T3KR robot can be elaborated [START_REF] Codourey | Dynamic modeling of parallel robots for computed-torque control implementation[END_REF]. The following assumptions are considered to simplify the PKM dynamic model while maintaining its relevance: Assumption 1: The masses of the forearms are smaller than the others parts of the robot, hence their inertia is neglected.

Assumption 2: The mass of each forearm is divided into two pointwise masses located at both extremities of the forearms. Assumption 3: Both dry and viscous frictions in all passive and active joints are neglected.

The dynamics of T3KR robot can be reduced to the analysis of the dynamics of its moving platform and those of the actuators with their corresponding reararms and forearms. Regarding the moving platform's dynamics, one can consider two kinds of forces acting on it produced by the gravity and Cartesian accelerations. The contributions of these two forces to actuator torques can be expressed as follows:

Γ G t p = -J T M t p G, Γ F t p = J T M t p Ẍ (7)
where M t p = diag{m t p , m t p , m t p , i t p } with m t p = m n + 4 m f 2 is the total mass of the mobile platform including the mass of the actuator integrated on the platform, the payload handled by the end-effector and the four half masses of the forearms. i t p is the total inertia of the platform. G ∈ R 4 is the gravity vector (G = [0, 0, g, 0] T , being g = 9.81 m/s 2 the gravity acceleration). Ẍ ∈ R 4 denotes the Cartesian acceleration vector. Regarding the dynamics of the reararms, three torques acting on them can be distinguished : (i) the contribution of the actuators input torque Γ ∈ R 4 , (ii) the torque due to the gravitational forces acting on the reararms Γ G arm ∈ R 4 , and (iii) the inertial contribution torque produced by the joint acceleration on the reararms Γ arm ∈ R 4 :

Γ G arm = -g M r Cos(q), Γ arm = I arm q (8)

where

M r = diag{m req , m req , m req , m req } with m req = m r L rG + L m f 2
, m r is the mass of each reararm, l rG is the distance from the axis of rotation of each reararm to its center of gravity, while L is the complete length of each reararm. Cos(q) = [cos(q 1 ), cos(q 2 ), cos(q 3 ), cos(q 4 )] T and q ∈ R 4 represents the accelerations in joint space. I arm ∈ R 4x4 is a diagonal inertia matrix that gathers the actuators inertia, the reamarms inertia and the inertial contribution of the forearms with respect to the actuators' rotation axes using Assumption 2. Following [START_REF] Codourey | Dynamic modeling of parallel robots for computed-torque control implementation[END_REF], the sum of all non-inertial forces should be equal to the sum of all inertial forces, then the inverse dynamic model of T3KR robot can be expressed in terms of the joint coordinates as follows:

M(q) q +C(q, q) q + G(q) = Γ(t) (9) where M(q) = I arm + J T M t p J denotes the total mass and inertia matrix of the robot, C(q, q) q = J T M t p J denotes the Coriolis and centrifugal forces matrix, with J being the time derivative of J, G(q) = -Γ G arm -Γ G t p represents the gravitational forces vector, and Γ(t) is the control input vector. For more details on the development of PKM dynamic model, the reader can refer to Bennehar et al. (2018). If the external disturbances are considered, the dynamic model of T3KR robot can be rewritten as in (1). The main geometric and dynamic parameters of T3KR parallel robot are summarized in Table 1.

Control application

Our main objective is to use the proposed RDCAL-AG controller in a Pick-and-Throw selective sorting task. In such an application, the robot is subject to payload variations since the mass of the mobile platform (including the payload) may vary continuously depending on the object being handled. Accordingly, in our case study, the adaptation algorithm of the parameters estimation only accounts for the mass of the moving platform, including the payload, while taking advantage of the known dynamic parameters of the other robot parts. Therefore, the dynamic model of T3KR ( 9) with the consideration of the external disturbances can be reformulated according to (2) as follows: M(q) q +C(q, q) q + G(q) + Γ d (t) = W u (q, q, q)Φ u + I arm q + gM r Cos(q) + Γ d (t)

The partial regression matrix W u only accounts for the dynamics of the mobile platform and is given by:

W u = J T (J q + J q + G) = J T ( Ẍ + G) (11)
Φ u is the mass m t p of the mobile platform including the mass of the carried payload, Φ u = m t p . Therefore, the original DCAL and the proposed control law are reformulated as follows:

Γ DCAL = W u (q d , qd , qd ) Φu + I arm qd + gM r Cos(q d ) + Λ p e(t) + Λ v e v (t) + σ ∥e(t)∥ 2 e v (t) (12)
Γ RDCAL-AG = W u (q d , qd , qd ) Φu + I arm qd + gM r Cos(q d ) + Λ p (t) e(t) + Λ v (t) e v (t) + σ ∥e(t)∥ 2 e v (t) + β sgn(e v (t)) ( 13
)
where Φu is the online estimation of Φ u .

NUMERICAL SIMULATION RESULTS

In this section, the obtained simulation results of the proposed RDCAL-AG, the original DCAL and a model-based adaptive robust control from the literature are presented and discussed.

Implementation issues

Model-based adaptive RISE control: The RISE-based adaptive control, developed in (Bennehar et al., 2018), will be used for comparison purposes as it is a robust model-based adaptive controller. Its control law is expressed as follows: (

) 14
where e 2 (t) = ė(t) + α 1 e(t) is the combined error. α 1 , α 2 , K s and β ∈ R n×n are positive-definite, diagonal gain matrices, I ∈ R n×n is identity matrix, t 0 is the initial time and sgn(.) is the sign function. n is equal to 4 in our case study. The standard PID controller will not be considered for comparison in this work because its performance has been shown in the literature to be less good than robust and model-based controllers [START_REF] Natal | Dual-space control of extremely fast parallel manipulators: Payload changes and the 100g experiment[END_REF]; [START_REF] Hassan | Rise feedback control of cable-driven parallel robots: design and real-time experiments[END_REF].

Reference trajectories generation:

The P&T reference trajectories, sketched in Fig. 2, are generated using a 3 rd order polynomial S-curve motion profile. In these trajectories, the robot has to pick and throw eight objects, located at different positions, towards a target position, P f , located outside of its workspace. The robot moves from the homing position P 0 to the first pick position P 1 . After picking the 1 st object, the robot accelerates while moving along a straight line towards the release position P r 1 at which it throws the object towards the target position P f . Once released, the object follows its free-flight ballistic trajectory to P f while the robot decelerates back to pick the second object. The same cyclic movement is repeated for the second, the third, and the fourth objects, located at P 2 , P 3 , and P 4 , respectively. After throwing the fourth object, the robot moves to P 5 , to pick the fifth object. The same throw motion is performed for the last fourth objects located at P 5 , P 6 , P 7 , and P 8 , respectively. After throwing the last object, the robot moves back to P 0 .

Performance Evaluation criteria: To evaluate the effectiveness of the proposed controller, a frequently used performance index, the Root Mean Square Error (RMSE) criterion, is adopted in our study. The RMSE for Cartesian translational positions RMSET C and joint positions RMSE J are expressed respectively as follows:

RMSET C = ( 1 N N ∑ i=1 (e 2 x (i) + e 2 y (i) + e 2 z (i)) (15)
RMSE J = ( 1 N N ∑ i=1
(e 2 q 1 (i) + e 2 q 2 (i) + e 2 q 3 (i) + e 2 q 4 (i)) ( 16)

where e x , e y and e z are the Cartesian position tracking errors along the x, y and z axes, respectively. While e q 1 , e q 2 , e q 3 and e q 4 denote the joint position tracking errors, and N is the total number of samples.

IFAC ALCOS 2022

June 29 -July 1, 2022. Casablanca, Morocco

Table 2. Summary of the tuned feedback gains. Tuning the control gains: For the controllers tested on T3KR, the gains were adjusted by the trial-and-error method. For the proposed RDCAL-AG controller, we first set a high value for λ and minimum possible values for Λ pm and Λ vm . Then, the gains Λp and Λv were adjusted, either increasing or decreasing, until the best performance is obtained. The K gain, responsible for the parameters' estimation, is increased gradually till obtaining a good convergence of the mass of the platform. The σ gain is then increased to improve the overall performance while keeping the control input values away from saturation. Finally, the gain β , responsible for the robustness of the controller, is increased progressively in order to obtain better performance while maintaining low chattering input signals. The resulting gains values of the proposed controllers are summarized in Table 2. It is worth to note that the gain parameters are adjusted while respecting the actuators limits. In addition, for a fair comparison, the common parameters between the three tested controllers on T3KR are set to be the same.

Standard DCAL Adaptive RISE Proposed RDCAL-AG Λ p = 1135 σ = 4.5×10 6 α 1 = 150 β = 2.5 Λ pm = 1135 σ = 4.5 × 10 6 Λ v = 22.7 K = 200 α 2 = 0.3 K = 180 Λ p = 4×10 7 λ = 100 λ = 100 K s = 22 Λ vm = 22.7 K = 300 Λ v = 105 β = 2.5

Obtained simulations results

Numerical simulations have been conducted on T3KR robot in a P&T task to demonstrate the performance of the proposed controller. A comparison between the standard DCAL, the RISE-based adaptive control and the proposed controller has been established, in Matlab/Simulink environment with a fixed step solver equal to 0.4 ms, using the P&T reference trajectories illustrated in Fig. 2. Two main scenarios have been conducted in this demonstration: 1) scenario 1: robustness towards payload changes, 2) scenario 2: robustness towards speed variations. To be more realistic, white noise has been added to the output joint positions in both scenarios as well as 10% of uncertainty on I arm . Therefore, the robot dynamic model and the one used in the adaptive feedforward term of the controllers are not exactly the same.

Scenario 1 (S1) -Robustness towards payload changes: This scenario has been performed with 4.2 G as maximum acceleration. It is the smallest value that allows the robot to throw an object outside of its workspace. Eight different objects have been considered for this P&T task. The red lines in Fig. 2 correspond to the trajectory portions where the robot carries a payload, while the green lines represent the portions after the release point where the robot is moving without a payload. The 1 st and 5 th objects located at P 1 and P 5 , respectively, are of 50 g of mass, the 2 nd and 6 th objects located at P 2 and P 6 , respectively, have a mass of 100 g (i.e. ∆ mass = +100 % w.r.t the 1 st object), the 3 rd and 7 th objects at P 3 and P 7 , respectively, have a mass of 150 g (i.e. ∆ mass = +200 % w.r.t the 1 st object), while the 4 th and 8 th ones located at P 4 and P 8 , respectively, are of 200 g of mass (i.e. ∆ mass = +300 % w.r.t the 1 st object). The Cartesian tracking errors for all the controllers are plotted in Fig. 3. It is clearly shown that the proposed controller outperforms the other controllers especially for z-axis. This can validate the robustness of the proposed RDCAL-AG controller, towards the effect of gravity, compared to the others controllers.

The RMSE performance indices, in both Cartesian and joint spaces, are evaluated for all controllers and the obtained results are summarized in Table 3. These indices show a significant improvement of 44.4 % in Cartesian space and 46.3 % in joint space w.r.t to DCAL. While compared to ARISE, the indices show an improvement of 38.6 % and 18.8 % in Cartesian and joint spaces, respectively. The evolution of the estimated parameter mt p , initialized to zero, is displayed in Fig. 4. It is worth to note that this mass includes both the mass of the carried payload and the mobile platform. This explains why the adaptive mass increases or returns to its nominal value depending on whether the robot is carrying a payload. Fig. 5 illustrates the evolution of the adaptive gains, Λ p (t) and Λ v (t), versus time. One can observe, on the left side of Fig. 5, that the minimum value of Λ p (t) is 1135 as defined in the value of Λ pm ; similarly, on the right side of Fig. 5, the minimum value taken by Λ v (t) is the one as established in Λ vm . Besides, it can be seen from this figure that, Λ p (t) can reach values of up to 5000 while Λ v (t) reaches 23.5. The evolution of the control input torques is depicted in Fig. 6. The control signals show, for all controllers, a good and smooth behavior within the admissible limits of the actuators of the robot (the maximum torque generated by the actuators of the T3KR robot is 28.9 N.m).

As a result, this scenario demonstrates the superiority of the proposed controller over the standard DCAL and ARISE controllers. The RDCAL-AG control scheme is more robust towards variations in payload, thus, it is more suitable for P&T applications.

Scenario 2 (S2) -Robustness towards speed variations: In this scenario, the operating acceleration of the robot endeffector is increased up to 9 G following the same P&T reference trajectories described above. The objective of this scenario is to evaluate the performance of the proposed control solution for high-speed motions, where the nonlinear effects of the parallel manipulator increase substantially. In this scenario, the robot carries the same objects used in the previous scenario. The Cartesian tracking errors for all controllers are recorded and depicted in Fig. 7. The tracking errors of the three control schemes increased notably on all translational axes. Nevertheless, the proposed controller shows noticeably better performance, compared to the others two controllers. Besides, it can be noticed that the high nonlinearities and disturbances induced by the increase in the accelerations are considerably compensated by the proposed RDCAL-AG controller.

CONCLUSION AND FUTURE WORK

In this work, we proposed to amend the original DCAL with adaptive gains function of the system errors to counteract perturbations and uncertainties. In addition, the controller has been extended by a nonlinear sliding-based term to further improve its robustness against external disturbances. The standard DCAL, the RISE-based adaptive controller and the proposed robust DCAL with adaptive gains (RDCAL-AG) have been implemented and compared through numerical simulations on T3KR parallel robot. The obtained results clearly show the superiority of the proposed control approach compared to the others two controllers, in terms of tracking accuracy and robustness towards payload and velocity changes. Future directions may focus on extending this work with the stability analysis of the proposed controller as well as its validation in real-time experiments. Moreover, theoretical approaches will be considered to justify how the gain of the sliding-based term can be chosen to maximize robustness.

  Fig. 1. A CAD view of the T3KR parallel robot with its main components.

  Fig. 2. 3D-view of the P&T reference trajectories of the robot (the red and green lines) with the ballistic motions of the thrown objects in Cartesian space. Γ ARISE = W u (q d , qd , qd ) Φu + I arm qd + gM r Cos(q d ) + (K s + I)e 2 (t) -(K s + I)e 2 (t 0 ) +

Fig. 3 .

 3 Fig. 3. Scenario 1: Evolution of the Cartesian tracking errors versus time for the tested controllers.

Fig. 4 .

 4 Fig. 4. Scenario 1: Evolution of the estimated mass versus time for the tested controllers.AG over DCAL with a sliding-based term is 46%. Fig.8shows the evolution of the estimated parameter mt p with respect to time. The adaptation law adjusts the mass, from an initial zero value, similarly for the three controllers. We can clearly observe the oscillations induced by the changes of the payload at each pick and throw cycle.The evolution of the adaptive gains, Λ p (t) and Λ v (t), versus time is depicted in Fig.9. As it can be seen, the behavior of the adaptive gains Λ p (t) and Λ v (t) is slightly modified by the increase in the operating acceleration; Λ p (t) manages to reach values close to 6000, while Λ v (t) reaches values up to 24. The control input signals for the four motors of the robot, for all the controllers, are displayed in Fig.10. As we can see, the control signals of all the control schemes are within the allowable capacities of the motors. They are continuous and chatteringfree since the sign function in the proposed controller is replaced by a continuous sigmoid function. Furthermore, it is clear that the values of the generated torques increase with the operating acceleration.Besides, it can be noticed that the high nonlinearities and disturbances induced by the increase in the accelerations are considerably compensated by the proposed RDCAL-AG controller.

Fig. 5 .

 5 Fig. 5. Scenario 1: Evolution of the adaptive gains, Λ p (t) and Λ v (t), versus time for the proposed controller.

Fig. 6 .

 6 Fig. 6. Scenario 1: Evolution of the control input torques versus time for the tested controllers.

Fig. 7 .

 7 Fig. 7. Scenario 2: Evolution of the Cartesian tracking errors versus time for the tested controllers.

Fig. 8 .

 8 Fig. 8. Scenario 2: Evolution of the estimated mass versus time for the tested controllers.

Fig. 9 .

 9 Fig. 9. Scenario 2: Evolution of the adaptive gains, Λ p (t) and Λ v (t), versus time for the proposed controller.

Fig. 10 .

 10 Fig. 10. Scenario 2: Evolution of the control input torques versus time for the tested controllers.

Fig. 11 .

 11 Fig. 11. Scenario 2: Evolution of the Cartesian tracking errors versus time for DCAL, DCAL with sliding-based term and RDCAL-AG controllers.

Table 1 .

 1 Summary of the main geometric and dynamic parameters of T3KR robot.

	Parameter	Value	Parameter	Value
	Rear-arm length (L)	400 mm	Nacelle mass (m n )	5.68 Kg
	Forearm length (l)	900 mm Actuator inertia (I act ) 0.000969 Kg.m 2
	Rear-arm mass (m r )	3.28 Kg	Rear-arm inertia (I r )	0.173723 Kg.m 2
	Forearm mass (m f )	0.8 Kg

Table 3 .

 3 Control performance evaluation

		RMSET C [mm]	RMSE J [deg]
		S1	S2	S1	S2
	Standard DCAL	0.0981	0.1252	0.0136	0.0149
	Adaptive RISE	0.0888	0.0877	0.0090	0.0101
	Porposed RDCAL-AG	0.0545	0.0629	0.0073	0.0074
	Improvements w.r.t DCAL	44.4 %	49.7 %	46.3 %	50.3%
	Improvements w.r.t ARISE 38.6 %	28.3 %	18.8 %	26.7%

Table 3

 3

	summa-

Fig. fig:IsolateData.eps depicts the Cartesian tracking errors for DCAL, DCAL extended by a sliding-based term and RDCAL-AG.

It is obvious that adding a sliding-based term to DCAL improves the tracking performance and its robustness. This improvement is about 8%, while the improvement of RDCAL-

ACKNOWLEDGEMENTS

This work is supported by Tecnalia "Pick-and-Throw" research project. The corresponding author acknowledged the Lebanese University and the social foundation AZM and SAADE for financial support.