Stéphane Abide
email: stephane.abide@univ-perp.fr

Mikaël Barboteu
email: barboteu@univ-perp.fr

Soufiane Cherkaoui
email: soufiane.cherkaoui@univ-perp.fr

Serge Dumont
email: serge.dumont@unimes.fr

A Semi-Smooth Newton and Primal-Dual Active Set Method for Non-Smooth Contact Dynamics

Keywords: AMS Subject Classification : 70F40, 70-08, 70E55, 35Q70 Granular media, Unilateral constraint, Friction, Rigid body, Discrete Element Method, Non-Smooth Contact Dynamics, Semi-Smooth Newton method, Primal-Dual Active Set, Bi-potential method, Augmented Lagrangian method, Multi-rigid-body contact, Numerical simulations

Multi-rigid-body dynamic contact systems, in other words Non Smooth Contact Dynamics (NSCD) problems, generate some inherent difficulties to multivocal laws, which results in non-linearities and non-smoothness associated to frictional contact models. Recently, Primal-Dual Active Set strategies (PDAS) have emerged as a promising method for solving contact problems. These methods are based on the following principle: the frictional contact conditions are restated as non-linear complementary functions for which the solution is provided by the iterative semi-smooth Newton method. Based on these prerequisites, this contribution aims to provide a generalization of the NSCD-PDAS for dynamic frictional contact problems. Several numerical experiments are reported for algorithm validation purposes and also to assess the efficiency and performances of PDAS methods with respect to the Newton/Augmented Lagrangian and the Bi-Potential methods.

Introduction

Predicting snow avalanches, railway ballast track deformations or designing effective process for C02 capture requires taking into account local interactions modeling to assess macroscopic properties. This research issue is particularly topical because of its wide range of applications and the computational resources constantly increasing. As a matter of fact, multi-rigid-body dynamic contact field, as granular materials, illustrates some inherent difficulties to multivocal laws, which results in non-linearities and non-smoothness associated to frictional contact models. Numerical simulations of discrete materials, and specifically granular ones, was first initiated within the framework of the Discrete Element Method (DEM) [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Dubois | The contact dynamics method: A nonsmooth story[END_REF][START_REF] Radjai | Contact dynamics as a nonsmooth discrete element method[END_REF]. DEM relies on appropriate regularisation techniques for non-linear contact laws while using an explicit time scheme to ease their numerical treatments. Currently, this approach is still widely used in the engineering process design (see for example the MFIX open source software [START_REF] Garg | Documentation of open-source mfix-dem software for gas-solids flows[END_REF]). Some DEM weaknesses, reported in [START_REF] Dubois | The contact dynamics method: A nonsmooth story[END_REF], include the small time steps inherent in the explicit time scheme or the damping introduced for stability matters.

Moreau established the theoretical concepts to properly handle the non-smooth mechanics featured by non-differentiable laws in its common definition. Thus, he introduced the second-order sweeping process as the framework of mathematical analysis and numerical schemes designed for granular materials. The cornerstone is the formulation for the Signorini problem in terms of velocities and impulses, leading to the Moreau-Jean time-stepping approach [START_REF] Jj | Some numerical methods in multibody dynamics: application to granular materials[END_REF]. This numerical scheme is implicit and the energy conservation property holds unlike the DEM framework. The local conditions of contact between rigid bodies are ensured by means of a Non-Linear Gauss-Seidel (NLGS) developed by M. Jean and J. J. Moreau [START_REF] Jj | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF][START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Jj | Numerical aspects of the sweeping process[END_REF]. Also, The numerical aspects and certain algorithmic developments for the Non-Smooth Contact Dynamics (NSCD) have been proposes in [START_REF] Acary | Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics[END_REF]. Some alternatives such as conjugate gradient solvers have been proposed [START_REF] Renouf | Conjugate gradient type algorithms for frictional multi-contact problems: applications to granular materials[END_REF], with a specific analysis for parallel computing [START_REF] Renouf | A parallel version of the non smooth contact dynamics algorithm applied to the simulation of granular media[END_REF][START_REF] Visseq | High performance computing of discrete nonsmooth contact dynamics with domain decomposition[END_REF]. As pointed out by Dubois et al. [START_REF] Dubois | The contact dynamics method: A nonsmooth story[END_REF], this non-linear formulation is computationally expensive compared to DEM, especially for compressed granular materials, but larger time steps can be considered. Another point to consider in NSCD is the way to ensure the local contact condition. Usual approaches are based on the bi-potential or the augmented Lagrangian theory [START_REF] De Saxcé | New inequality and functional for contact with friction: the implicit standard material approachâ[END_REF][START_REF] Fortin | Numerical simulation of granular materials by an improved discrete element method[END_REF][START_REF] Feng | The bi-potential method applied to the modeling of dynamic problems with friction[END_REF]. Efficiency of such methods strongly depends on the penalty coefficients [START_REF] Fortin | Numerical simulation of granular materials by an improved discrete element method[END_REF][START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF][START_REF] Fortin | An improved discrete element method based on a variational formulation of the frictional contact law[END_REF]. Some improvements have been proposed to overcome this issue [START_REF] Dumont | On enhanced descent algorithms for solving frictional multicontact problems: application to the discrete element method[END_REF][START_REF] Joli | Uzawa and newton algorithms to solve frictional contact problems within the bi-potential framework[END_REF], but they are more time consuming.

Numerical simulations carried out for hyper-elastic bodies with frictional contact have led to specific approaches [START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF][START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and coulombâs friction[END_REF][START_REF] Acary | Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact[END_REF][START_REF] Chabrand | Various numerical methods for solving unilateral contact problems with friction[END_REF][START_REF] Chouly | An overview of recent results on Nitsche's method for contact problems[END_REF][START_REF] Khenous | Hybrid discretization of the signorini problem with Coulomb friction. theoretical aspects and comparison of some numerical solvers[END_REF][START_REF] Lebon | Contact problems with friction: models and simulations[END_REF][START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF][START_REF] Raous | Numerical methods for solving unilateral contact problem with friction[END_REF][START_REF] Sofonea | Mathematical models in contact mechanics[END_REF][START_REF] Wriggers | Computational contact mechanics[END_REF]. Recently, Primal-Dual Active Set strategies (PDAS) have emerged as a promising method for solving contact problems (see [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF][START_REF] Hintermüller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF][START_REF] Hüeber | A primal-dual active set algorithm for threedimensional contact problems with coulomb friction[END_REF][START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF]). These methods are based on the following principle: the frictional contact conditions are restated as non-linear complementary functions for which the solution is provided by the iterative semi-smooth Newton method [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF][START_REF] Hintermüller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF]. In practice, contact with Coulomb friction conditions can be formulated in terms of a fixed point problem related to a quasi-optimisation problem.

Therefore, the whole elastic body's problem, including the frictional contacts, consists of successive solutions of many elastic problems with simpler boundary conditions: Dirichlet, Neumann or Robin boundary conditions depending on the pseudo-transient status of the contact nodes [START_REF] Abide | Analysis of two active set type methods to solve unilateral contact problems[END_REF][START_REF] Abide | Inexact primal-dual active set method for solving elastodynamic frictional contact problems[END_REF]. This method is also known as the stabilized or Nitsche's methods [START_REF] Chouly | An overview of recent results on Nitsche's method for contact problems[END_REF][START_REF] Chouly | An adaptation of Nitsche's method to the tresca friction problem[END_REF][START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. space semi-discretization and time-marching schemes[END_REF].

According to our knowledge, few references on Active Set methods for solving muti-rigid-body dynamic contact problems are available. Based on a pseudo-rigid body assumption, Koziara and Bicanic [START_REF] Koziara | Semismooth Newton method for frictional contact between pseudorigid bodies[END_REF] modified a semi-smooth Newton technique to effectively deal with the frictional contact problem. Within a framework of optimization, Sharaf [START_REF] Im | An active set algorithm for a class of linear complementarity problems arising from rigid body dynamics[END_REF] proposes an Active Set method for solving positive definite and positive semi-definite Linear Complementary Problems (LCP). The effectiveness of its approach has been tested on the toy problem of static vertical rigid spheres assuming frictionless contacts. The resulting discrete model for the Newton-Euler dynamics equations with Signorini's conditions is formulated as a LCP. For multi-rigid-body dynamic contact problems, Bar-boteu and Dumont [START_REF] Barboteu | A primal-dual active set method for solving multi-rigid-body dynamic contact problems[END_REF] have developed a PDAS type method to address the local treatment of the contact conditions using the Non-Smooth Contact Dynamics (NSCD) formalism [START_REF] Jj | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF][START_REF] Jj | Numerical aspects of the sweeping process[END_REF][START_REF] Jj | Application of convex analysis to some problems of dry friction[END_REF][START_REF] Jj | Unilateral contact and dry friction in finite freedom dynamics[END_REF]. Comparisons with the well-established methods bi-potential and augmented Lagrangian outline the efficiency of the NSCD-PDAS methods for the considered validation tests [START_REF] Barboteu | A primal-dual active set method for solving multi-rigid-body dynamic contact problems[END_REF], dealing with frictionless model of muti-rigid-bodies contacts. The noted efficiency in [START_REF] Barboteu | A primal-dual active set method for solving multi-rigid-body dynamic contact problems[END_REF] must be confirmed for problems involving friction. This contribution aims to provide a generalization of the NSCD-PDAS for dynamic frictional contact problems and to assess its efficiency with respect to the bi-potential and augmented Lagrangian theories.

The article is organized as follows. First, the mathematical model for the frictional contact conditions as formulated in contact dynamics is stated in Sec. 2. The complementary functions for both contact and friction are formulated in Sec. 3. In addition, their generalized derivatives are derived to provide the core of the Semi-Smooth Newton method. Section 4 presents the NSCD framework, especially the equations of motion, the local-global mapping and the NLGS iterative solver. Based on these prerequisites, Sec. 3 introduces two versions of the PDAS to tackle the frictional contact for solid rigid bodies. Several numerical experiments are reported in Sec. 6 for algorithm validation purposes and also to assess the efficiency and performances of PDAS methods with respect to the Newton/Augmented Lagrangian and the Bi-Potential methods. Finally, some perspectives and future works are presented in Sec. 7.

Frictional Contact conditions for Contact Dynamics

Let us consider, for the rest of this paper, a dynamic collection of rigid bodies involved in several simultaneous contacts. We focus exclusively on unilateral contact and dry friction interactions. In order to ensure reliable dynamics of such systems, one have to take into account the energy conservation properties. Hence, the choice of the contact model in the context of Non-Smooth Contact Dynamics (NSCD) depends on whether the nature of the physics to deal with is a collection of rigid bodies or not. On the assumption that a multi-rigid body system is concerned, the contact conditions are formulated in terms of velocity, unlike deformable bodies (cf [START_REF] Abide | Analysis of two active set type methods to solve unilateral contact problems[END_REF][START_REF] Abide | Inexact primal-dual active set method for solving elastodynamic frictional contact problems[END_REF]), relating thus the impulse forces to the velocities. Then, the issue is to predict the velocities of the bodies and the impulse forces acting on the simultaneous multi-contacts.

Denoting by α a potential contact between two rigid particles of our system (α being a contact of the set of contact nodes S), we recall in detail the unilateral contact with friction within the framework of rigid bodies. For that purpose, and for the sake of conciseness, we denote by n the unit inner normal vector from a particle to another and τ the associated tangent vector. We use the notation u and p for the local displacement and the local impulse tensor at the contact point α, respectively (a dot superscript represents the time derivative with respect to the time variable t, e.g. u = ∂u/∂t). Also, we denote by un and uτ the normal and tangential components of velocity u given by un = u • n, uτ = uu n n. Finally, p n and p τ will represent the normal and the tangential impulse forces on the contact point, defined by

p n = (p • n)n and p τ = p -p n n.
First, as long as the normal distance u n between the two particles (corresponding to a gap) remains positive, there is no impulse force acting on each particle, and therefore, the normal impulse force p n remains zero. When the gap u n = 0, a repulsive normal impulse force p n is mobilized at the contact point and it can take indefinitely large values depending on the forces acting on each particle. Furthermore, as we consider an unilateral contact with dry friction, the discrete conditions on the contact point with friction define two complementary relations, the first one called Signorini's conditions [START_REF] Signorini | Sopra alcune questioni di elastostatica[END_REF], between the normal velocity un and the normal impulse force p n , and the second one called Coulomb's conditions [START_REF] Desplanques | Amontons-Coulomb friction laws, a review of the original manuscript[END_REF], relating the tangential component of velocity uτ and the tangential component of impulse force p τ . The inelastic shock law's conditions can be written as follows:

if g α > 0 then p α n = 0, (2.1) if g α ≤ 0 then uα n ≥ 0, (2.2)
p α n ≥ 0, (2.3)
uα n p α n = 0, (2.4)
where g α is the signed distance between the particles (the gap). Relations (2.1) -(2.4) lead to the conditions of a complete contact law formulated by Moreau (see [START_REF] Jj | Some numerical methods in multibody dynamics: application to granular materials[END_REF], [START_REF] Jj | Numerical aspects of the sweeping process[END_REF]). The last condition (2.4), is a complementary relation that ensures impenetrability between bodies (see [START_REF] Dubois | The contact dynamics method: A nonsmooth story[END_REF] for more details).

On the other hand, for energy conservation purposes, (2.4), also called persistency condition, means that normal contact reaction only appear during persistent contact, and has to be added in order to vanish the work of normal contact impulses at time t and preserve the amount of energy before and after the shock. J.-J. Moreau has proven that these conditions ensure impenetrability between bodies; see J.-J. Moreau viability lemma (cf [START_REF] Jj | Some numerical methods in multibody dynamics: application to granular materials[END_REF], [START_REF] Jj | Numerical aspects of the sweeping process[END_REF]).

The Coulomb's friction law reads as,

p α τ = 0. (2.5)        p α n = 0 =⇒ uα n ≥ 0, p α n > 0 and ||p α τ || < µp α n =⇒ uα τ = 0, p α n > 0 and |p α τ || = µp α n =⇒ ∃β ≥ 0, uα τ = β p α τ ||p α τ || . (2.6)
where µ is the friction coefficient.

The non-smooth motion in the event of a collision implies velocity jumps. These velocity jumps are the result of the velocity formulation of the non-inter-penetrability condition. The time derivative u is not unique and then, the left-limit velocity un and the right-limit velocity u+ n have to be distinguished. By analogy with a binary collision between two rigid particles, one can compute the right-limit velocity u+ n on the basis of the left-limit one just before the contact in a time-stepping scheme. Thus, the computation of un involved in Signorini's conditions is motivated by a simple choice which requires the use of a weighted average between un and u+ n : un = η un + (1η) u+ n , where η is a material parameter characterizing the contact (for more details about the computation of un , see [START_REF] Radjai | Contact dynamics as a nonsmooth discrete element method[END_REF]). Then, for η = 0, a binary shock between two rigid particles involves

-u+ n u- n = η (1-η)
. By identifying this ratio with normal restitution coefficient of the material, we obtain η = en (1+en) . Therefore, we define:

un = u+ n + e n u- n (1 + e n ) (2.7)
Identically, the computation of the tangential velocity un related to the frictional Coulomb's conditions follows the same way. The tangential impulse force p τ represents the average effect of static and repulsive forces relative to the contact during the time lapse ∆t. The tangential velocity is then an average velocity (see [START_REF] Radjai | Contact dynamics as a nonsmooth discrete element method[END_REF]). By the same token as normal impulse forces, a simple consistent model with the tangential restitution coefficient is to define uτ as follows:

uτ = u+ τ + e τ u- τ (1 + e τ ) ( 2
) = 0 C p n ( uα n , p α n ) = p α n -[p α n -γ n uα n ] + ∀α ∈ S, (3.1)
where S is the set of all potential contact nodes, and γ n the normal Active Set parameter. Now, let us prove this result.

Proposition 3.1. Let γ n > 0, the unilateral contact conditions expressed in velocity (2.2)-(2.4) for each potential contact α of the set of nodes S are equivalent to C p n ( uα n , p α n ) = 0, where p α n is the normal impulse force between two particles in contact.

Proof. Let us assume that (2.2)-(2.4) hold. We consider successively the cases uα n > 0 and uα n = 0. First, if uα n > 0, the condition uα n p α n = 0 implies that p α n = 0. Thus,

C p n ( uα n , p α n ) = -[-γ n uα n ] + = 0,
since γ n > 0. We suppose now that uα n = 0 and p α n > 0; therefore

C p n ( uα n , p α n ) = p α n -[p α n ] + = 0.
Conversely, we assume now that C p n ( uα n , p α n ) = 0 holds; it implies that p α n ≥ 0. Next, if p α n = 0, we have

C p n ( uα n , p α n ) = -[-γ n uα n ] + = 0, which leads to uα n ≥ 0, since γ n > 0. Finally, if p α n > 0, we have C p n ( uα n , p α n ) = p α n -[p α n -γ n uα n ] + = 0 =⇒ p α n = p α n -γ n uα n
and since γ n > 0, it implies that uα n = 0, which concludes the proof. A similar proof is available in [START_REF] Barboteu | A primal-dual active set method for solving multi-rigid-body dynamic contact problems[END_REF].

Complementary function for the friction

The frictional Coulomb's conditions (2.6) are represented by the following non-linear complementary function

C p τ ( uα n , uα τ , p α n , p α τ ) = 0 C p τ ( uα n , uα τ , p α n , p α τ ) = max (µp α n , ||p α τ -γ τ uα τ ||)p α τ -µp α n (p α τ -γ τ uα τ ) ∀α ∈ S, (3.2)
where S is the set of all potential contact nodes, and γ τ the tangential Active Set parameter. Now, let us prove this result.

Proposition 3.2. let γ τ > 0, the frictional contact conditions expressed in velocity (2.6) for each contact α of the set of nodes S are equivalent to C p τ ( uα n , uα τ , p α n , p α τ ) = 0, where p α τ is the tangential impulse force between two particles in contact.

Proof. Let us assume that (2.6) hold. We consider successively the cases p α n = 0 and p α n > 0. First, if p α n = 0, it implies that uα n ≥ 0, and thus,

C p τ ( uα n , uα τ , p α n , p α τ ) = ||p α τ -γ τ uα τ ||p α τ , Since there is no contact, p α τ = 0. Finally, C p τ ( uα n , uα τ , p α n , p α τ ) = 0.
We suppose now that p α n > 0 and ||p α τ || < µp α n ; it implies that uα τ = 0. Therefore, and since γ τ > 0, we have

C p τ ( uα n , uα τ , p α n , p α τ ) = max (µp α n , ||p α τ ||)p α τ -µp α n p α τ , Then, C p τ ( uα n , uα τ , p α n , p α τ ) = µp α n p α τ -µp α n p α τ = 0. We suppose now that p α n > 0, ||p α τ || = µp α n and uα τ = β p α τ ||p α τ || with β ≥ 0; therefore C p τ ( uα n , uα τ , p α n , p α τ ) = max (µp α n , ||p α τ -γ τ β p α τ ||µp α n || ||)p α τ -µp α n (p α τ -γ τ β p α τ ||µp α n || ) = 0.
Conversely, we assume now that C p τ ( uα n , uα τ , p α n , p α τ ) = 0 holds; depending on the value of p α τ and uα τ , one can have

µp α n = max (µp α n , ||p α τ -γ τ uα τ ||), (3.3
)

||p α τ -γ τ uα τ || = max (µp α n , ||p α τ -γ τ uα τ ||), (3.4)
By combining (3.2) and (3.3), we get

µp α n p α τ -µp α n p α τ + µγ τ p α n uα τ = 0, which means that p α n uα τ = 0, since γ n > 0 and µ > 0. If p α n = 0, condition (2.4) implies that uα n ≥ 0. If not, p α n > 0, uα τ = 0 and from (3.3), µp α n > ||p α τ ||.
At last, we combine (3.2) and (3.4) to obtain

||p α τ -γ τ uα τ ||p α τ -µp α n p α τ + µγ τ p α n uα τ = 0. It is trivial that p α τ = µγ τ p α n µp α n -||p α τ -γ τ uα τ || uα τ , let β = µγτ p α n µp α n -||p α τ -γτ uα τ || . If p α n > 0 and ||p α τ || = µp α n
, and implies that uα n > 0, since γ n > 0, which concludes the proof.

Generalized derivative of complementary functions

Now, we provide the generalized derivative of the complementary functions in the gap, stick and slip cases.

• Gap case : p α n -γ n uα n ≤ 0 According to the complementary functions C p n ( uα n , p α n ) = p α n and C p τ ( uα n , uα τ , p α n , p α τ ) = p α τ -γ τ uα τ p α τ , we have the following derivatives d uα n C p n = 0, (3.5)
d p α n C p n = dp α n , (3.6)
d uα n C p τ = 0, (3.7)
d uα τ C p τ = -γ τ p α τ (p α τ -γ τ uα τ ) T p α τ -γ τ uα τ d uα τ = 0, (3.8)
d p α n C p τ = 0, (3.9)
d p α τ C p τ = p α τ (p α τ -γ τ uα τ ) T p α τ -γ τ uα τ + p α τ -γ τ uα τ I 2 dp α τ . (3.10) • Stick case : µp α n ≥ p α τ -γ τ uα τ > 0 Given the complementary functions C p n ( uα n , p α n ) = γ n uα n and C p τ ( uα n , uα τ , p α n , p α τ ) = µγ τ p α n uα τ , we have d uα n C p n = γ n d uα n , (3.11)
d p α n C p n = 0, (3.12
)

d uα n C p τ = 0, (3.13)
d uα τ C p τ = µγ τ p α n d uα τ , (3.14)
d p α n C p τ = µγ τ uα τ dp α n , (3.15)
d p α τ C p τ = 0. (3.16) • Slip case : p α τ -γ τ uα τ > µp α n > 0 From C p n ( uα n , p α n ) = γ n uα n and C p τ ( uα n , uα τ , p α n , p α τ ) = p α τ -γ τ uα τ p α τ -µp α n (p α τ -γ τ uα τ ), it comes d uα n C p n = γ n d uα n , (3.17)
d p α n C p n = 0, (3.18)
d uα n C p τ = 0, (3.19)
d uα τ C p τ = -γ τ p α τ (p α τ -γ τ uα τ ) T p α τ -γ τ uα τ + µγ τ p α n I 2 d uα τ , (3.20)
d p α n C p τ = -µ(p α τ -γ τ uα τ )dp α n , (3.21)
d p α τ C p τ = p α τ (p α τ -γ τ uα τ ) T p α τ -γ τ uα τ + p α τ -γ τ uα τ I 2 -µp α n I 2 dp α τ . (3.22)
D C p n ( uα n , p α n )(δ uα n , δp α n ) = γ n (1 Stick + 1 Slip )δ uα n + 1 Gap δp α n , (3.23)
D C p τ ( uα n , uα τ , p α n , p α τ )(δ uα n , δ uα τ , δp α n , δp α τ ) = 1 Gap p α τ -γ τ uα τ δp α τ (3.24) + 1 Stick µγ τ p α n δ uα τ + µγ τ uα τ δp α n + 1 Slip -γ τ p α τ (p α τ -γ τ uα τ ) T p α τ -γ τ uα τ + µγ τ p α n I 2 δ uα τ -µ(p α τ -γ τ uα τ )δp α n + p α τ (p α τ -γ τ uα τ ) T p α τ -γ τ uα τ + p α τ -γ τ uα τ I 2 -µp α n I 2 δp α τ where 1 Gap = 1, 1 Stick = 0, 1 Slip = 0 if p α n -γ n uα n ≤ 0, 1 Gap = 0, 1 Stick = 1, 1 Slip = 0 if µp α n ≥ p α τ -γ τ uα τ > 0, 1 Gap = 0, 1 Stick = 0, 1 Slip = 1 if p α τ -γ τ uα τ > µp α n > 0.
Using now the semi-smooth Newton formalism at the current ( uα,(k)

n , uα,(k) τ , p α,(k) n , p α,(k) τ
), one can derive the new iterate ( uα,(k+1) ).

n , uα,(k+1) τ , p α,(k+1) n , p α,(k+1) τ ) D C p n ( uα,(k) n , p α,(k) n )(δ uα,(k+1) n , δp α,(k+1) n ) = -C p n ( uα,(k) n , p α,(k) n ), (3.25)
D C p τ ( uα,(k) n , uα,(k) τ , p α,(k) n , p α,(k) τ )(δ uα,(k+1) n , δ uα,(k+1) τ , δp α,(k+1) n , δp α,(k+1) τ ) (3.26) = -C p τ ( uα,(k) n , uα,(k) τ , p α,(k) n , p α,(k) τ ), ( uα,(k+1)
• Gap case:

1 Gap = 1, 1 Stick = 0, 1 Slip = 0
From the equations (3.25) and (3.26) we have

p α,(k+1) n -p α,(k) n = -p α,(k) n , (3.28)
p α,(k) τ -γ τ uα,(k) τ (p α,(k+1) τ -p α,(k) τ ) = -p α,(k) τ -γ τ uα,(k) τ p α,(k) τ . (3.29)
Next, the gap conditions of the semi-smooth Newton formalism are as follows

p α,(k+1) n = 0, (3.30) p α,(k+1) τ = 0, (3.31) since p α,(k) τ -γ τ uα,(k) τ > 0.
• Stick case :

1 Gap = 0, 1 Stick = 1, 1 Slip = 0
From the equations (3.25) and (3.26) we have

γ n ( uα,(k+1) n -uα,(k) n ) = -γ n uα,(k) n , (3.32
)

µγ τ p α,(k) n ( uα,(k+1) τ -uα,(k) τ ) + µγ τ uα,(k) τ (p α,(k+1) n -p α,(k) n ) = -µγ τ p α,(k) n uα,(k) τ . (3.33) Next, uα,(k+1) n = 0, (3.34) uα,(k+1) τ -uα,(k) τ
= -uα,(k) τ p α, (k+1) n p α,(k) n . (3.35)
For a given contact node α, the fundamental principle of dynamics can be written as follows:

uα = uα,free + W αα p α + β =α W βα p β , (3.36)
then, from (3.27), we have

   δ uα,(k+1) = uα,(k+1) -uα,(k) , δp α,(k+1) = p α,(k+1) -p α,(k) , (3.37)
using (3.36), one can write

δ uα,(k+1) = W αα δp α,(k+1) , (3.38)
and more particularly, • Slip case:

δ uα,(k+1) τ = W αα τ τ δp α,(k+1) τ , (3.39
1 Gap = 0, 1 Stick = 0, 1 Slip = 1
For C p n , we get once again uα,(k+1)

n = 0. (3.44)
For C p τ , we obtain

-γ τ p α,(k) τ (p α,(k) τ -γ τ uα,(k) τ ) T p α,(k) τ -γ τ uα,(k) τ + µγ τ p α,(k) n I 2 ( uα,(k+1) τ -uα,(k) τ ) (3.45) -µ(p α,(k) τ -γ τ uα,(k) τ )(p α,(k+1) n -p α,(k) n ) + p α,(k) τ (p α,(k) τ -γ τ uα,(k) τ ) T p α,(k) τ -γ τ uα,(k) τ + p α,(k) τ -γ τ uα,(k) τ I 2 -µp α,(k) n I 2 (p α,(k+1) τ -p α,(k) τ ) = -p α,(k) τ -γ τ uα,(k) τ p α,(k) τ + µp α,(k) n (p α,(k) τ -γ τ uα, (k) τ ).
For recall,

F (k) = p α,(k) τ (p α,(k) τ -γ τ uα,(k) τ ) T p α,(k) τ -γ τ uα,(k) τ , E (k) = 1 p α,(k) τ -γ τ uα,(k) τ .
Therefore, after an elementary computation

-γ τ E (k) F (k) -µp α,(k) n I 2 ( uα,(k+1) τ -uα,(k) τ ) -µE (k) (p α,(k) τ -γ τ uα,(k) τ )p α,(k+1) n + E (k) (F (k) -µp α,(k) n I 2 ) + I 2 p α,(k+1) τ -E (k) F (k) -µp α,(k) n I 2 p α,(k) τ = 0. Now, let M * (k) α = E (k) (F (k) -µp α,(k) n I 2 ), h (k) α = E (k) F (k) p α,(k) τ -γ τ uα,(k) τ . Then -γ τ M * (k) α uα,(k+1) τ -µE (k) (p α,(k) τ -γ τ uα,(k) τ )p α,(k+1) n + I 2 + M * (k) α p α,(k+1) τ = h (k) α -µE (k) (p α,(k) τ -γ τ uα,(k) τ )p α,(k) n .
In order to simplify even further the notations, we introduce the following operators:

L * (k) α = -γ τ (I 2 + M * (k) α ) -1 M * (k) α , r * (k) α = (I 2 + M * (k) α ) -1 h (k) α , v (k) α = µ(I 2 + M * (k) α ) -1 E (k) (p α,(k) τ -γ τ uα,(k) τ ).
And, at last

L * (k) α uα,(k+1) τ -v (k) α p α,(k+1) n + p α,(k+1) τ = r * (k) α -v (k) α p α,(k) n .
In this specific problem, and for a two dimensional case, one can obtain a simplified equivalent version of the algorithm. Let D slip C p τ be the generalized derivative of C p τ in the slip case

D slip C p τ ( uα τ , p α τ )(δ uα τ , δp α τ ) = p α τ (p α τ -γ τ uα τ ) T p α τ -γ τ uα τ (δp α τ -δ uα τ ) (3.46) -µp α n (δp α τ -γ τ uα τ ) + p α τ -γ τ uα τ δp α τ .
Denoting by τ , the unit slip vector, we have

p α τ = µp α n τ , (3.47)
(p α τ -γ τ uα τ ) p α τ -γ τ uα τ = τ , (3.48
)

δp α τ -γ τ uα τ = ητ . ( 3
D slip C p τ ( uα τ , p α τ )(δ uα τ , δp α τ ) = µp α n η(τ τ T -I 2 )τ + p α τ -γ τ uα τ δp α τ . (3.50)
Since τ τ T + nn T = I 2 in the 2D case, we have (τ τ T -I 2 )τ = nn T τ = 0. Using (3.26), we obtain

p α,(k) τ -γ τ uα,(k) τ (p α,(k+1) τ -p α,(k) τ ) (3.51) = -p α,(k) τ -γ τ uα,(k) τ p α,(k) τ + µp α,(k) n (p α,(k) τ -γ τ uα, (k) τ )
p α,(k+1) τ = µp α,(k) n (p α,(k) τ -γ τ uα,(k) τ ) p α,(k) τ -γ τ uα, (k) τ .
(3.54)

Non Smooth Contact Dynamics (NSCD)

As seen previously, the unilateral contact between rigid bodies involves a non-smoothness in law between the impulse force and the local relative velocity due to the shock because of friction, and a temporal non-smoothness because of the velocity jumps before and after shock. The Non-Smooth Contact Dynamics (NSCD) approach, based on two main computational level tasks, enables to solve the full problem. Indeed, the global level is devoted to solve the equations of motion, whereas the local level treats the contacts separately. Therefore, it is possible to solve, over a time step, many simultaneous contacts. The main purpose of this section, after presenting the equations of motion, is to come up with a general algorithm for NSCD for which the unilateral contact conditions with dry friction seen in (2.2)-(2.6) are treated numerically by several numerical methods.

Equations of motion

To describe the motion of a multi-contact system between rigid bodies, we use the following writing conventions. Assuming that a particle P among N p particles is described by the position of its center of gravity, we will denote by q the generalized coordinate describing its position in space, (q ∈ R d×Np , where d = 6 for a 3D problem and d = 3 for a 2D problem). As a consequence of the possible shocks between particles, we introduce the generalized velocity denoted by q as a function of bounded variations, and its associated differential d q. According to the fundamental principle of dynamics, the equations of motion formulated in terms of differential measures can be written as follows:

Md q + F int (t, q, q)dt = F ext (t, q, q)dt + dR (4.1) where

• M represents the genralized mass matrix;

• F int et F ext represent the internal and external forces respectively;

• dR is a non-negative real measure, representing the reaction forces and impulses between particles in contact.

For the sake of simplicity and without lost of generality, only the external forces are considered in the following.

Considering that the time interval of interest is [0, T ], we discretize for numerical purposes the previous equation (4.1). For that, we introduce uniform time instants t k defined by t k+1 = t k + ∆t pour k = 0, ..., N T -1, where

• ∆t = T /N T is the time step;

• N T is the number of time steps.

The equation (4.1) is then integrated on each time interval [t k , t k+1 ] and approximated using a θscheme, with θ ∈ [ 1 2 , 1] for stability reasons (see [START_REF] Jj | Numerical aspects of the sweeping process[END_REF], [START_REF] Renouf | Conjugate gradient type algorithms for frictional multi-contact problems: applications to granular materials[END_REF]). Therefore, the classical approximation of equation (4.1) yields

   M( qk+1 -qk ) = ∆t(θF k+1 + (1 -θ)F k ) + P k+1 , q k+1 = q k + ∆tθ qk+1 + ∆t(1 -θ) qk . (4.2)
where

• P k+1 represents the value of the total impulsion over the time step;

• F k (resp. F k+1 ) is the external force computed at time t k (resp. t k+1 ).

Denoting by qfree

k = qk + M -1 ∆t(θF k+1 + (1 -θ)F k )
, the velocity when the contact impulses vanish, the first equation of (4.2) becomes

qk+1 = qfree k + M -1 P k+1 (4.3)

Local-global mapping

In Non-Smooth Contacts Dynamics NSCD, the contact impulses are not explicit functions which define the state of equilibrium of the studied system. Therefore, contact impulses and velocities must be determined at the same time. As the contact law is expressed using the contact variables, we have to express the equations of motion using these same variables. Thus, a local-global mapping is defined from the global frame to the local one for each contact node α ∈ [1, N α ] (where N α is the total number of potential contacts between two rigid particles) to write the local contact law related. Local relative velocities u α n and u α τ between the two rigid bodies involved in the contact α can be grouped in a column vector u ∈ R d×Nc (where d = 3 for a 3D problem and d = 2 for a 2D problem). Likewise, the local contact impulses p α are represented by a vector p ∈ R d×Nc . We recall that p α can be decomposed into the sum of a normal component p α n and a tangential component p α τ as follows: p α = p α n n + p α τ . From there, the idea would be to transform the equations of motion seen in (4.2) from P and q to p and u.

For the sake of simplicity and to lighten the notations, the equations of motion (4.2) are converted into a single matrix equation whose expression is as follows:

M( q+ -q-) = ∆t(F + F ext ) (4.4)
The contact and particles velocities ( u and q) being linearly related, the local-global mapping to compute u to q at a contact node α is:

uα = H * (q, α) q (4.5)
where H * (q, α) local-global mapping matrix dN c × dN p carrying the information on the geometry of the contact network. Identically, we have

P = H(q, α)p α (4.6)
where H(q, α) is the transpose of matrix H * (q, α). Since the local-global mapping is computed for each contact node α, the total local-global mapping H(q, α) allows to compute all velocities and contact impulses from the global to the local frame:

   u = H * (q) q, P = H(q)p. (4.7)
By combining (4.2) and (4.7), the discretization of the multi-contact system's motion, with contact between rigid bodies can be written:

           ũk+1 = ũfree k + Wp k+1 , contact law( ũα n,k+1 , p α n,k+1 ) = .true. ∀α ∈ [1, ..., N α ], f riction law( ũα k+1 , p α k+1 ) = .true. ∀α ∈ [1, ..., N α ].
(4.8)

General algorithm for NSCD: Non-Linear Gauss Seidel Method (NLGS)

This paragraph is devoted to the detailed description of the algorithm used at the global level to solve the multi-rigid-body dynamic contact problems ( [START_REF] Fortin | Numerical simulation of granular materials by an improved discrete element method[END_REF]). Following the ideas of Jean and Moreau (see for example [START_REF] Jean | The non-smooth contact dynamics method[END_REF], [START_REF] Jourdan | A Gauss-Seidel like algorithm to solve frictional contact problems[END_REF] , [START_REF] Jj | Unilateral contact and dry friction in finite freedom dynamics[END_REF]), we use the NLGS algorithm which consists in considering successively each contact until the convergence. The time-stepping method combined with the NLGS algorithm takes the following form:

• Loop on the step time k -Prediction of a position (for the computation of the local-global mapping): • Actualization of the generalized displacement:

q k+ 1 2 = q k + ∆t 2 qk ; (4.
qj+1 k+1 = qj k+1 + M -1 P (q k+ 1 2 , α)p α,j+1
. * End of the loop on contacts α.

-End of the loop on j of NLGS. When the convergence is reached, actualization of the velocity: qk+1 = qj+1 k+1 -Actualization of the generalized displacements: q k+1 = q k+ 1 2 + ∆t 2 qk+1 • End of the loop on the step time k. Remark 1. The numerical solutions of both "Exact" and "Iterative" methods end the NLGS loop at each time step when convergence is reached. In other words, when the contact impulses are computed on each contact node at the current time step, a stopping criterion is used to exit the non-linear loop. This criterion can be stated as follows:

max α |p α,j+1 n | + max α |p τ α,j+1 | < N LGS (4.15)
where:

• p α,j+1
n is the normal component of the contact impulse computed at the end of the NLGS loop,

• p τ α,j+1 is the tangential component of the contact impulse computed at the end of the NLGS loop.

A Primal-Dual Active Set method for NSCD

This section is devoted to the numerical treatment of the contact conditions using a Primal-Dual Active Set method within the framework of multi-rigid-body dynamic contact problems. After defining the active and inactive subsets of all nodes that are currently in contact, we compute the contact conditions on each subset only in terms of contact impulses, using the local general equations of motion (4.8) in the form of (3.1) and (3.2).

"Exact" Primal-Dual Active Set method (EPDAS)

Let us denote by S the set of potential contact between particles of the granular material and α a potential contact node between two particles belonging to S. The frictional contact condition (2.2)-(2.6) are realized by applying an active set strategy which derives directly from the computation of the fixed point on the non-linear complementary functions C p n and C p τ based on the Newton semi-smooth scheme derived in part 3. The active and inactive sets are defined as follows

A j+1 n = {α ∈ S : p α,j n -γ n uα,j n ≥ 0}, I j+1 n = {α ∈ S : p α,j n -γ n uα,j n < 0}, A j+1 τ = {α ∈ S : p α,j τ -γ τ uα,j τ -µp α,j n ≥ 0}, I j+1 τ = {α ∈ S : p α,j τ -γ τ uα,j τ -µp α,j n ≤ 0}.
The status of a given potential α at the non-linear iteration j depends on the set it belongs to. According to part 3, it can be either in the non-contact, slip or stick status.

Hereafter, the numerical computation of the local contact step inside the NLGS iteration loop of index j leads to the following Primal-Dual Active Set algorithm.

(i) Choose ( u(0) , p (0) ), γ n > 0, γ τ > 0 and set j = 0.

(ii) Set the active and inactive sets:

A j+1 n = {α ∈ S : p α,j n -γ n uα,j n ≥ 0}, I j+1 n = S \ A j+1 n , A j+1 τ = {α ∈ S : p α,j τ -γ τ uα,j τ -µp α,j n ≥ 0}, I j+1 τ = S \ A j+1 τ .
(iii) Find ( u(j+1) , p (j+1) ) such that

"Iterative" Primal-Dual Active Set method (IPDAS)

Since we have noticed that (5.4) can lead to numerical instabilities when p α,j n is small, we introduce a variant of the previous method. Then, the condition (5.4) leads to p α,j+1 τ + uα,j τ Wττ = p α,j τ . Hereafter, the numerical computation of the local contact step inside the NLGS iteration loop of index j leads to the following Primal-Dual Active Set algorithm.

(i) Choose ( u(0) , p (0) ), γ n > 0, γ τ > 0 and set j = 0.

(ii) Set the active and inactive sets:

A j+1 n = {α ∈ S : p α,j n -γ n uα,j n ≥ 0}, I j+1 n = S \ A j+1 n , A j+1 τ = {α ∈ S : p α,j τ -γ τ uα,j τ -µp α,j n ≥ 0}, I j+1 τ = S \ A j+1 τ .
(iii) Find ( uj+1 , p j+1 ) such that (5.8)

(iv) If ( uj+1 , p j+1 ) -( uj , p j ) ≤ , A j+1 n = A j n and A j+1 τ = A j τ stop, else goto (ii).
Before concluding this section, we may underline one of the most characteristic feature of the active set type method, which is to enforce directly and exactly the contact conditions that are found.

Remark 2. In the simulations which follow, we will use two well-known and powerful methods to calculate the local contact impulses in the previous algorithm. These are methods based on the work of Alart and Curnier (Standard Augmented Lagrangian (SAL)), Saxcé and Feng (The Standard Bi-Potential (SBP)) and an improvement of the Bi-Potentiel Method (IBP, see [START_REF] Feng | The bi-potential method applied to the modeling of dynamic problems with friction[END_REF]), which will allow us to evaluate and compare with the Active Set method. For more detail about these methods, we refer to [START_REF] Dumont | On enhanced descent algorithms for solving frictional multicontact problems: application to the discrete element method[END_REF].

Remark 3. It should also be noted that only one iteration of the Active set method is computed for each contact at each global iteration of the NLGS loop, as in the case of the IBP method, and contrary to the other methods.

Numerical experiments

The main purpose of this section is to provide several numerical experiments on academic cases to assess the PDAS methods in solving granular media flows. To do so, we consider for the first two simulations the following reference configurations: a single spherical steel ball sliding and rolling on a conveyor belt and another one inside a fixed drum. These two configurations are addressed to outline the ability of the PDAS in enforcing basic mechanical properties as conservation properties, and, since the analytical solution is available for the first case, to compare it with the numerical one. The second part of this section brings the light on the robustness and accuracy of the PDAS methods, by carrying out more complex simulations like the sedimentation of a collection of balls. The computational performances of such method are illustrated by comparing it with two other methods (SAL, IBP). The final numerical experiment focuses on the granular flow in a 2D rotating drum which allows comparisons with experiments.

One sliding and rolling steel ball on a conveyor belt

As mentioned before, the aim of this first example is to carry out comparisons between numerical and analytical solutions. A spherical steel ball, with a non-vanishing initial horizontal velocity (v 0 = 1m.s -1 ) is placed on a conveyor belt with a constant velocity (V 1 = v 1 x with v 1 = 2m.s -1 ) as shown in Figure . We provide below a description of the physical settings ρ = 8000 Kg/m 3 , r = 2.7 10 -3 m, x 0 = 0 m, y 0 = 2.7 10 -3 m, g = -9.80665 m/s 2 Since there is friction, the steel ball slides until the sliding velocity is zero, then, starts rolling, rolling time t being the reference time from which the status of the punctual contact between the ball and the conveyor belt changes. Therefore, we go from sliding to rolling when

t = |v 0 -v 1 | µg(1 + r 2 m I ) (6.1)
with

I = 2mr 2 5
the moment of inertia of the steel ball. The duration of the simulation is T = 1s, and the time step is taken as a sub-multiple of t. For the numerical experiment, the parameter values are as follows: This simulation is computed with a PDAS method. The numerical solution obtained is then compared to the analytical one in terms of position, angular velocity and rolling velocity. Graphs in Figure 2 clearly show that numerical and analytical curves are the same, both in rolling and sliding steps. Moreover, one can notice (see Table 1) that the absolute error between the numerical and analytical solution for position angular velocity and rolling velocity is is zero, which allows us to validate our model in this case. Concerning the computing time (see Table 2), PDAS, SAL and IBP methods need a similar total CPU time to compute the solution. Moreover, the PDAS methods need only 2 NLGS iterations to converge during the sliding step, 3 for the rolling step, which means that the good status for the contact is directly found. For the SAL method, 20 iterations of NLGS are necessary, and finally, IBP method needs 2 NLGS iterations in both sliding and rolling steps to converge. Let us consider in this part another representative example, the same steel ball as the previous one, posed inside a fixed drum (see Figure 3), with a non-vanishing initial horizontal velocity (v 0 = 0.5m.s -1 ). We provide below a description of the physical settings: The dynamic behavior of the steel ball is similar to a pendulum. Indeed, the steel ball reaches a maximum position, position for which the ball's translational velocity is zero. Moreover, since there is friction, the ball has an angular velocity that varies over time. For the numerical experiment, the parameter values are as follows: Figure 4a describes the evolution of the ball's angular velocity over time during the sliding-rolling process inside the fixed drum. Since there is friction, the steel ball slides until the angular velocity reaches a value of -110 rad/s, then, starts rolling. One can observe that the angular velocity oscillates over time and reaches it's maximum in absolute value. This maximum value does not decrease during the rolling step and reaches the same maximum value, since there is no more friction and the normal restitution coefficient is equal to 1.

T = 1 s, dt ≈ 1.
T = 5 s, dt = 2.
Figure 4b shows the evolution of the system's mechanical energy over time during the slidingrolling process inside the fixed drum. It decreases brutally during the sliding step because of the friction, then it is conserved during the rolling step. On the other hand, it is notable (see Table 3) that all the methods need a similar total CPU time to compute the solution, except from SAL method which is more expensive in terms of computing time. Besides, the PDAS methods need only 3 NLGS iterations to converge, 2 for IBP method, while SAL method need more than 30 NLGS iterations to converge.

Sedimentation of a collection of balls

We introduce for this third numerical example a representative configuration to show the performances of PDAS methods. The sedimentation of a collection of rigid balls in a box is considered, with particular emphasis on the comparison between the PDAS methods and the standard methods used to compute the local contact impulses. Each rigid ball of the collection has an initial random velocity. Figure 5 illustrates the sedimentation process of 100, 200, 400, 800 and 1600 particles. For each configuration, two different radii are considered for the particles. We provide below the value of the physical settings used for the computation: The simulation is computed over 1s, the numerical parameters related to this experiment are the following: This example is representative of granular flow simulations due to the large number of rigid bodies involved and the large amount of CPU time needed to compute the solution. According to the graphs in Figure 6, one can observe that both EPDAS and IPDAS methods provide the best results in terms of computing time during the whole sedimentation process of 100 and 200 rigid balls in a box. Indeed, the cumulative time required to compute the contact impulses using PDAS method increases slightly over time compared to SAL and IBP methods, which need a larger CPU time to compute the solution. For a large number of particles involved, On the other hand, this numerical experiment assesses not only the performances of PDAS methods, but also the robustness of such methods. Graphs in Figure 7 show the evolution of the cu- mulative NLGS iterations required to the convergence of the contact impulses' computation during the sedimentation process, and it can be noted that both EPDAS and IPDAS curves are quite similar and tend to converge faster than the SAL method, and provide comparable results with IBP method, especially at the end of the sedimentation process. For a large number of particles involved, Table 6 resumes the total NLGS iterations values needed to compute the solution for each numerical method.

ρ = 2600 Kg/m 3 , r = (see table 4),
T =
Remark 4. It should be pointed out that results marked with an asterisk in Table 3 correspond to the simulations conducted with SAL and IBP methods where we noticed some large interpenetrations between rigid particles at the end of the sedimentation process.

Granular materials in a 2D rotating drum

The final numerical example concerns the well-known rotating drum in which we study the behavior of a collection of steel balls in motion. The physical and numerical parameters of the following simulation have been chosen so that it corresponds to the one performed in [START_REF] Maione | DEM investigation of granular flow and binary mixture segregation in a rotating tumbler: Influence of particle shape and internal baffles[END_REF], in order to provide a comparison between the experiment achieved and the numerical results. Figure 8 presents screenshots of the experiment produced with 1.5 Kg of steel balls in a rotating drum at 40 rpm and the DEM simulation related. Steel balls tend to organize themselves in layers: one layer in contact with the wall and the rest within the bed. In this study, PDAS methods are performed to evaluate the average rising height and the mean angle of repose during 1s of simulation. Beyond the numerical comparison, the interest of this example lies in the fact that such experiment assesses the efficiency and the reliability of PDAS methods.

Several simulations were carried out with PDAS methods in order to measure the corresponding mean angle of repose and the rising height. In that respect, we put 100, 200, 400, 800 and then 1600 steel balls in the rotating drum. First, we provide below a description of the physical settings used for this experiment: These simulations are also representative of multi-rigid body contact problems due to the large number of rigid balls considered and the large amount of CPU time needed to compute the solution.

In Figure 9, one can observe that both of repose angle and rising height achieved by PDAS methods are quite similar to the ones observed in Figure 8, and this regardless of the number of particles involved. According to [START_REF] Maione | DEM investigation of granular flow and binary mixture segregation in a rotating tumbler: Influence of particle shape and internal baffles[END_REF], the resulting angle of repose is equal to 26 • . With regard to the convergence of the numerical methods, Graphs in Figure 10 assess once again the efficiency of PDAS methods, that provide very comparable results with IBP method. It is to highlight that both EPDAS and IPDAS curves are quite similar and tend to converge faster than SAL method. For instance, when it comes to deal with 1600 steel balls (see Table 7), the total NLGS iterations needed to compute the whole simulation varies from 8 million to 10 million for PDAS and IBP methods, whereas it takes more than 30 million NLGS iterations to the SAL method to converge. method. Regardless of the number of rigid balls involved in the process, PDAS methods remain the least expensive in terms of computing time. Indeed, as we can see from the graphs in Figure 11, the cumulative CPU time required to compute the contact impulses using PDAS methods increases slightly over time compared to SAL and IBP methods. Furthermore, one can notice according to Table 8 that the CPU time strongly depends on the number of rigid balls involved in the simulation. Indeed, the more rigid balls we have in the rotating drum, the more the gap in CPU time widens, making the PDAS methods the most relevant in this example.

Conclusion

The main contribution of this work consisted of the treatment of frictional contacts within the framework of Active Set method. An adaptation of these methods to solve multi-rigid-body dynamics contact problems has been proposed in the first part of this paper. Indeed, after defining the frictional contact conditions through the NSCD formalism, a general algorithm has been detailed to treat those conditions numerically. A series of academic test cases covering a broad range of granular materials applications has been carried out to assess the relevance of PDAS methods. As a first result, it turns out that their efficiency and speed are better compared to standard and improved methods, the saving in computing time is considerable. Moreover, the numerical experiments reported above demonstrate that computing times widen and are no longer comparable when the number of rigid bodies involved becomes larger, but still faster for PDAS methods compared to other ones. These promising results push us to extend this work to 3D in order to cover a wide range of industrial applications with complex behaviors. Treating these applications in 3D by means of Active Set type method remains, as far as we know, an open work that will be worth discussing in more details. Furthermore, as the number of rigid bodies will increase considerably, we must give serious thought to parallel computing and show the relevance of PDAS methods in terms of computing time. Continuous effort will be devoted to study the parallelization performances of such methods in future works.

By combining ( 3 . 5 )

 35 -(3.22), with D C p n and D C p τ the generalized derivative of C p n and C p τ , respectively, we obtain

  uα,(k) n , uα,(k) τ , p α,(k) n , p α,(k) τ ) + (δ uα,(k+1) n , δ uα,(k+1) τ , δp α,(k+1) n , δp α,(k+1) τ

  stick conditions of the semi-smooth Newton formalism are as follows (in the case of a diagonal matrix W αα τ τ ) uα,(k+1)

  slip 2D conditions of the semi-smooth Newton formalism are as follows uα,(k+1)

9 )-( 4 . 11 )••

 9411 Initialization of the motion: q0 k+1 = qfree k (initialization of the contact impulses with P = 0). -Loop on j ≥ 0 (NLGS), until convergence * Loop on the contacts α: • Computation of the local-global mapping u-= H * (q k+ 1 2 , α) qk ; (4.10) uα,j + = H * (q k+ 1 2 , α) qj k+1 Newton shock law (using formal Moreau velocity) Computation of the frictional contact law:    contact law( ũα,j+1 n , p α,j+1 n ) = .true. f riction law( ũα,j+1 , p α,j+1 ) = .true. (4.14)

  If ( uj+1 , p j+1 ) -( uj , p j ) ≤ , A j+1 n = A jn and A j+1 τ = A j τ stop, else goto (ii).

Figure 1 :

 1 Figure 1: Rolling-Sliding of a steel ball on a plane conveyor belt.

  32 10 -5 s, e n = 1.0, e τ = 1.0, µ = 0.22, γ n = 10, γ τ = 10, r SAL = r IBP = m dt stopping criterion : N LGS = 10 -8

Figure 2 :

 2 Figure 2: Comparison between numerical and analytical solution for the sliding-rolling steel ball. (a) Position; (b) Velocity; (c) Angular velocity.

Figure 3 :

 3 Figure 3: Rolling-Sliding of a steel ball inside a fixed drum.

Figure 4 :

 4 Figure 4: Evolution of two characteristic variables over time for the sliding-rolling steel ball inside a fixed drum. (a) Particle's angular velocity; (b) Particle's mechanical energy.

  Domain size : [0, 0.012] × [0, 0.022] m 2 , g = -9.80665 m/s 2 , v 0 = random

  1 s, dt = 10 -4 s, e n = 1.0, e τ = 1.0, µ = 0.2, γ n = 10 -4 , γ τ = 10 -9 , r SAL = r IBP = m dt stopping criterion : N LGS = 10 -6

Figure 5 :

 5 Figure 5: Screenshots of the sedimentation process of a collection of rigid balls in a box over time. [(a), (b), (c)] 100 rigid balls, [(d), (e), (f)] 400 rigid balls, [(g), (h), (i)] 1600 rigid balls.

Figure 6 :

 6 Figure 6: Evolution of the cumulative CPU time required to compute the contact impulses during the sedimentation process. (a) 100 rigid balls and (b) 200 rigid balls.

Figure 7 :

 7 Figure 7: Evolution of the cumulative NLGS iterations required to the convergence of the contact impulses computation during the sedimentation process. (a) 100 rigid balls and (b) 200 rigid balls.

N

  umber of particles : 100 / 200 / 400 / 800 / 1600, Rotational speed = 40 rpm, ρ = 8000 Kg/m 3 , r ≤ 2.7 10 -3 m, Rotating drum radius : 0.05 m, g = -9.80665 m/s 2 ,

Figure 8 :

 8 Figure 8: Comparison between the experiment and the DEM simulation performed in [46] with steel balls inside a rotating drum.The numerical parameters of the rotating drum experiment are as follows:

Figure 9 :

 9 Figure 9: Perspective views of the rotating drum performed with PDAS methods with (a) 100 steel balls, (b) 400 steel balls, 1600 steel balls

Figure 10 :

 10 Figure 10: Evolution of the cumulative time required to compute the contact impulses during the the rotating drum simulation. (a) 100 steel balls, (b) 400 steel balls and (c) 1600 steel balls.

Figure 11 :

 11 Figure 11: Evolution of the cumulative CPU time required to compute the contact impulses during the rotating drum simulation. (a) 100 steel balls, (b) 400 steel balls and (c) 1600 steel balls.

Table 1 :

 1 Time step Absolute error on position Absolute error on angular Absolute error on velocity 1.32 10 -7 Absolute error between numerical and analytical solution on position, angular velocity and rolling velocity for differents time steps

	6.62 10 -6 1.32 10 -6 2.65 10 -6 1.32 10 -5	7.83 10 -11 2.21 10 -11 1.59 10 -12 1.73 10 -12 5.18 10 -13	1.97 10 -8 4.70 10 -9 1.44 10 -9 1.41 10 -9 2.03 10 -10	5.32 10 -11 1.27 10 -11 3.88 10 -12 3.52 10 -12 3.83 10 -13

Table 2 :

 2

	Numerical method	NLGS it. (Sliding step)	NLGS it. (Rolling step)	CPU time (s)
	EP DAS	2	3	0.664
	IP DAS	2	3	0.692
	SAL	20	20	1.584
	IBP	2	2	0.744

Number of NLGS iterations needed by each method (second column) and total CPU time (third raw) devoted to the computation 6.2. One sliding and rolling steel ball in a fixed drum

Table 3 :

 3 Number of NLGS iterations needed by each method (second column) and total CPU time (third raw) devoted to the computation

Table 4 :

 4 Radii ranging for each sedimentation's configuration

Table 5

 5 resumes the total CPU time values for each numerical method.

			Total CPU time (s)
	Numerical method	100 balls	200 balls	400 balls	800 balls	1600 balls
	EPDAS	2.588	8.412	28.347	107.919	12298.419
	IPDAS	2.562	8.392	28.347	105.408	11515.381
	SAL	22.294	44.214	87.293	182.144 ( * )	18673.187 ( * )
	IBP	13.751	33.533	65.856	186.395 ( * )	16798.548 ( * )

Table 5 :

 5 Total CPU time devoted to the computation of the contact impulses during the sedimentation process for each numerical method and different number of rigid balls.

Table 6 :

 6 Total NLGS iterations needed to the computation of the contact impulses during the sedimentation process for each numerical method and different number of rigid balls.

			Total NLGS iterations
	Numerical method	100 balls	200 balls	400 balls	800 balls	1600 balls
	EPDAS	67726	89347	95995	147445	25578301
	IPDAS	67121	85663	99765	142802	23836195
	SAL	326121	363596	442943	144305 ( * )	2648147 ( * )
	IBP	106194	120975	138970	83404 ( * )	1978658 ( * )

Table 7 :

 7 Total NLGS iterations needed to the computation of the rotating drum simulation for each numerical method and different number of steel balls.

			Total NLGS iterations
	Numerical method	100 balls	200 balls	400 balls	800 balls	1600 balls
	EPDAS	1680882	2020349	4069672	5114280	9422476
	IPDAS	1674177	1933869	4036740	5036011	9999666
	SAL	7836547	8566285	19607909	15600032	32656132
	IBP	1101849	1340917	4121136	3688553	8375264

Table 8 :

 8 Total CPU time devoted to the computation of the rotating drum simulation for each method and different number of steel balls.

			Total CPU time (s)
	Numerical method	100 balls	200 balls	400 balls	800 balls	1600 balls
	EPDAS	69.482	163.918	677.212	1817.413	7492.300
	IPDAS	68.813	153.254	657.011	1804.862	7867.878
	SAL	337.16	755.553	3663.254	6055.488	27425.824
	IBP	124.918	259.370	1941.158	3232.919	16740.516

Declaratio if ioterettt

The authors declare that they have no known competng fnancial interests or personal relatonships ☐ that could have appeared to infuence the work reported in this paper.

The authors declare the following fnancial interests/personal relatonships which may be considered ☒ as potental competng interests:

Barboteu reports a relatonship with European Commission Horizon 2020 Research and Innovaton Programme under the Marie Sklodowska-Curie grant agreement No. 823731 -CONMECH. that includes: funding grants.


Sabri Boubaker
email: sboubaker@em-normandie.fr

Zhenya Liu
email: zhenya.liu@ruc.edu.cn

Shanglin Lu
email: lushanglin@ruc.edu.cn

Yifan Zhang
email: zhangyifan@ruc.edu.cn

Brent Mar

Trading Signal, Functional Data Analysis, and Time Series Momentum $

Keywords: Asset Pricing, Futures Pricing, Time Series Momentum, Trading Signal, Functional Data Analysis G12, G13, G15

Prior empirical results show that the time series momentum portfolio outperformed the buy-andhold benchmark well from 1985 to 2009, but this profitable pattern unexpectedly vanishes after 2009. In this paper, we reconstruct the time series momentum portfolio by applying new trading rules derived from the functional data analysis approaches. Using a dataset that contains 24 commodities from January 2010 to December 2018, our daily-based strategy documents an improvement in the Sharpe ratio of 0.75 compared to 0.07 in terms of the original time series momentum portfolio. This finding offers an alternative strategy for trend-following investors in the commodity futures market.

Introduction

Since the seminal work of [START_REF] Moskowitz | Time series momentum[END_REF], time series momentum (TSM, hereafter) effect has been documented as a universal phenomenon across international financial markets. It was later known as a widely implemented trend-following strategy in the CTA (commodity trading advisor) industry [START_REF] Baltas | Momentum strategies in futures markets and trend-following funds[END_REF][START_REF] Baltas | Demystifying time-series momentum strategies: Volatility estimators, trading rules and pairwise correlations[END_REF]. Based on the dataset of 58 futures contracts from January 1985 to December 2009, [START_REF] Moskowitz | Time series momentum[END_REF] proved that the TSM strategy could generate an impressive alpha, yielding more than 1% monthly excess return. However, adding the data sample from 2010 to 2018, we find that the previous significant performance of the TSM strategy can not be verified after 2009. In this paper, we explore possibilities of improving the performance of TSM strategy through implementing new rules of generating trading signals.

The theory of functional data analysis (FDA, hereafter) provides a data analyzing methodology when observations are curves or are generated from functional processes [START_REF] Ramsay | Some tools for functional data analysis[END_REF]. It has been widely applied in many fields such as financial market modeling. [START_REF] Ramsay | Functional Data Analysis[END_REF] and [START_REF] Horváth | Inference for Functional Data with Applications[END_REF] provide comprehensive details on statistical tools and approaches involved in the theory. In this paper, we firstly apply the FDA approaches to obtain the functional patterns implied in discretely observed asset price sequences. Based on the extracted pattern, we then capture the cyclical and seasonal characteristics of the price trend over a given look-back window. Our proposed trading rules for reconstructing the TSM strategy rely on more information derived by the FDA approaches. We also compare the performance of various strategies under different settings.

The main finding of [START_REF] Moskowitz | Time series momentum[END_REF] is that the sign of cumulative return over the past twelve months predicts future price movement in the next month. However, this predictive power in terms of the TSM portfolio became unclear after 2009. This paper contributes to the literature with modeling price dynamics over TSM look-back window as Hilbert-space-valued random functions rather than scalar random variables, therefore identifying the cyclical and seasonal patterns of price trend involved in the TSM portfolio construction. From an academic perspective, our proposed portfolio constructing method provides insights into understanding the TSM strategy in the CTA industry. From the strategy evaluation perspective, our proposed strategy shows a significant improvement in the Sharpe ratio from 0.07 to 0.75 during the period 2010 to 2018 compared to the original TSM strategy.

The rest of this paper is organized as follows. Section 2 describes our dataset. Section 3 introduces the functional data analysis approach and the trading rules of the improved time series momentum strategy. Section 4 presents the empirical performance of the proposed strategies. Section 5 concludes this paper.

Data

Our dataset is constructed from Bloomberg. It consists of the price series of 24 commodities, which are among the most liquid futures contracts in the world. Following [START_REF] Bessembinder | Systematic risk, hedging pressure, and risk premiums in futures markets[END_REF], [START_REF] Moskowitz | Time series momentum[END_REF], and [START_REF] Asness | Value and momentum everywhere[END_REF], we collect daily settlement prices and adjust for rollovers so that the portfolios are always composed of the most active contracts, which are typically the nearest or next nearest-to-delivery contracts. Likewise, the excess returns over the risk-free rate are computed for all instruments. The sample covers more than fifty years from January 1966 to December 2018, with all 24 commodities became available after 1995. The summary

Methodology

Functional Data Analysis

The FDA theory has been demonstrated to provide sufficient statistical tools for exploring underlying features among discrete observations. Using functional smoothing method, we see a sequence of discrete observations as a curve, i.e., a functional observation. Let p s t,u , tk ≤ u ≤ t denote the prices of security s over latest k observable periods at time t, and the price is considered as an index compiled using the excess returns of corresponding security. Then, discrete prices p s t,u , tk ≤ u ≤ t can be functionalized as a smooth curve in the form of p s t (u), 0 ≤ u ≤ 1. According to the theory, p s t (u) is assumed to be a square-integrable function in L 2 [0, 1], of which the fourth moment is finite. Therefore, it can be approximated by a group of basis functions. In this case, the excess kurtosis and skewness of excess returns presented in Table 1 are satisfied in the functional modeling.

Following previous financial applications of the FDA, we employ the B-Spline basis system to transform the discrete data into smooth functions [START_REF] De Boor | A Practical Guide to Splines[END_REF][START_REF] Cao | A study of data-driven momentum and disposition effects in the Chinese stock market by functional data analysis[END_REF]. The core of functional smoothing is to minimize the sum of squared errors between a linear combination of basis functions and target observations. As a result of dealing with the trade-off between fitting capacity and function smoothness, the roughness penalty is used to generate better results. Moreover, the fourth-order derivative is chosen in the penalty term, considering that we apply the second derivative of the smoothed functions in the latter analysis. 1

Trading Rules

This paper improves the TSM strategy performance in terms of its trading signals generated from the asset past performance over k look-back periods. As proposed in [START_REF] Moskowitz | Time series momentum[END_REF], the long or short trading signal is determined by the sign of cumulative return over past k periods:

sign TSM s,t =        +1, r s t-k,t > 0, -1, r s t-k,t < 0,
where r s t-k,t indicates the cumulative return of security s from tk to t. Consequently, the portfolio return of TSM strategy with 12-months look-back period and 1month holding period across S t securities from time t to t + 1 is given by:

r TSM t,t+1 = 1 S t S t s=1 sign TSM s,t r s t-12,t σ target σ s t r s t,t+1 , (1)
where S t is the number of securities available at time t. We suggest interpreting the economic intuition behind the success of TSM formula, also known as trend-following strategy, through its three decomposed components: i) trading signal, identifying the price trend by comparing the latest price to its counterpart 12 months ago; ii) constant risk parity, scaling the time-varying volatility of each security by a consistent target volatility; 2 iii) equally weighted portfolio, distributing an equal amount of funds to each security within the portfolio following the "naive 1/N rule".3

1 For more details about functional smoothing methodology, we refer to [START_REF] Ramsay | Functional Data Analysis[END_REF] and Horváth & Kokoszka (2012).

2 The annual target volatility is chosen to be 40% due to its similarity to the risk of an average individual stock, as stated in [START_REF] Moskowitz | Time series momentum[END_REF]. We keep the same 40% target volatility to compare our proposed methods with them. Likewise, the ex-ante annualized variance σ 2 t for each security is calculated as follows:

σ 2 t = 261 ∞ i=0 (1 -δ)δ i (r t-1-i -r t ) 2 ,
where δ is chosen to satisfy ∞ i=0 (1δ)δ i = 60 days, and r t is the exponentially weighted average return calculated similarly.

Assuming the asset price follows a continuous-time diffusion process with a stochastic trend, [START_REF] He | Asset allocation with time series momentum and reversal[END_REF] suggest more complicated buy/sell decisions if the goal is to maximize the expected future wealth. Importantly, the FDA method enables us to capture more features implied in price dynamics within the look-back window [tk, t], not merely concentrating on the sign of the partial sum of past returns. Although the partial sum is helpful for filtering noises in asset returns, considering its arithmetic sign as an upward or downward price trend indicator is insufficient. Its deficiency comes in the sense that the aggregated value fails to capture the informative signals of stochastic price trends, particularly in the fragile post-crisis recovery period of the global economy.

We suggest new rules for generating trading signals using additional information derived from the smoothed curve of past price dynamics. On the one hand, the smoothed curve keeps the long-memory pattern of the path-dependent asset price process while filtering out noises. On the other hand, the first-order and second-order derivatives of the smoothed function capture cyclical and seasonal characteristics of the price movements. As an analogy of velocity and acceleration in Physics, the sign consistency of both derivatives identifies a segment of a strong price trend, whereas the inconsistency indicates a segment of a trend reversal. We expect to see improvement in portfolio returns considering that more information is involved in this case when identifying the ongoing price trend, and refer to the new strategy briefly as ITSM.

More specifically, we maintain a strong long (short) trend signal only when the signs of both the first-order and second-order derivatives are positive (negative). The inconsistent cases are suggested to represent a relative weak trend signal. Therefore, we define the trading rules of a conservative strategy as:

sign ITSM c s,t (u) =              +1, Dp s t (u) ≥ 0 and D 2 p s t (u) ≥ 0, 0, Dp s t (u) × D 2 p s t (u) < 0, -1, Dp s t (u) < 0 and D 2 p s t (u) < 0, (2)
and an aggressive strategy as:

sign ITSM a s,t (u) =                    +1, Dp s t (u) ≥ 0 and D 2 p s t (u) ≥ 0, +0.5, Dp s t (u) ≥ 0 and D 2 p s t (u) < 0, -0.5, Dp s t (u) < 0 and D 2 p s t (u) ≥ 0, -1, Dp s t (u) < 0 and D 2 p s t (u) < 0, (3)
where 0 ≤ u ≤ 1, and D(•) and D 2 (•) are the first and second derivative operators. The conservative strategy is designed to clear out position when the trend signal becomes weak, while the aggressive strategy still keeps a half position.

In the empirical studies, we consider testing the performance of ITSM strategies on both monthly and daily rebalancing basis. For the monthly rebalancing strategy, we collect the price index series over the past 12 months and smooth it into function p s t (u), 0 ≤ u ≤ 1 at the end of year t. Then we apply the trading signals generated from Eq. ( 2) and (3) for the following 12 months from t to t + 12. As for the daily rebalancing strategy, we similarly functionalize the daily price indices in each of the past 12 months into 12 smoothed curves at the end of month t. By averaging these 12 curves, we see the trading signals generated from the mean function with Eq.

(2) and (3) as the prediction of price trend pattern over the trading days in the next holding month t to t + 1.

Strategy Performance

Time Series Momentum

Following Eq. ( 1), we calculate the TSM return for all commodities futures contracts with available data from January 1985 to December 2018. Table 2 reports some evaluation statistics of the strategy over the whole sample and two sub-periods. As shown in the last row of Table 2, the well-documented TSM effect in [START_REF] Moskowitz | Time series momentum[END_REF] vanishes after 2009, with a slump in the Sharpe ratio from 0.95 to 0.07.

Imporved Time Series Momentum

As we explained earlier, the FDA theory provides sufficient statistical tools for investigating the past price movements as a curve, not just a scalar variable. We evaluate the strategy performance of applying proposed trading rules instead of the original TSM trading signals in this subsection.

Fig. 1 describes the cumulative returns of the daily and monthly rebalancing conservative ITSM strategies: ITSM d c and ITSM m c , and aggressive ITSM strategies: ITSM d a and ITSM m a from 2010 to 2018. The buy-and-hold (BAH) and TSM strategies are also plotted as comparable counterparts. We further present the strategy performance evaluating statistics, including Sharpe ratio (SR) and maximum drawdown (MDD) for each strategy in Table 3.

As illustrated in Table 3, the TSM strategy significantly outperformed the BAH benchmark from 1985 to 2009. The Sharpe ratio was 0.95 relative to 0.12 in this period. Nevertheless, an almost zero profit TSM investment was demonstrated after 2009 in Fig. 1 suffering long and deep equity drawdowns, even though it was still better than the BAH benchmark. Meanwhile, we see from the subgraphs of Fig. 1 that both the conservative and aggressive ITSM strategies exhibited better profitability relative to the BAH and TSM strategies from 2010 to 2018. Despite no holding positions when the trend signals are weak, it turned out that the conservative ITSM strategies still capture the price trend patterns. The aggressive trading signals worked slightly better than the conservative ones at the daily frequency, achieving a higher terminal wealth in December 2018. Moreover, we report the evaluation statistics according to different commodity classes over two subsamples in Table 3. Generally, all ITSM strategies reported higher SR and smaller MDD from 2010 to 2018, not only in the class of all commodities, but also in four commodity classes of Energy, Agriculture, Livestock, and Metal. Furthermore, the daily rebalanced strategies performed much better than the monthly rebalanced strategies in all commodity classes except for the Metal Sector, since more detailed information of price dynamics is taken into consideration on a daily basis. Additionally, we note that the daily rebalanced conservative strategy ITSM d c and aggressive strategy ITSM d a obtained a Sharpe ratio of 0.75 and 0.73 from 2010 to 2018. Both of them turned out to be a significant improvement relative to 0.07 that of the TSM strategy in the same testing sample.

Performance in Crisis Periods

Since our data sample covers several financial crisis periods, we are also interested in assessing whether the proposed ITSM strategies are resistant to these financial crises. We consider three crisis periods, namely, the Asian Financial Crisis from 1997 to 1999, the Global Financial Crisis from 2007 to 2009, and the European Debt Crisis from 2010 to 2012. Table 4 reports the Sharpe ratio and maximum drawdown statistics for different ITSM strategies in terms of the class of all commodities.

MDD statistics show substantial differences between the proposed ITSM strategies and the TSM strategy in each crisis period. All ITSM strategies present significantly smaller drawdowns compared to the TSM strategy. Except for the Asian Financial Crisis period, the ITSM strategies also document a higher Sharpe ratio than TSM in periods of the Global Financial Crisis and the European Debt Crisis. We note that a smaller drawdown in the crisis period is of higher practical meaningfulness for market participants, particularly fund managers. ‡SR indicates the gross Sharpe ratio of portfolio returns, and MDD represents the maximum drawdown of investment strategy.

Conclusion

The construction of a time series momentum portfolio, usually characterized as the trendfollowing strategy, can be decomposed into three components: trading signal, constant risk parity, and equally weighted portfolio. Based upon this decomposition, we have focused on improving the profitability of the time series momentum portfolio in terms of trading signals over a sample period from January 2010 to December 2018, posterior to that of [START_REF] Moskowitz | Time series momentum[END_REF]. The dataset containing the most liquid contracts of 24 commodities futures from Bloomberg enabled us to investigate the patterns of commodities price movements over a given look-back window.

With the implementation of functional data analysis approaches, we propose new trading rules for identifying the price trend patterns and how the improved time series momentum strategies are constructed. Our trading signals are defined by the signs of the first-order and second-order derivatives of the smoothed functional observations of past price dynamics. As an analogy of velocity and acceleration in Physics, these signals capture the cyclical characteristics of time series momentum effect to some extent. Finally, two differently rebalanced strategies were investigated during the period from 2010 to 2018. All suggested strategies outperform the original TSM strategy in the sense of Sharpe ratio and maximum drawdown of portfolio returns. These findings are consistent for the class of all commodities and four commodities sectors of Energy, Agriculture, Livestock, and Metal. We suggest further research on exploring more financial markets, such as equity index futures, bond futures, currency forwards, and crypto currencies.

Figure 1 :

 1 Figure 1: Equity Curves of the TSM and ITSM investments.

Table 2 :

 2 Performance of the Original Time Series Momentum Strategy.

			1985-2009			2010-2018			1985-2018
		t-stat	SR	D-SR t-stat	SR	D-SR t-stat	SR	D-SR
	Energy
	BRENT	2.46	0.66	0.72	1.21	0.40	0.46	2.65	0.55	0.62
	GASOIL	2.34	0.63	0.73	1.70	0.57	0.55	2.89	0.60	0.65
	WTI	2.16	0.43	0.49	0.24	0.08	0.10	1.99	0.34	0.40
	GASOLINE	2.23	0.47	0.55	0.63	0.21	0.22	2.23	0.40	0.45
	HEATINGOIL	2.39	0.51	0.61	0.75	0.25	0.24	2.44	0.44	0.49
	NATURALGAS 1.47	0.34	0.45	1.20	0.40	0.48	1.87	0.36	0.46
	Agriculture
	COTTON	1.87	0.38	0.41	0.08	0.03	0.03	1.61	0.28	0.32
	COFFEE	0.49	0.10	0.11	0.49	0.16	0.13	0.66	0.11	0.11
	COCOA	0.72	0.15	0.16 -2.56	-0.85	-0.91 -0.78	-0.13	-0.14
	SUGAR	0.17	0.03	0.04	0.21	0.07	0.07	0.26	0.04	0.05
	CORN	1.02	0.20	0.25	0.89	0.29	0.31	1.32	0.23	0.27
	SOYBEAN	0.48	0.10	0.10 -1.64	-0.55	-0.52 -0.54	-0.09	-0.09
	SOYMEAL	0.10	0.02	0.02 -1.74	-0.58	-0.51 -0.93	-0.16	-0.16
	SOYOIL	0.70	0.14	0.14	1.05	0.35	0.40	1.11	0.19	0.20
	WHEAT	1.62	0.32	0.41 -0.57	-0.19	-0.19 1.06	0.18	0.21
	Livestock
	LIVECATTLE	2.33	0.47	0.53 -0.90	-0.30	-0.33 1.60	0.27	0.31
	LEANHOGS	1.59	0.33	0.47 -0.10	-0.03	-0.04 1.28	0.23	0.31
	Metal
	ALUMINUM	0.86	0.19	0.18	0.93	0.31	0.38	1.22	0.22	0.23
	COPPER	2.55	0.54	0.69	0.37	0.12	0.12	2.41	0.43	0.51
	NICKEL	2.38	0.51	0.59 -1.27	-0.42	-0.48 1.40	0.25	0.29
	ZINC	1.93	0.43	0.46 -0.73	-0.24	-0.27 1.26	0.24	0.26
	GOLD	2.01	0.40	0.40	0.97	0.32	0.40	2.23	0.38	0.40
	SILVER	0.98	0.20	0.20	0.50	0.17	0.18	1.09	0.19	0.19
	PLATINUM	0.92	0.18	0.17 -0.07	-0.02	-0.03 0.77	0.13	0.13
	ALL	4.73	0.95	1.08	0.21	0.07	0.06	4.02	0.69	0.72
	Notes: †The time series momentum strategy returns are computed following Eq. (1). For comparison with Moskowitz et al. (2012), we divide the whole sample into two subsamples: 1985 to 2009 and 2010 to 2018.
	‡SR indicates the gross Sharpe ratio of portfolio returns, and D-SR indicates the downside-risk Sharpe ratio introduced by Ziemba (2005).

Table 3 :

 3 Performance Evaluating Statistics of the Improved Time Series Momentum Strategies.

				1985-2009			2010-2018
			SR MDD(%) t-stat	SR MDD(%) t-stat
		BAH	0.12	53.25	0.61	-0.38	54.71	-1.14
		TSM	0.95	24.08	4.73	0.07	31.07	0.21
	All Commodities	ITSM d c ITSM m c ITSM d a ITSM m a	0.42 0.89 0.46 0.63	21.02 10.35 21.17 15.39	2.12 4.14 2.32 3.17	0.73 0.33 0.75 0.20	12.32 15.57 12.54 17.01	2.19 0.89 2.25 0.55
		BAH	0.24	76.19	1.20	-0.35	74.33	-1.06
		TSM	0.55	66.37	2.73	0.49	38.29	1.46
	Energy Sector	ITSM d c ITSM m c ITSM d a ITSM m a	-0.01 0.27 -0.06 0.23	70.93 67.30 77.81 68.86	-0.05 1.34 -0.32 1.15	0.92 0.62 0.98 0.50	25.75 17.31 29.37 26.40	2.74 1.85 2.92 1.50
		BAH	-0.15	73.12	-0.76	-0.29	60.24	-0.86
		TSM	0.36	57.31	1.78	-0.34	58.30	-1.02
	Agriculture Sector	ITSM d c ITSM m c ITSM d a ITSM m a	0.42 0.65 0.42 0.54	41.88 21.09 39.34 20.13	2.11 3.26 2.10 2.69	0.46 -0.26 0.45 -0.28	27.16 34.82 27.91 39.95	1.36 -0.78 1.36 -0.83
		BAH	-0.12	83.03	-0.59	-0.24	50.68	-0.71
		TSM	0.55	76.54	2.72	-0.22	64.58	-0.65
	Livestock Sector	ITSM d c ITSM m c ITSM d a ITSM m a	0.39 0.40 0.29 0.22	48.27 41.46 55.24 51.56	1.93 1.96 1.38 1.08	-0.06 0.17 0.53 -0.04	68.39 16.68 46.96 26.47	-0.17 0.51 1.58 -0.12
		BAH	0.16	-64.43	0.82	-0.18	54.43	-0.53
		TSM	0.60	44.71	3.01	0.07	47.60	0.20
	Metal Sector	ITSM d c ITSM m c ITSM d a ITSM m a	0.13 0.31 0.10 0.25	45.75 18.61 52.15 21.59	0.65 1.56 0.51 1.24	0.01 0.35 0.13 0.25	28.62 15.32 17.00 22.02	0.02 1.05 0.38 0.74
	Notes: †BAH is short for the buy-and-hold portfolio, which is equal-weighted and long-only. TSM is the 12-1 time series momentum strategy in line with Eq. (1). ITSM d c , ITSM m c , ITSM d a , and ITSM m a are the conservative and aggressive improved time series momentum strategies with daily and monthly rebalancing.
	‡For the convenience of comparing with Moskowitz et al. (2012), we present the strategy performance on two divided subsamples: 1985 to 2009 and 2010 to 2018.
	§SR indicates the gross Sharpe ratio of portfolio returns, and MDD represents the maximum drawdown of investment strategy.

Table 4 :

 4 Performance of ITSM Strategies during the Crisis Periods.

		Asian Financial Crisis	Global Financial Crisis	European Debt Crisis
			(1997-1999)			(2007-2009)			(2010-2012)
	Strategy	SR MDD(%) t-stat	SR MDD(%) t-stat	SR MDD(%) t-stat
	BAH	-0.79	40.13	-1.35	-0.13	53.25	-0.23	0.11	26.28	0.20
	TSM	1.03	13.48	1.75	0.18	24.08	0.31	-0.40	31.07	-0.68
	ITSM d c ITSM m c ITSM d a ITSM m a	0.47 0.98 0.43 0.57	7.97 4.73 9.21 6.76	0.81 1.68 0.74 0.98	0.59 0.31 0.09 0.22	6.98 9.95 13.07 15.39	1.00 0.53 0.15 0.38	1.16 0.02 0.88 0.10	11.31 15.57 12.54 17.01	1.98 0.03 1.49 0.17
	Notes: †BAH is short for the buy-and-hold portfolio, which is equal-weighted and long-only. TSM is the 12-1 time series momentum strategy in line with Eq. (1). ITSM d c , ITSM m c , ITSM d a , and ITSM m a are the conservative and aggressive improved time series momentum strategies with daily and monthly rebalancing.

Note that we only focus on exploring the possibilities of changing the trading signal component in this study.

Appendix A. Commodity Futures Dataset

We provide additional details about the commodities futures contracts involved in this paper. Our dataset from Bloomberg consists of 24 commodities: six of them are in the energy sector (Brent Crude, Gas Oil, WTI Crude, RBOB Gasoline, Heating Oil, Natural Gas), nine of them are in the agriculture sector (Cotton, Coffee, Cocoa, Sugar, Corn, Soybeans, Soy Meal, Soy Oil, Wheat), two of them are in the livestock sector (Live Cattle, Lean Hogs), and seven of them are in the metal sector (Aluminum, Copper, Nickel, Zinc, Gold, Silver, Platinum). For all commodities futures contracts, the exchanges at which they are traded, and tickers from Bloomberg are given in the following table. Note: The first and second nearest future contract can be obtained from Bloomberg by substituting 1 and 2 with 'x' in the contract ticker. For example, 'CO1 Comdty' and 'CO2 Comdty' are the tickers of the first and second nearest active contract of brent crude oil.

Declaration of Competing Interest

None.

Highlights:

1. Empirical results present that the well-documented time series momentum effect on the commodity futures market vanishes after 2009. 2. The past asset price dynamics are employed to generate new trading signals, other than the sign of past asset returns. 3. The discrete sample observations of the asset prices are functionalized to curves using the functional smoothing method. 4. The first-order and second-order derivatives of the functionalized curves are designed to generate new trading signals for the return formula of time series momentum. 5. Compared with the original time series momentum, the empirical application of the newly generated trading signals shows a significant improvement in the Sharpe ratio from 0.07 to 0.75 after 2009.


Marine Poulain
email: marine.poulain@hopital-foch.com

Camille Fossard

Achraf Benammar

Eric Farfour

Stephanie Ranga

Meryem Filali

Jessica Vandame

Paul Pirtea

Marc Vasse

Catherine Racowsky

Jean-Marc Ayoubi

Mcu-Ph

Validation of the French IVF guidelines during the COVID-19 pandemic

Keywords: Covid-19, French guidelines, IVF, follicular fluid

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Infertility is a disease [START_REF] Zegers-Hochschild | on behalf of ICMART and WHO. The International Committee for Monitoring Assisted Reproductive Technology (ICMART) and the World Health Organization (WHO) revised glossary on ART terminology[END_REF], it is not elective, and requires treatment in a timely manner. In vitro fertilization (IVF) centers must therefore balance pressures to treat against potential risks of introducing SARS-CoV-2 virus into their practices.

All three major societies of reproductive medicine, ASRM, ESHRE and IFFS have independently published recommendations for the management of patients planning to undergo IVF treatment during the pandemic. More recently, these societies published a joint statement reaffirming the importance for continued reproductive care during this unprecedented time [START_REF] Veiga | Assisted reproduction and COVID-19: A joint statement of ASRM[END_REF]. Although specific to administration of reproductive care, each society emphasizes the importance of strict adherence to official government guidelines, given differences in national, regional, and local regulations and variations in viral penetration, transmission rates and changing conditions. Shortly after WHO declared the Covid-19 outbreak a pandemic, the French authorities issued guidelines for the safe administration of IVF treatment. Several updates have been released. All recommend use of screening symptom questionnaires, and while none mandate testing for SARS-CoV-2, they advise testing in symptomatic patients.

Little is known regarding whether use of a symptom questionnaire without testing is reliable for triage of asymptomatic patients. This question is particularly important for IVF laboratories, that must maintain the safety not only of laboratory personnel and patients, but also the gametes and embryos in their care. Therefore, the present study was undertaken to assess efficacy of using a symptom questionnaire for screening IVF patients for COVID-19 infection. Secondarily, we investigated the prevalence of SARS-CoV-2 RNA detection in the follicular fluid (FF).

METHODS

Patients planning IVF from June 2020 through February 2021 were stratified into two groups: Group 1 with planned IVF, June through September 2020; and Group 2 with planned IVF, October 2020 through February 2021). Both groups completed the symptom questionnaire for Covid-19 according to the mandatory prevailing French guidelines. However, Group 2 also underwent RT-PCR testing for SARS-CoV-2, 72h before the retrieval, which was not specifically recommended by the guidelines.

Questionnaires were administered two weeks before stimulation, at each appointment, on the day before retrieval and two weeks after using the Wistim application, as recommended by the French biomedicine agency (French Biomedicine Agency, May 2020). Group 1 patients consented for: 1) FF testing for SARS-CoV-2, which was performed after RNA extraction from 1 ml samples (Nimbus, Seegene*) and analysis using a SARS-CoV-2 multiplex RT-PCR using CFX96 (Biorad*); and 2) participation in a comprehensive epidemiological evaluation by phone interview more than 6 months post-cycle to determine Covid-19 exposure during 6 months before and after retrieval.

Cycle cancellation rates between the two groups were compared and incidence of Covid-19 infection among IVF personnel recorded. The percentage of FF samples positive for SARS-CoV-2 was determined, as was the proportion of Group 1 patients interviewed post-IVF who had COVID-19 symptoms (anosmia and ageusia) or a positive RT-PCR test after their IVF cycle.

That none of the 183 FF samples from Group 1 were positive for SARS-CoV-2 gave us reassurance that despite not having a PCR test before retrieval, the safety of our patients' oocytes was ensured.

Recently, [START_REF] Rajput | Absence of SARS-CoV-2 (COVID-19) virus within an active IVF laboratory using strict patient screening and safety criteria[END_REF] also reported absence of SARS-CoV-2 viral RNA in FF, although their patients had tested negative for Covid-19, 3-4 days prior to the procedure. While neither observation negates FF being a viral reservoir, a single case report failing to detect viral RNA in FF of a SARS-CoV-2 positive woman suggests that viral transmission in FF may not occur [START_REF] Demirel | Failure to detect viral RNA in follicular fluid aspirates from a SARS-CoV-2-positive woman[END_REF].

In summary, our results provide reassurance that a symptom questionnaire completed by patients at the intervals recommended by the French guidelines, is adequate for effectively triaging IVF patients.

Figure 1 :

 1 Figure 1: Design and results of study assessing efficacy of using a symptom questionnaire for triaging IVF patients

Screened with questionnaire and RT-PCR test

				Oct 2020 -Feb 2021
				N=577
	Responded	Consented for FF	Cycles cancelled	Cycles cancelled	Cycles cancelled
	to the	testing for	from questionnaire	from	from RT-PCR
	interview	SARS-CoV-2 RNA	responses	questionnaire	Test
	N=171	N=183	N=0	responses	results
	(93.4%)	(85.5%)	(0%)	N=1/577 (0.2%)	N=7/577 (1.2%)
	Patients with Covid-19 symptoms or positive RT-PCR test	FF samples positive for	Total with	Total with
	within 6 months pre-	SARS-CoV-2	cycle	cycle
	or post IVF N=16/171 (9.3%)	RNA N=0/183 (0%)	cancellation N=0/214 (0%)	cancellation N=8/577 (1.4%)

RESULTS

The results are shown in Figure 1. Of the 791 patients who planned IVF during the study period, 214 comprised Group 1 and 577 Group 2. Cycle cancellation rates based on questionnaire responses were 0% (0/214) and 0.2% (1/577) in Groups 1 and 2, respectively (p=1.000; Fisher's Exact Test); the one patient in Group 2 had a fever and confirming PCR positive test. Seven other patients in Group 2 were cancelled: all 7 were asymptomatic for COVID-19 but tested positive, giving an overall cancellation rate in Group 2 of 1.4% (8/577).

183 patients in Group 1 consented to donate samples for FF testing (183/214; 85.5%), none of which tested positive for SARS-CoV-2 (0/183; 0%). 171 patients were interviewed post-cycle (171/214; 93.4%), sixteen of whom (16/171; 9.3%), declared they had Covid-19 symptoms or a positive RT-PCR test, but none within 2 months pre-or post IVF.

DISCUSSION

Upon resumption of activities after the first lock-down early in 2020, French IVF centers set up a triage system using a symptom questionnaire based on the recommendations of the French biomedicine agency, ESHRE and the prevailing literature. Subsequently, some centers, including ours, have set up a systematic RT-PCR test for SARS-CoV-2, performed 72 hours preceding the retrieval. However, accuracy of this test is limited by viral load, specimen handling and the patient's window of positivity, which together may lead to unreliable results of patient status [START_REF] Demirel | Failure to detect viral RNA in follicular fluid aspirates from a SARS-CoV-2-positive woman[END_REF]. The test is also expensive, not all IVF centers have access to testing, and during a peak epidemic, testing can add complexity to the routine of a laboratory. We therefore considered it important to study the relevance of such testing in the setting of clinical IVF.

In the present study, the cycle cancellation rate of patients who complied with the guidelines by answering a symptom questionnaire was insignificantly different from those who also had a RT-PCR test for SARS-CoV-2 (0% versus 1.4%; p=0.116). Moreover, no IVF personnel developed Covid-19 during the study period. We therefore conclude that use of a symptom questionnaire is effective in screening patients planning to undergo IVF and that adding a PCR test does not improve patient triage. Finally, none of the patients we interviewed post-IVF in Group 1 developed pathognomonic symptoms of Covid-19 or had a positive PCR test in the two months before or after their retrieval, which is once again reassuring about our practices without testing.