G Deaconu

C Louembet

A Théron

Constrained periodic spacecraft relative motion using non-negative polynomials

Keywords: impulsive orbital rendezvous, periodic relative motion, non-negative polynomials, semi-definite programming

A new method for obtaining constrained periodic relative motion between spacecraft on Keplerian orbits is presented. The periodic relative trajectory is required to evolve autonomously inside a tolerance box centered in a specified position. Unlike the classical time-sampling approaches, our method guarantees continuous satisfaction of the constraints on infinite horizon. This is done by reformulating the tolerance box constraints on the relative trajectory as conditions of non-negativity of some polynomials. The resulting problem is solved using semi-definite programming.

I. INTRODUCTION

In the recent years, new rendezvous missions and proximity operations between spacecraft have been considered, including sample returns, repairing, refueling, upgrading, and other on-orbit servicing [START_REF] Irvin | Fuel-Optimal Maneuvers for Constrained Relative Satellite Orbits[END_REF]. The ability to design fuel efficient maneuvers in order to reach naturally periodic relative orbits, confined in a specified region in space, will be an important step towards the success of these missions. We propose a new method for obtaining an open loop maneuvers plan leading to this kind of bounded autonomous periodic motion.

The relative motion between spacecraft on Keplerian orbits is naturally bounded [START_REF] Gurfil | Relative Motion between Elliptic Orbits: Generalized Boundedness Conditions and Optimal Formationkeeping[END_REF]. Moreover, the relative motion can be initialized to result as periodic when viewed from a local reference frame. In the unperturbed case, the necessary and sufficient condition for periodicity is that the semi-major axis of the spacecraft orbits are equal [START_REF] Alfriend | Dynamics and Control of Spacecraft Formations: Challenges and Some Solutions[END_REF]. The period of the resulting relative trajectory is then equal to the orbital period of the spacecraft.

Using Cartesian coordinates and a local frame attached to one of the spacecraft, Inalhan obtained in [START_REF] Inalhan | Relative Dynamics and Control of Spacecraft Formations in Eccentric Orbits[END_REF] a periodicity condition at perigee, for the linearized relative motion and for arbitrary eccentricity. This condition was then used for spacecraft formation initialization. Inalhan's condition was later explicitly generalized for any true anomaly in [START_REF] Sengupta | Relative Motion and the Geometry of Formations in Keplerian Elliptic Orbits[END_REF], where the effect of the eccentricity on the shape of the periodic relative trajectory was analyzed. Gurfil wrote the energy-matching condition for periodicity in [START_REF] Gurfil | Relative Motion between Elliptic Orbits: Generalized Boundedness Conditions and Optimal Formationkeeping[END_REF], obtaining a sixth degree polynomial equation, valid for the non linear dynamics.

The previous conditions guarantee that the resulting relative trajectory will be periodic, without giving any insight on its shape or its size. In [START_REF] Gurfil | Manifolds and Metrics in the Relative Spacecraft Motion Problem[END_REF] and [START_REF] Kholshevnikov | Natural Metrics in the Spaces of Elliptic Orbits[END_REF], Gurfil and Kolshevnikov analytically studied the minimum and maximum distances between two spacecraft on elliptical orbits. Their method amounts to solving an eighth degree trigonometric polynomial and obtaining the true anomalies corresponding to the extremal distances, in a worst case analysis. However, the dimensions of the relative trajectory are highly dependent on the initial conditions of the motion, which are not considered in the aforementioned studies.

In [START_REF] Breger | Cooperative Spacecraft Formation Flying: Model Predictive Control with Openand Closed-Loop Robustness[END_REF] model predictive control (MPC) and linear programming (LP) are used to design a fuel optimal maneuvers plan, leading to autonomous periodic relative motion confined inside a specified tolerance region. The method is based on the propagation of the autonomous relative motion starting at the end of the plan for one additional orbital period. Periodicity is ensured by imposing the initial and the final state of the autonomous motion to match. To control the dimensions of the periodic trajectory, its inclusion in the tolerance box is explicitly checked in a finite number of time samples along the propagation. Because of the nature of the orbital dynamics, an autonomous relative trajectory that remains inside the box for one orbit and begins and terminates at the same state is a nominal invariant set [START_REF] Breger | Cooperative Spacecraft Formation Flying: Model Predictive Control with Openand Closed-Loop Robustness[END_REF]. Thus attaining this kind of trajectory prevents any future constraints violation.

The LP formulation has the advantage of providing a general framework for minimizing fuel consumption while including various types of state and actuator constraints [START_REF] Tillerson | Co-ordination and control of distributed spacecraft systems using convex optimization techniques[END_REF]. The drawback is that it cannot guarantee continuous satisfaction of the constraints since the behavior between the time samples is not controlled.

A new method for designing a fuel optimal maneuvers plans leading to this type of bounded periodic relative motion is presented in this paper. This approach no longer requires the explicit propagation of the autonomous relative trajectory at specified time samples. It is based on the fact that the periodic motion can be represented by rational expressions. We show that, by using these rational expressions, the dimension constraints for the periodic trajectory can be written in terms of non-negativity of some polynomials. The main advantage of this method is that it guarantees constraints satisfaction continuously in time.

The paper is organized as follows. Section II first gives a brief description of the spacecraft relative dynamics and their closed-form solutions. Then the linear condition used to ensure the periodicity of the relative trajectory is presented.

Finally, the rational expressions for the propagation of the autonomous periodic relative motion are derived. Our main contribution is detailed in Section III, where the previously obtained rational expressions are used to write the problem of designing a periodic relative trajectory contained in a tolerance box as a polynomial optimization problem. This polynomial optimization problem is transformed into a convex semi-definite program (SDP), for which efficient solvers are available. An always feasible formulation of the problem is also presented, by considering the dimensions of the tolerance box as part of the optimization variables. The methodology is illustrated through some impulsive rendezvous scenarios in Section IV.

II. RELATIVE DYNAMICS AND PERIODICITY

A. Relative Motion Modeling

Let us consider the relative motion between two spacecraft on elliptic Keplerian orbits. The reference orbit corresponds to the leader spacecraft M 1 while the other orbit corresponds to the follower spacecraft M 2 .

Consider an Earth-centered inertial (ECI) frame with the orthonormal basis B 0 = ( X0 , Ȳ0 , Z0 ). The relative motion is described in the rotating Cartesian local-vertical/localhorizontal (LVLH) frame attached to the leader. The corresponding basis is B 1 = ( x 1 , y 1 , z 1 ), with z 1 lying along the radius vector from the satellite to the center of the Earth, y 1 normal to the plane defined by the position and velocity vectors of M 1 in the ECI frame, in opposite direction to the angular momentum, and x 1 completing the basis (Fig. 1).

Let r = x y z T be the position of the follower spacecraft in the LVLH frame of the leader. The linearized relative dynamics with respect to the orbit of the leader are described by the well-known Tschauner-Hempel equations [START_REF] Tschauner | Elliptic orbit rendezvous[END_REF]:

ẍ = 2 ν ż + νz + ν2 x -µ R 3 x + u x ÿ = -µ R 3 y + u y z = -2 ν ẋ -νx + ν2 z + 2 µ R 3 z + u z (1)
where ν and R are the true anomaly and the norm of the radius of the leader, µ is the gravitational parameter of the

x 1 z 1 ν a O X 0 Y 0 Z 0 M 1 M 2 r
Fig. 1. The local frame and the relative position of the spacecraft Earth:

ν = µ a 3 (1-e 2 ) 3 (1 + e cosν) 2 , R = a(1 -e 2 ) 1 + e cosν , (2)
and a and e are the semi-major axis and the eccentricity of the leader. After replacing time as the independent variable with the true anomaly:

d(•) dt = d(•) dν dν dt = (•) ′ ν, d 2 (•) dt 2 = ν2 (•) ′′ + ν(•) ′ (3)
and scaling the variables by:   x ỹ z

 = (1 + e cosν)   x y z     x′ ỹ′ z′   = (1 + e cosν)   x ′ y ′ z ′   -e sin ν   x y z   , (4)
equations ( 1) become:

x′′ = 2z ′ + ũx ỹ′′ = -ỹ + ũy z′′ = 3 1+e cos ν z -2 x′ + ũz (5)
In the sequel, the ˜sign is used to mark the variables after the variables change ( 4). If we define the state vector

X(ν) = x(ν) ỹ(ν) z(ν) x′ (ν) ỹ′ (ν) z′ (ν) T and the input vector ũ(ν) = ũx (ν) ũy (ν) ũz (ν) T , system (5)
can be written as a periodic state-space model:

d X(ν) dν = Ã(ν) X(ν) + B ũ(ν) (6)
with:

Ã(ν) =         0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 2 0 -1 0 0 0 0 0 0 3 1+e cos ν -2 0 0         , B =         0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1         (7)
Starting from system (6), Yamanaka and Ankersen give in [START_REF] Yamanaka | New state transition matrix for relative motion on an arbitrary elliptical orbit[END_REF] a transition matrix Φ which enables the propagation of the relative motion starting from an initial condition X(ν 0 ), under impulsive control:

X(ν) = Φ ν ν 0 X(ν 0 ) + ∑ i Φ ν ν i B ũi . (8)

B. Periodic motion propagation

An arbitrary state X(ν 0 ) is an initial state for a naturally periodic motion if it is equal to the autonomously propagated state after one period. This condition can easily be expressed using the Yamanaka-Ankersen transition matrix:

X(ν 0 + 2π) = Φ ν 0 +2π ν 0 X(ν 0 ) = X(ν 0 ) (9)
After some calculations, the following periodicity condition is obtained [START_REF] Sengupta | Relative Motion and the Geometry of Formations in Keplerian Elliptic Orbits[END_REF]:

x′ (ν 0 ) = (2 + 3e cosν 0 + e 2 ) (1 + e cosν 0 ) 2 z(ν 0 ) + e sin ν 0 (1 + e cosν 0 ) z′ (ν 0 ). (10)
Considering u i = 0 in [START_REF] Breger | Cooperative Spacecraft Formation Flying: Model Predictive Control with Openand Closed-Loop Robustness[END_REF] and assuming that the initial state X(ν 0 ) satisfies [START_REF] Tschauner | Elliptic orbit rendezvous[END_REF], the equations for the propagation of the periodic motion can be deduced:

x(ν) = (2 + e cosν)(d 1 sin ν -d 2 cos ν) + d 3 ỹ(ν) = d 4 cos ν + d 5 sin ν z(ν) = (1 + e cosν)(d 1 cos ν + d 2 sin ν)
, ν ≥ ν 0 [START_REF] Yamanaka | New state transition matrix for relative motion on an arbitrary elliptical orbit[END_REF] where d i , i = 1...5 are constants linearly depending on the initial state X(ν 0 ).

Let D = d 1 d 2 d 3 d 4 d 5
T , then:

D = C(ν 0 ) X(ν 0 ) (12)
where the matrix C(ν 0 ) is defined by:

C(ν 0 ) =         0 0 2ec 2 0 +c 0 -e (1+ec 0 ) 2 0 0 -s 0 1+ec 0 0 0 s 0 (1+2ec 0 ) (1+ec 0 ) 2 0 0 c 0 1+ec 0 1 0 es 0 (2+ec 0 ) (1+ec 0 ) 2 0 0 2+ec 0 1+ec 0 0 c 0 0 0 -s 0 0 0 s 0 0 0 c 0 0         (13)
with c 0 = cos ν 0 and s 0 = sin ν 0 . By replacing the true anomaly ν with w ∈ R such that:

w = tan( ν 2 ), cos ν = 1 -w 2 1 + w 2 , sin ν = 2w 1 + w 2 (14)
a rational form for equations ( 11) is obtained:

x(w) = p x 1 w 4 + p x 2 w 3 + p x 3 w 2 + p x 4 w + p x 5 (1 + w 2 ) 2 ỹ(w) = p y 1 w 2 + p y 2 w + p y 3 1 + w 2 z(w) = p z 1 w 4 + p z 2 w 3 + p z 3 w 2 + p z 4 w + p z 5 (1 + w 2 ) 2 (15)
Denoting P x = [p x 1 p x 2 p x 3 p x 4 p x 5 ] T , P y = [p y 1 p y 2 p y 3 ] T , P z = [p z 1 p z 2 p z 3 p z 4 p z 5 ] T , it comes:

P x = C x D P y = C y D P z = C z D ( 16
)
where:

C x =      0 2 -e 1 0 0 4 -2e 0 0 0 0 0 2e 2 0 0 4 + 2e 0 0 0 0 0 -2 -e 1 0 0      C y =   0 0 0 -1 0 0 0 0 0 2 0 0 0 1 0   C z =      e -1 0 0 0 0 0 2 -2e 0 0 0 -2e 0 0 0 0 0 2 + 2e 0 0 0 e + 1 0 0 0 0      (17)
A similar approach is used in [START_REF] Jiang | Study on Relative Orbit Geometry of Spacecraft Formations in Elliptical Reference Orbits[END_REF] in order to obtain algebraic expressions for the propagation of the periodic motion, in the case of a periodic motion starting at perigee. The objective of [START_REF] Jiang | Study on Relative Orbit Geometry of Spacecraft Formations in Elliptical Reference Orbits[END_REF] was to determine the type of quadratic surface on which the relative orbit laid. As far as we are concerned, the rational expressions (15) will be used in the purpose of designing maneuvers plans leading to bounded periodic motion after control.

III. TRAJECTORY OPTIMIZATION FOR AUTONOMOUS

RENDEZVOUS

The fixed time impulsive rendezvous problem consists in determining fuel optimal impulsive maneuvers that steer the follower spacecraft from an initial state X(ν 1 ) to a desired state

X f = x f y f z f ẋ f ẏ f ż f
T . An impulsive maneuver involves an ideally instantaneous change ∆V in the spacecraft velocity. The number N of impulsive thrusts is known a priori and so are the thrusting instants ν 1 , ..., ν N .

Hence the optimization variables are ∆V 1 , ..., ∆V N ∈ R 3 , the coordinates of the thrusts in the local basis of the leader.

A. Proximity periodic relative motion

For our purpose, the terminal condition is relaxed by no longer requiring the follower spacecraft to reach an exact state X(ν N ) = X f . The objective is to minimize the fuel consumption necessary to reach a final state that can guarantee periodic movement after control, on a trajectory contained in a tolerance region X tol , while taking into account actuators saturation. The optimization problem can be written as: min

∆ Ṽ ∑ N i=1 ∆ Ṽi 1 s.t.      ∆ Ṽi ≤ ∆ Ṽi , ∀ i = 1...N X(ν N ) satisfies (10) ( x(ν), ỹ(ν), z(ν)) ∈ Xtol , ∀ν ≥ ν N ( 18
)
where

∆ Ṽ = ∆ Ṽ1 • • • ∆ ṼN T ∈ R 3N
and ∆ Ṽi ∈ R 3 are the thrusters saturation levels. X(ν N ) is the state at the end of control, which considering (8) can be expressed as:

X(ν N ) = A N + B N ∆ Ṽ (19)
with:

A N = Φ ν N ν 1 X(ν 1 ) B N = Φ ν N ν 1 B • • • Φ ν N ν N-1 B B . (20)
Since the 1-norm criterion is only piecewise linear, slack variables Z ∈ R 3N are introduced in order to write J in a linear form. Z are such that:

-Z i ≤ ∆ Ṽi ≤ Z i Z i ≤ ∆ Ṽi , ∀ i = 1, • • • , N (21)
Let M(ν N ) be:

M(ν N ) = 0 0 -2+3e cosν N +e 2 (1+e cosν N ) 2 1 0 -e sin ν N 1+e cosν N (22)
Then the periodicity constraint (10) can be written as:

M(ν N ) X(ν N ) = 0 (23)
The tolerance region X tol is a box centered on the final position X f , defined by X tol = x tol y tol z tol . The variable change (4) must be taken into consideration when writing the tolerance box constraints on the autonomously propagated trajectory:

(1 + e cosν)(x f -x tol ) ≤ x(ν) ≤ (1 + e cosν)(x f -x tol ) (1 + e cosν)(y f -y tol ) ≤ ỹ(ν) ≤ (1 + e cosν)(y f -y tol ) (1 + e cosν)(z f -z tol ) ≤ z(ν) ≤ (1 + e cosν)(z f -z tol ) (24)
Constraints (24) must be respected for all instants starting from the end of the plan (∀ν ≥ ν N ). The optimization problem (18) becomes: min

Z,∆ Ṽ N ∑ i=1 Z i s.t.                    -Z i ≤ ∆ Ṽi ≤ Z i Z i ≤ ∆ Ṽi , ∀ i = 1...N M(ν N ) X(ν N ) = 0 x f -xtol ≤ x(ν) ≤ x f + xtol ỹ f -ỹtol ≤ ỹ(ν) ≤ ỹ f + ỹtol z f -ztol ≤ z(ν) ≤ z f + ztol , ∀ν ≥ ν N (25)
The LP approach for solving (25) presented in [START_REF] Breger | Cooperative Spacecraft Formation Flying: Model Predictive Control with Openand Closed-Loop Robustness[END_REF] requires the discretization of constraints (24) over a specified interval [ν N+1 , ν N+q ] after control:

x f -xtol ≤ x(ν k ) ≤ x f + xtol ỹ f -ỹtol ≤ ỹ(ν k ) ≤ ỹ f + ỹtol z f -ztol ≤ z(ν k ) ≤ z f + ztol , ∀ k ∈ {N + 1, ..., N + q} (26)
with q the number of discrete points where constraints are explicitly checked. This approach cannot guarantee that the box constraints (24) are not violated between the discretization points.

The variable change ( 14) can be used to define rational bounds for each coordinate. Please note that the methodology is detailed only for the x axis. The y and z axis follow the same procedure.

xmax (w) = (1-e)w 2 +1+e 1+w 2 (x f + x tol ) xmin (w) = (1-e)w 2 +1+e 1+w 2 (x f -x tol ) (27)
Using equations ( 15) and ( 27), the tolerance box constraints (24) can be written as rational inequalities:

xmin (w) ≤ x(w) ≤ xmax (w), ∀w ∈ R (28)
Let x and x be:

x(w) = -x(w) + xmax (w) = 1 (1+w 2 ) 2 Γ x (w) x(w) = x(w) -xmin (w) = 1 (1+w 2 ) 2 Γ x (w) (29)
where Γ x and Γ x are the following polynomials:

Γ x (w) = 5 ∑ i=1 γ x i w i-1 Γ x (w) = 5 ∑ i=1 γ x i w i-1 (30) whose coefficients γ x = γ x 1 γ x 2 γ x 3 γ x 4 γ x 5 T and γ x = γ x 1 γ x 2 γ x 3 γ x 4 γ x 5
T are given by:

γ x = -P x + T (x f + x tol ) γ x = P x -T (x f -x tol ) (31)
and T = 1e 0 2 0 1 + e T . Using ( 12), ( 16) and ( 19) it comes:

γ x = -C x C(ν N )B N ∆ Ṽ -C x C(ν N )A N + T (x f + x tol ) γ x = C x C(ν N )B N ∆ Ṽ + C x C(ν N )A N -T (x f -x tol ) (32)
Hence the coefficients γ x and γ x of the polynomials Γ x and Γ x depend linearly on the optimization variables ∆ Ṽ .

The box constraints (28) can be expressed as nonnegativity conditions on x and x. However, since the denominator (1 + w 2 ) 2 is always strictly positive, these conditions are equivalent to the non-negativity of the polynomials Γ x and Γ x :

Γ x (w) ≥ 0 Γ x (w) ≥ 0 , ∀w ∈ R. (33)
As showed in [START_REF] Nesterov | Squared functional systems and optimization problems[END_REF], the cone of non-negative univariate polynomials can be seen as the linear image of the cone of positive semi-definite matrices. Thus, finding non-negative polynomials over an infinite interval becomes a semi-definite programming problem.

Let P be a polynomial of degree 2n. Searching for p i , i = 1..2n + 1, the coefficients of the polynomial P, such that P(w) = ∑ 2n+1 i=1 p i w i-1 ≥ 0, ∀ w ∈ R is equivalent to searching for a symmetrical positive semi-definite matrix Y ∈ R (n+1)×(n+1) 0, such that [START_REF] Nesterov | Squared functional systems and optimization problems[END_REF]:

p i = tr(Y H n,i ), i = 1..2n + 1 ( 34
)
where H n,i ∈ R (n+1)×(n+1) are the Henkel matrices:

H n,i ( j, k) = 1, if j + k = i + 1 0, otherwise (35)
and tr denotes the trace of the matrix. The same principle is applied for the polynomial constraints (33). If there exist two symmetric positive semi-definite matrices Y x ∈ R (n+1)×(n+1) 0 and Y x ∈ R (n+1)×(n+1) 0 such that the coefficients of the polynomials Γ x and Γ x are equal to:

γ x = tr(Y x H n,1 ) ... tr(Y x H n,2n+1 ) T γ x = tr(Y x H n,1 ) ... tr(Y x H n,2n+1 ) T (36)
then constraints (33) are guaranteed to be satisfied continuously (∀w ∈ R).

The optimization problem (25) becomes a semi-definite program (SDP):

min Z,∆ Ṽ N ∑ i=1 Z i s.t.                                            -Z i ≤ ∆ Ṽi ≤ Z i Z i ≤ ∆ Ṽi , ∀ i = 1 . . . N M(ν N ) X(ν N ) = 0 Y x 0 Y y 0 Y z 0 Y x 0 Y y 0 Y z 0 γ x = [tr(Y x H 2,1 ) ... tr(Y x H 2,5 )] T γ x = [tr(Y x H 2,1 ) ... tr(Y x H 2,5 )] T γ y = [tr(Y y H 1,1 ) ... tr(Y y H 1,3 )] T γ y = [tr(Y y H 1,1 ) ... tr(Y y H 1,3 )] T γ z = [tr(Y z H 2,1 ) ... tr(Y z H 2,5 )] T γ z = [tr(Y z H 2,1 ) ... tr(Y z H 2,5 )] T (37)

B. Always feasible formulation

When solving an optimization problem like problem (37) as part of an autonomous rendezvous algorithm, it is important to guarantee at each step the feasibility of the solution. Infeasibility could arise from the fact that the objective cannot be reached in the N steps of the plan (because the control action is bounded by saturation constraints) or from the fact that no periodic trajectory could be found in the given tolerance region.

Infeasibility can be avoided by considering the dimensions of the tolerance box as optimization parameters. For the LP formulation of the problem, feasibility is guaranteed by using a scale factor for the tolerance box [START_REF] Tillerson | Co-ordination and control of distributed spacecraft systems using convex optimization techniques[END_REF]. In our case, each dimension of the tolerance box is considered as a separate parameter. As illustrated by equations (32), the coefficients of the polynomials that constitute the constraints depend linearly on X tol . So transforming X tol into a decision variable does not change the nature of the optimization problem and does not increase its complexity.

The optimization criterion must be modified to include the new variables:

min Z,∆ Ṽ ,X tol N ∑ i=1 Z i + ρ 3 ∑ i=1 X tol (i) (38)
where ρ > 0 is large enough to ensure that the tolerance box is modified only when needed to guarantee feasibility. Bounds should be given for the minimal box X m , to ensure that the optimization algorithm will not minimize the box and as a consequence increase the fuel consumption:

X tol ≥ X m (39)
IV. NUMERICAL EXAMPLES To illustrate our approach for finding the optimal solution to (18), we used a Prisma rendezvous mission, whose scenario information is given in Table . I. The problem was solved in Matlab using Yalmip [START_REF] Lofberg | YALMIP : a toolbox for modeling and optimization in MATLAB[END_REF] and SeDuMi [START_REF] Sturm | Using SeDuMi 1.02, A Matlab toolbox for optimization over symmetric cones[END_REF]. For comparison, we used the solution obtained from the LP approach, as presented in [START_REF] Breger | Cooperative Spacecraft Formation Flying: Model Predictive Control with Openand Closed-Loop Robustness[END_REF]. The LP problem was resolved with the Matlab's linear solver linprog, using 15 discretization points over 1 period for the tolerance box constraints.

The maneuvers plans obtained with each method yielded very close values for the optimization criterion. But the advantage of our methodology can be noticed when looking at the trajectories obtained after control, in Fig. 2. For the LP formulation, the box constraints are violated between some of the discretization instants. Our methodology produced trajectories that evolve very close to the bounds, without crossing them. It guarantees continuous satisfaction of the constraints, without adding extra inequalities to reach a certain precision.

For better understanding, Table II shows the influence on the performances of the LP based algorithm of the number of instants where constraints are explicitly checked. The same simulation data as for the first example is used. The discretization is done over 1 orbital period (T = 5843s), taking 15, 30 and then 50 instants. The given solver time (ST) is the mean time for 10 runs of each algorithm. The cost is the sum of ∆V (passing back to the original variables).

The costs obtained with each method are very similar. It can be easily observed that for the LP based algorithms increasing the number of discretization points increases the solver time and the number of constraints. On the other hand, it decreases the time spent outside the box (TOB), where the tolerance box constraints on the trajectory are violated. However, a large number of points would be necessary to get closer to the performances of the SDP approach (zero violation of the constraints).

A slightly different Prisma mission is used to illustrate the advantages of the always feasible formulation. The changes to the rendezvous scenario are given in Table III.

In the always feasible case, the obtained solution suggests that a box larger than X m is needed to ensure feasibility. As shown in Fig. 3, X tol = [1, 1, 33.0664] so an expansion on the z axis of the given tolerance region is necessary.

Used in an open-loop scheme, the always feasible formulation can provide some insight to the size of the smallest tolerance box that ensures feasibility. Integrated in a closedloop scheme, it can help avoiding unwanted situations where no solution is available.

V. CONCLUSIONS

In this paper a new method for obtaining an open-loop impulsive maneuvers plan leading to periodic proximity relative motion, between two spacecraft on Keplerian orbits was presented. The periodic relative trajectory is confined inside a specified tolerance region. By using the rational expressions for the propagation of the periodic motion and the result linking the cone of non-negative polynomials to the cone of semi-definite positive matrices, the requirement to define a prediction horizon where constraints are explicitly checked for a finite number of instants was removed. The method guarantees constraints satisfaction all along the path, without increasing the computational load or the fuel consumption. An always feasible formulation of the problem was also given, by considering the dimensions of the tolerance region as optimization parameters.

Fig. 2 .

 2 Fig. 2. xz axis projections of the periodic relative trajectory after control


David Skålid Amundsen

Julien Scheibert
email: julien.scheibert@ec-lyon.fr

Kjetil Thøgersen

Jørgen Trømborg

Anders Malthe-Sørenssen

D Skålid

Amundsen J K Scheibert

J Thøgersen

A Trømborg

Malthe-Sørenssen

D model of precursors to frictional stick-slip motion allowing for robust comparison with experiments

Keywords: Sliding friction, Stick-slip, Precursors, Spring-block model, Interfacial stiffness, Numerical simulation, Friction coefficient

We study the dynamic behaviour of 1D springblock models of friction when the external loading is applied from a side, and not on all blocks like in the classical Burridge-Knopoff-like models. Such a change in the loading yields specific difficulties, both from numerical and physical viewpoints. To address some of these difficulties and clarify the precise role of a series of model parameters, we start with the minimalistic model by Maegawa et al. (Tribol. Lett. 38 313, 2010) which was proposed to reproduce their experiments about precursors to frictional sliding in the stick-slip regime. By successively adding (i) an internal viscosity, (ii) an interfacial stiffness and (iii) an initial tangential force distribution at the interface, we manage to (i) avoid the model's unphysical stress fluctuations, (ii) avoid its unphysical dependence on the spatial resolution and (iii) improve its agreement with the experimental results, respectively. Based on the behaviour of this improved 1D model, we develop an analytical prediction for the length of precursors as a function of the applied tangential load. We also discuss the relationship between the microscopic and macroscopic friction coefficients in the model.

Introduction

The dynamics of frictional interfaces are crucial to many situations in mechanical engineering [1,[START_REF] Bhushan | Nanotribology and Nanomechanics[END_REF], geosciences [START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF] or biology [START_REF] Scheibert | [END_REF]5]. Today, after decades of studies, the science of contacts under time-invariant loading conditions, e.g. static contacts or steady sliding contacts [START_REF] Johnson | Contact Mechanics[END_REF][START_REF] Persson | Sliding Friction[END_REF][START_REF] Baumberger | [END_REF], has reached a high level of advancement which, in many instances, enables quantitative reproduction of global [START_REF] Baumberger | [END_REF]9] or local [START_REF] Scheibert | [END_REF]10,11] measurements. In contrast, the dynamics of contacts under rapidly evolving loads or during fast unstable motion like stick-slip [12][13][14] is far less understood. In particular, recent experiments on the transition from static to kinetic friction of side-driven poly(methyl methacrylate) (PMMA) rough samples forming line contacts with a rough PMMA substrate have revealed unexpected features [15][16][17][18]. The transition occurs through the fast (comparable to the speed of sound) propagation of micro-slip fronts through the contact [15,18]. It can also be preceded by a series of fronts that arrest before having ruptured the whole contact, thus denoted as precursors [16,17].

These results, which may have important implications for e.g. the study of earthquakes, have triggered an active modelling activity. Braun et al. [19], using a one-dimensional (1D) spring-block model with a complex time-dependent friction law, produced three types of micro-slip front velocities, analogous to that observed in [15]. Maegawa et al. [17], using a 1D spring-block model with a simple Amontons-Coulomb (AC) friction law, showed that the length of precursors is modified when the external normal load is made asymmetric. Scheibert and Dysthe [20], using a quasi-static 1D model with AC friction, showed how the increasing tangential load itself induces an increasing pressure asymmetry which influences the precursors' series. Due to the intrinsic limitations of 1D models, all these studies yielded only a qualitative agreement with experiments. Very recently, Trømborg et al. demonstrated, using a 2D spring-block model with AC friction, that quantitative agreement with the kinematics (i.e. the properties of the states in which no micro-slip front is propagating) of the experiments requires an accurate description of the interfacial stresses, and therefore the use of realistic boundary conditions on the sample [21].

Compared to such a 2D model, the strongest advantage of 1D models is that their results are much easier to analyse and understand, so that they provide opportunities for theoretical approaches. They also require simpler implementation and lower computational power. In most situations, in which very accurate results are not needed and/or a thorough qualitative understanding of the behaviour of the system is desired, 1D models are preferable. From the pioneering work of Burridge and Knopoff [22], spring-block models of friction have been extensively studied (see e.g. [23][24][25][26]). These models have mainly been used to describe the statistical properties of the series of earthquakes at a seismic fault. Inertial blocks are connected in series via internal springs that model the crust's elasticity. A homogeneous tectonic loading is modelled by coupling, via loading springs, all individual blocks to the same rigid driving body. Such statistical analysis of homogeneously driven systems contrast with the recent 1D studies by Braun et al. [19] and Maegawa et al. [17], in which the time evolution of side-driven systems is analysed deterministically, in order to produce data comparable to the experimental measurements. The two models by Braun et al. and Maegawa et al. are actually very different. Both consider an array of blocks connected by internal springs, but the model by Braun et al. [19] also considers viscous dissipation and a complex time-dependent friction law emerging from the collective behaviour of interfacial springs with random stiffness, breaking threshold and reattachment time. In contrast the model by Maegawa et al. [17], which only considers blocks and springs and the minimalistic AC friction law, is probably the simplest possible model.

On the one hand, we will see that, due to its extreme simplicity, the model by Maegawa et al. [17] yields results that are strongly resolution dependent, which prevents robust comparison with experiments. On the other hand, it is difficult to disentangle the respective roles of the many parameters of the model by Braun et al. [19]. The scope of this article is therefore to (i) construct, step by step, a minimal 1D side-driven spring-block model, the results of which are essentially resolution independent and (ii) qualify the capabilities of this model to reproduce the main qualitative features of recent experimental observations. This Letter is organised as follows. In section 2, we describe the model by Maegawa et al. [17] and show its limitations. In section 3, we improve this model by introducing successively an internal viscosity, an interfacial stiffness and an initial tangential force distribution. In section 4, we de-scribe an analytic prediction for the length of precursors in both Maegawa et al.'s and our improved model.

The model by Maegawa et al. and its limitations

2.1 The model by Maegawa et al. [17] In the model developed by Maegawa et al. [17], the slider is modelled as a chain of blocks connected by springs (Fig. 1), with material spring constant k and block mass m = M/N, where M is the total mass of the slider and N is the number of blocks. In experiments, the base (also called track) is fixed on a very stiff support, and it is therefore modelled as a rigid surface for simplicity. The tangential force is applied at the trailing edge of the system through a loading spring with stiffness K. One end of this spring is attached to the trailing edge block (block 1), while the other end of the spring moves at a constant velocity V . The normal force p i is imposed on each block, satisfying the criterion ∑ N i=1 p i = F N , where F N is the total applied normal force.

f 1 f 2 f N m m k K m V x p 2 p N p 1
Fig. 1 Schematics of the model system. The sample is modeled by N blocks of mass m connected in series through springs of stiffness k. The trailing edge of the system (block 1 ) is slowly driven through a loading spring of stiffness K. Each block is also submitted to a normal force p i and to a friction force f i .

The equations of motion are given by

m ün =    k(u 2 -u 1 ) + F T + f 1 , n = 1 k(u n+1 -2u n + u n-1 ) + f n , 2 ≤ n ≤ N -1 k(u N-1 -u N ) + f N , n = N, (1 )
where u n = u n (t) is the position of block n as a function of time relative to its equilibrium position and ¨denotes the double derivative with respect to t. F T = F T (t) is the driving force (or tangential force/load) given by

F T = K(Vt -u 1 ). (2)
A local friction law giving the friction forces f n is imposed between the blocks and the base. AC friction is used with local kinetic and static friction coefficients µ k and µ s , respectively. The resulting global friction coefficients are similarly denoted by µ K and µ S . The friction force on block n, f n , is then given by

f n = ≤ µ s p n , un = 0 -sgn( un )µ k p n , un = 0, (3)
where, when un = 0, equilibrium of block n imposes that f n balances all other forces acting on block n. The material spring constant k is chosen such that the elastic deformation of the model is similar to that of a linear elastic medium with Young's modulus E, which yields:

k = (N -1 )ES/L, (4)
where S and L are the cross-section area and the length of the slider, respectively.

In the experiments by Maegawa et al. [1 7], an asymmetric normal loading was used, leading to the following linear model for p n :

p n = F N N 1 - 2n -N -1 N -1 θ , (5)
where θ ∈ [-1 ,1 ] is a measure of the non-uniformity in the normal loading. The values of the parameters are chosen to be in agreement with their respective values in the experiments by Maegawa et al. [1 7], which are: K = 0.8 MN/m, V = 0.1 mm/s, F N = 400 N, M = 0.01 2 kg, L = 1 00 mm, S = 1 00 mm 2 , E = 2.5 GPa, µ s = 0.7 and µ k = 0.45.

Results and limitations of the model by Maegawa et al.

We have implemented the model by Maegawa et al. and tested our code by comparing it to (i) two other similar codes and (ii) the analytical solution of the equation of motion for a one-block system. The simulation starts with each spring at its equilibrium length, i.e. u n = 0 for all blocks. From t = 0, the force from the driving spring on block 1 is increased. Figure 2a shows that the resulting time evolution of F T (t) exhibits stick-slip behaviour. Each time a drop in the driving force is observed, some part of the slider moves relative to the base. We call the short time intervals during which movement occurs events.

Between events, no block is moving. Since the driving force is applied only at the trailing edge block, only block 1 is loaded and eventually reaches its static friction threshold, so that all events must nucleate at the trailing edge. The movement of block1 then loads block 2, which itself reaches its threshold and so on. This succession of blocks starting to move defines a micro-slip front, which propagates towards the leading edge, in analogy with the fronts observed in experiments [1 5, 1 7]. The distance from the trailing edge to the micro-slip front as a function of time, x f (t), is shown in Fig. 2b. If this front reaches the leading edge, the event is a global event, and the whole slider moves relative to the base.

From Fig. 2 it is evident that not all events are global: smaller events are observed between global events. In addition, a series of events with increasing maxima of x f (t) is seen to precede the first global event. These events occur for To perform a quantitative comparison between their experimental and numerical results, Maegawa et al. focused on the relationship between the normalized length L p /L of the series of precursors (Fig. 2b) and the normalized tangential force F T /F N at which they are triggered (Fig. 2a). To do this, they discarded all simulated precursors having a length smaller than any of the previous events, with the justification that no such smaller event was observed in the experiments. Note that 2D models produce series of precursors of monotonically increasing length, so that practically no event has to be discarded [21 ].

x f (t)/L t [s]
Figure 3 shows the tangential force τ normalised by the normal force p on each block at the time of initiation and arrest of a precursor of length L p = 0.6L. The tangential force is here defined as the total force on a block excepting the friction force. Each event is initiated when the tangential force on block 1 reaches the local static friction threshold. As block 1 moves, the tangential force on block 2 increases, eventually reaching the local static friction threshold, and starts to move. The precursor event arrests when the tangential force built on a block by its left neighbour is not sufficient to make it reach its static friction threshold. The slow loading of block 1 continues, and will eventually trigger a new event nucleating at the trailing edge.

The spatial resolution used in the above results, i.e. N = 1 0 as used by Maegawa et al. [1 7], is rather low. This is especially evident in Fig. 3, which only contains 1 0 data points for τ/ p. We expected to improve these results by simply increasing the resolution to N =1 00. The corresponding loading and front position curves, F T (t) and x f (t), are shown in Fig. 4. Fig. 3 Tangential force distribution at initiation of a microslip-event (•). Block 1 has reached its slip threshold (τ/p = µ s ). The event in- volves all blocks from 1 to 6 (i.e. a precursor of length L p /L = 0.6), leading to a modified tangential force distribution at arrest ( ). Results obtained using N = 1 0 (as in [1 7]) and θ = 0. One of the limitations of the model by Maegawa et al. is now clearly evident. By only changing the spatial resolution, the loading and front position curves are changed significantly. Some aspects are unchanged: precursors precede the first global event and then stick-slip behaviour is observed. The final average level of F T /F N also appears to be conserved. However, the amplitude of the drops in the loading curve is reduced, while the number of events, both global and precursory, is seen to increase significantly.

x f (t)/L t [s] F T (t)/F N b a 0 1 2 3
To illustrate this scaling with respect to the model resolution, we have plotted in Fig. 5 the evolution with N of (i) the total number of events and the number of global events (Fig. 5a), and (ii) the total number of precursors and the number of precursors longer than any previous ones (Fig. 5b). An approximately linear increasing trend is observed in all four curves. This behaviour is problematic as soon as one wants to compare to experiments, in which the size of the drops in F T /F N and the number of events are well-defined experimental measurements. Maegawa et al. [1 7] used N = 1 0, which produced a number of precursors similar to that observed in their experiments, but this agreement appears to be casual. A robust model should produce almost identical numbers of events whatever the spatial discretization of the slider. In this respect, note that 2D models do produce a resolution-independent number of events [21].

Another problematic behaviour of the model affects the tangential force spatial distribution, as shown in Fig. 6. The tangential force has been plotted at three different times: at t = 0.5 s, t = 3 s, and at the arrest of the precursor of length L p /L = 0.7. Strong oscillations are observed, with a half spatial period of the order of the lattice spacing, whatever the number of blocks N used. Again this unphysical dependence on the model resolution impedes a robust comparison with experimental measurements of the tangential stress distribution τ( x), like e.g. those of [1 8].

In their quantitative comparison between model and experimental results, Maegawa et al. [1 7] focused on the length of precursors L p as a function of the tangential load F T , L p (F T ). This relationship has also been studied experimentally by Rubinstein et al. [1 6]. While similar behaviours are found, we will only compare our model results to the experimental results of Maegawa et al..

The experimental setup by Maegawa et al. allowed for non-uniformities in the normal loading, and they studied its consequences on the length of precursors. The non-uniformity in the normal loading is modelled as an asymmetric distribution of the normal loads p n by using Eq. ( 5). The value of the parameter θ = ±0.833 is chosen to be in agreement with its corresponding experimental value.

The length of precursors L p is plotted as a function of the tangential force F T at event arrest in Fig. 7. Three dif- ferent values of θ are used, and the results for both N = 1 0 and N = 1 00 are included. As discussed above, the number of precursors increases with N. However, the shape of the curves does not change significantly with N, thus enabling comparison with the shape of the experimental curves. The qualitative behaviour of L p as a function of F T is in agreement with experiments. For θ = 0.833, the normal force on the trailing edge is reduced, leading to longer precursors for the same tangential load compared to θ = 0. On the other hand, θ = -0.833 leads to an increased normal load on the trailing edge, and precursors are shorter for a given F T . All curves converge to the same point at L p /L = 1 , meaning that the global static friction threshold µ S (the value of F T /F N when the first global event occurs) is independent of the normal loading distribution. The value of µ S appears to be approximately 0.45, which incidentally is the value of the local kinetic friction coefficient µ k . The reason for this will become clear in Sec. 4. A quantitative comparison between the experimental results in [1 7] and the model results in Fig. 7, however, reveals large discrepancies: all three experimental curves are found way below their simulated counterpart, meaning that the simulation strongly overestimates L p for any given F T ; the rapid increase in precursor length after L p /L ≈ 0.5 that is observed experimentally has no equivalent in the model.

L p /L = 0.7 t = 3 s t = 0.
Summing up, three main limitations of the model by Maegawa et al. have been identified: (i) the tangential force shows large oscillations, the wavelength of which is controlled by the lattice spacing, (ii) the number of all kinds of events is an increasing function of N and (iii) the quantitative agreement with the L p (F T )-curve between model and experiments is poor. In the following section, we propose improvements of the model that contribute to overcome these limitations.

3 Improvements of the model by Maegawa et al.

Introducing a relative viscous damping

Resolution dependent oscillations are known to occur in Burridge-Knopoff-like models and more generally in dynamic rupture models involving AC friction at the interface between dissimilar elastic media [27]. Classical ways to reduce them significantly are either to regularize the AC friction law (see e.g. [28] and references therein) or to introduce a viscous damping in the system [1 9, 21 , 29-31]. Both are physically sound, but we choose to adopt the second approach. Physically, such viscosity is a way to model the energy dissipation that any material undergoes during deformation. After Knopoff and Ni [31], we consider the following form for the viscous force F η n :

F η n =    η ( u2 -u1 ) , n = 1 η ( un+1 -2 un + un-1 ) , 2 ≤ n ≤ N -1 η ( uN-1 -uN ) , n = N, (6)
which is a damping on the relative motion of neighbouring blocks. As in e.g. [1 7, 21 ], we assume that energy dissipation due to the motion of a block relative to the substrate is satisfactorily included in f n , and therefore do not, in contrast to e.g. [1 9, 32], include any viscous damping at the interface in our system. This also serves to keep the model as simple as possible. The equations of motion are then given by

m ün =    k(u 2 -u 1 ) + F T + F η 1 + f 1 , n = 1 k(u n+1 -2u n + u n-1 ) + F η n + f n , 2 ≤ n ≤ N -1 k(u N-1 -u N ) + F η N + f N , n = N. ( 7
)
L p /L = 0. The tangential force τ n is still defined as the sum of all forces on a block excepting the friction force, i.e. now including F η n . It remains to choose a reasonable value of the damping coefficient η.

In App. A, the value of η for critical damping of waves of wavelength λ = 2a, corresponding to the cut-off wavelength, is calculated. The result is η c = √ km. It is also shown that the values of η corresponding to critical damping of higher wavelength oscillations are always larger than √ km. We want to damp out oscillations of wavelength λ = 2a, and using η = √ km is then a possibility. However, waves of other wavelengths close to λ = 2a will also be highly damped, causing significant changes to the dynamics. Since this is an undesirable effect, a compromise has to be made. As suggested in [31 ], the value

η = √ 0.1 √ km ≈ 0.32 √ km (8)
is used in the following. Note that since k ∝ N -1 and m ∝ 1 /N, η is N-independent for N ≫ 1 .

Figure 8 shows the tangential force distribution obtained when relative viscous damping is included. The improvement with respect to Fig. 6 is clear: the short wavelength oscillations have almost disappeared, resulting in a physically reasonable smooth tangential force profile. τ/p also appears to be on average approximately equal to µ k in rup- tured regions, a fact that will be utilised below to predict the precursor lengths. Some small one-or two-node spikes remain. They have a different origin, as they are mainly caused by the discreteness of the friction law: at the tip of the rupture, one block is moving and therefore increasing the force on its neighbour, which is still stationary. The region beyond this stationary block is therefore not affected by the approaching rupture front, and a spike will therefore appear at the rupture front. Spikes may also appear as an event arrests, also caused by one part of the interface slipping while another part is stuck.

In all the following, the viscous damping introduced in this section will be used.

Introducing a tangential stiffness of the interface

As discussed above, the model by Maegawa et al. exhibits an unphysical scaling with N. It is possible to understand this scaling by considering how the system is tangentially loaded. As stated above, the driving force only acts on block 1 . In order for an event to nucleate, this block has to reach its static friction threshold, which is proportional to the normal force: µ s p 1 ∝ p 1 ∝ 1 /N. Since the added driving force per time is independent of N, the time between two events will be proportional to 1 /N, and the frequency of events consequently scales as N, an argument which is fully consistent with Fig. 5. The origin of the odd N-dependence of the model by Maegawa et al. is therefore the decreasing size of the loading region as N is increased.

In a physical system, like in the experiments by Maegawa et al., the loading region has a well defined spatial extension, which is a combination of various effects. First, the tangential loading is applied at the trailing edge of the slider at some effective height h above the interface. As discussed experimentally in [1 6] and modelled in [21], such a loading condition makes the tangential stress very high in the region near the trailing edge, the extension of which is of order h. Second, the interface between slider and base is not rigid. Both surfaces are rough and the multi-contact layer between them has a finite tangential stiffness associated to the tangential deformation of the microasperities involved in the contact. Such stiffness can be measured experimentally [33] and is found much smaller than the slider's bulk stiffness. Such a low interfacial stiffness can be responsible for deviations with respect to AC friction in static contacts [34]. The interfacial stiffness also results in a physical finite size of the loading region, since a localized tangential force at the interface will induce tangential strains within the rough layer not only at the loaded point but also in its neighbourhood.

The first effect relates to the elastic coupling of points of the interface through the slider's bulk, which cannot be accounted for explicitly by a 1 D model like the one considered here. In contrast, the effect of the interfacial stiffness can be introduced [32] in a 1 D model in the following way. Each block is initially attached to the track by a spring with stiffness k t as seen in Fig. 9a. The spring connecting block n to the track has a breaking strength µ s p n (Fig. 9b). When the spring to the track is detached the block is subject to the kinetic friction force ±µ k p n (Fig. 9c). As the velocity of this block reaches zero, the spring reattaches such that the total force on the block is 0 at the time of reattachment (Fig. 9d).

The friction force is now given by

f n = -k t u n -u stick n if attached, -sign( un ) µ k p n if detached, (9)
where u stick n is the attachment point of the spring between block n and the track. It is given by

u stick n = u n - τ n k t , ( 1 0)
where u n and τ n are the position of and tangential force on block n at the instant of its last reattachment to the track. This causes the total force on block n to be zero at the time of reattachment. The spring then detaches at the time step in which one finds

-k t u n -u stick n > µ s p n . (11 )
The equations of motion are still given by Eqs. [START_REF] Persson | Sliding Friction[END_REF].

As the system is loaded tangentially, a finite region around the driving point is affected. This is illustrated in Fig. 1 0 where the tangential force at the time of nucleation of the first precursor is plotted. The length of this region depends on the stiffness k t of the springs between the blocks and the track.

Assuming (i) N ≫1 , (ii) the length of the loading region to be much smaller than L and (iii) slow loading compared to the internal dynamics of the system, we calculate analytically in App. B the tangential force profile at the time of nucleation of an event. No assumption is made for the tangential force profile at the time of arrest of the previous event, and the calculation is therefore valid for all events, not only for the first precursor shown in Fig. 1 The first block has just reached it threshold for slip (τ/p = µ s ). The tangential force decays exponentially according to Eq. ( 46). Results obtained using N = 1 00, θ = 0 and k t = 1 0 7 N/m.

is an exponential decay of the tangential force with a characteristic length l 0 , which consequently is a measure of the size of the loading region. From App. B, l 0 is related to k t by

l 0 = ESL Nk t . ( 1 2)
In this expression, we recognize Nk t to be the total stiffness, k tot t , of the interface (N springs of individual stiffness k t in parallel), which is a measurable quantity in a given experimental setup. We then obtain the relation

l 0 = ESL/k tot t , ( 1 3)
which shows that the size l 0 of the loading region is now independent of N, i.e. of the model's spatial resolution.

In the model by Maegawa et al., the simple rigid-plasticlike AC friction ruled the behaviour of the interface. Now that the interfacial stiffness is introduced, the friction law is elasto-plastic-like, an improvement that has often been defended as a necessary extension of AC [1 9, 33-36]. Our model is the simplest improvement of the model by Maegawa et al. that accounts for interfacial stiffness.

The envelopes of the time dependence of both the tangential load and length of precursors in this improved model are not changed significantly compared to those obtained using AC friction, as seen by comparing Figs. 2 and11 . However, comparing Figs. 11 to 4 shows a significant decrease in the number of events and a consequent increase in the amplitude of the drops in F T /F N . The shear force profiles are also very similar to that of Fig. 8, except that spikes now decay exponentially with a characteristic length l 0 . We will now check the model's behaviour with respect to its scaling with N. To do this, one has to choose the interfacial stiffness k t (or equivalently a loading zone size l 0 ). The value of k t could be calibrated using an experimental measurement of the total interfacial stiffness k tot t . Here, we do not have access to such a measurement, so we rather exploit the fact that the number of precursors longer than any previous one is controlled by the parameter l 0 . By trial and error, we found that l 0 = 5 mm, corresponding to k t = 1 0 7 N/m produced a number of such events similar to that observed in the experiments by Maegawa et al.. We will use this value of l 0 in the following.

x f (t)/L t [s]
In Fig. 1 2 the number of events using the elasto-plastic friction law is plotted in the same way as in Fig. 5 for AC friction. Both the number of global events and the number of precursors longer than any previous one are now seen to be approximately constant. The total number of events and the total number of precursors are, however, still increasing with N. This is mainly due to events involving one block, i.e. events in which no real front propagation occurs. In other words, interfacial friction is found to satisfactorily solve the resolution dependence of the model's results, provided one considers events that have a measurable length.

Introduction of an initial tangential force distribution

The introduction of an internal viscosity and an interfacial stiffness in the model by Maegawa et al. allowed us to obtain force distributons and numbers of micro-slip events that were physically sound. However, the predicted L p vs. F T curves still follow the same shapes as those shown in Fig. 7 which, as already mentioned, deviate significantly from those obtained experimentally. In an effort to further improve the 1 D model, we note that one of the main differences between the model and the experiment is the initial tangential force and total number of events (dashed line). (b) Number of precursors longer than any previous one (solid line) and total number of precursors (dashed line). Comparison with Fig. 5 shows that introduction of an interfacial tangential stiffness suppresses the resolution dependence of the numbers of global events and of precursors longer than any previous one. distribution (when no external tangential load has been yet applied). In the model, such forces are assumed to be zero all along the contact. However, in the experiments, both the slider and the base undergo different expansion rates during application of the normal loading. As discussed in [37], the associated slip at the interface is impeded by friction, thus yielding a significant tangential force distribution at the interface. Such distributions have been measured to be antisymmetric [1 8], in agreement with basic contact mechanics calculations [START_REF] Johnson | Contact Mechanics[END_REF], thus ensuring F T (0) = ∑ N n=1 τ n (0) = 0. This effect is again a bulk effect which is quantitatively reproduced in 2D models [21 ]. Here, in 1 D, we will only study the qualitative influence of an initial tangential force distribution on the length of precursors. We consider the simple linear distributions shown in Fig. 1 3a.

Implementation of an initial tangential force requires an initial relative displacement of the blocks. The initial tangential forces are given by

τ n (0) =    k(u 2 (0) -u 1 (0)) + F T (0), n = 1 k(u n+1 (0) -2u n (0) + u n-1 (0)), 2 ≤ n ≤ N -1 k(u N-1 (0) -u N (0)), n = N, ( 1 4)
and by choosing F T (0) = 0 and u 1 (0) = 0 the above equation can be rewritten to

u n (0) =    0, n = 1 u 1 + τ n /k, n = 2 2u n-1 -u n-2 + τ n-1 /k, n = 3, 4, . . . , N, ( 1 5)
thus enabling calculation of the initial positions of all blocks given τ n . The initial attachment position of the interfacial springs is then calculated using Eq. (1 0), which ensures that the total force on each block is zero at t = 0.

The length of precursors corresponding to the three initial tangential force profiles shown in Fig. 1 3a is shown in Fig. 1 3b, using θ = 0. With respect to a zero initial force distribution, the stronger the asymmetry the lower the L p vs. F T curve, i.e. the shorter the precursors for the same tangential force. Moreover, the slope of the curves is significantly increased at large forces. Both effects lead to a significantly improved qualitative agreement with the experimental results by both Maegawa et al. and Rubinstein et al. [1 6]. This shows that the initial force distribution, which arises naturally in 2D models [21 ], is a crucial parameter for the kinematics of precursors to sliding. Note that the number of precursors, which is closely related to the choice of l 0 (as mentioned in the previous subsection) is only weakly affected by the introduction of an initial tangential force distribution.

Analytical prediction of precursor lengths

In order to complete this study we derive an analytic prediction for the precursor length as a function of the tangential force at event arrest, L p vs. F T . We will first look at the simpler case of the model by Maegawa et al., and then extend the prediction to our improved model using both an interfacial stiffness and an initial tangential force profile.

Prediction in the model by Maegawa et al.

Assume that a precursor has reached block n p and has the length L p = (n p /N)L. We want to calculate the tangential force F T at the time of arrest of this event. At that time, all blocks are stuck, so that

F T = N ∑ n=1 τ n .
(

This means that, given the tangential force distribution at event arrest, the corresponding tangential force is found using Eq. (1 6). According to Figs. 6 and 8 the tangential force is observed to be approximately equal to the kinetic friction level from block 1 to n p , and 0 elsewhere. Using this assumption, Eq. (1 6) yields

F T = µ k n p ∑ n=1 p n = µ k F N N n p ∑ n=1 1 - 2n-N-1 N -1 θ , ( 1 7)
where Eq. ( 5) has been inserted for p n . If N ≫1 , the sum can be approximated by an integral and n replaced by x = nL/N, which yields

F T ≈ µ k F N N N L L p 0 1 - 2( xN/L) -N -1 N -1 θ dx, ( 1 8)
and approximating N ±1 ≈ N yields

F T ≈ µ k F N L L p 0 [1 -(2( x/L) -1) θ ] dx ( 1 9
)

F T ≈ µ k F N L p L 1 + θ 1 - L p L . ( 20
)
As seen in Fig. 7, this prediction is in very good agreement with our simulation results. The deviation between the actual precursors and the analytical curve is the result of a slightly incorrect assumed tangential force profile. Note that a similar good prediction scheme, numerical rather than analytical, was previously developed in the 2D study by Trømborg et al. [21 ].

According to Eq. ( 20) the global static friction coefficient µ S (µ S ≈ F T F N is easily evaluated for L p = L from Eq. ( 20)) is independent of θ and almost equal to the local kinetic friction coefficient µ k . This is in agreement with the model result in Fig. 7 and with the results of the 2D model by Trømborg et al. [21 ].

Prediction in our improved model

Prediction of the precursor length in our improved model including both an interfacial stiffness and an initial tangential force profile follows the same line as that for the model by Maegawa et al.. However, blocks may now move even though no event is occurring. Despite this, it is expected that block accelerations are small when all track springs are attached, which leads to the approximate validity of Eq. ( 1 6).

An approximate tangential force profile at the time of arrest of an event (Fig. 1 4) has to be found for a given L p . Again, blocks in [0, L p ] are assumed to have a tangential force equal to the kinetic friction level. Blocks in the interval [L p , L], however, now needs to be taken into account for two reasons: both the initial tangential force profile and the springs to the track lead to a non-zero tangential force for x > L p . The form of this profile in a static situation has been calculated in App. B, and is given in Eq. ( 46), where τ 0 ( x) now is the initial tangential force profile. However, it has to be modified to take into account that the loaded block is not located at x = 0, but at x = L p , and that this block does not have to be loaded up to the static friction level, but may take some other value, say α p n p , with α a coefficient to be defined. Our assumed tangential force profile at the arrest of a precursor of length L p is therefore given by τ

( x) = µ k p( x), x ∈ [0, L p ] α p(L p ) -τ 0 ( x) e - x-Lp l 0 + τ 0 (x), x ∈ [L p , L] (21 )
We have considered that all values of the amplitude of the peak at x = L p have the same probability to occur between µ k p and µ s p, so that we have given α its average value α = ( µ s + µ k )/2 in the predictions seen in Fig. 1 3b. Fig- ure 1 4 shows the assumed tangential force profile and the actual tangential force profile at the arrest of an event, and the agreement is seen to be satisfactory. From Eqs. (21), one can calculate F T (L p ) in the same way as was done for the model by Maegawa et al. This calculation is provided in Appendix C. As seen in Fig. 1 3b, the prediction scheme works well. Deviations between the actual precursors and the analytical curve have two contributions: incorrectly assumed tangential force profile and inertial effects, where the former gives the largest contribution. Again, the global static friction coefficient µ S is seen to be approximately equal to the local kinetic friction coefficient µ k . The prediction curves appear to bend slightly backwards at F T /F N ∼ 0.45. In our continuous prediction scheme, this corresponds to micro-slip fronts that are so close to the leading edge that all blocks at x > L p have a tangential force above the kinetic friction level, causing F T to get smaller as the front moves further. In reality, however, no precursor ever arrests in such a state but propagates all the way to the trailing edge.

Conclusion

Recent experimental results about the transition from static to kinetic friction in line contacts have triggered the study of the deterministic dynamics of 1 D spring-block friction models in which driving is applied at one extremity of the chain of blocks. In this Letter, we have improved the simplest of such models [1 7] in order to solve its intrinsic unphysical resolution-dependence and to ameliorate its qualitative agreement with experimental results on the kinematics of micro-slip fronts along the contact. In particular, the introduction of a tangential stiffness of the interface, by introducing a new length scale in the model, practically suppresses its resolution dependence and allows for reproduction of realistic numbers of precursory micro-slip fronts. The additional introduction of an initial tangential force distribution at the interface significantly improves the agreement with the evolution of the precursor length with the external tangential load obtained in experiments. Our improved model is intented to serve as a framework for more complex friction models when robust comparisons with experiments are desired.

We focused on1 D models because they are simple enough to enable deep insights into the qualitative effects of the model's parameters. However, it is known that 2D models [21 ] are required to provide quantitative agreement with experiments. In this respect, the improvements brought to the 1 D model are effective ways to account for intrinsically 2D effects: First, the length scale introduced through the interfacial stiffness enables coupling between remote points along the interface, analogous to the coupling through the slider's bulk; Second, the initial tangential force distribution accounts for the shear stress arising form the differential Poisson expansion of two bodies pressed together.

In analogy with 2D results we developed, based on the well defined force distribution left by an arrested precursor, an efficient analytical prediction for the precursors' length as a function of the external tangential load applied. We also find that, like in 2D, the macroscopic static friction coefficient of a side-driven contact is approximately equal to its microscopic kinematic friction coefficient.

Provided l 0 /L ≪ 1, the trailing edge will not be affected by the loading. The latter of these two boundary conditions yields τ(L) = ka 2 l 2 0 L 2 Ae L/l 0 + Be -L/l 0 + τ 0 (L) (42)

≈ ka 2 l 2 0 L 2 Ae L/l 0 + τ 0 (L) (43)
= τ 0 (L),

i.e. A = 0. The first boundary condition yields

τ(0) = ka 2 l 2 0 L 2 B + τ 0 (L) = µ s p 1 , (45)
and the shear force is therefore given by τ(x) = µ s p 1 -τ 0 (x) e -x/l 0 + τ 0 (x).

The characteristic length l 0 is given by Eq. ( 37), and inserting for k given by Eq. ( 4) and a yields

l 0 = k k t a = ELS (N -1)k t , (47)
and hence Eq. ( 12) for N >> 1.

Note that in a 3D situation, the exponential decay of the tangential stress with x would be replaced by a power law [19,[START_REF] Landau | Theory of Elasticity[END_REF].

C Derivation of the prediction of precursor lengths in our improved model

We start with Eq. ( 16) and use the assumed shear force profile in Eq. ( 21), shown in Fig. 14. Again we go to the limit N → ∞, resulting in the substitution (50)

We limit ourselves to predicting the precursors in Fig. 13b, i.e. using a tangential interfacial stiffness and a linear initial shear forces as depicted in Fig. 13a, but with θ = 0. The normal and initial shear force are then given by p(x) = p = F N /N = constant, (51)

τ 0 (x) = β p 2(x -L/2) L , (52)
where the parameter β determines the slope in the initial shear force profile. Inserting Eqs. The above integrals can be calculated easily, and the result is the tangential load F T as a function of the precursor length L p :

F T (L p ) = F N µ k L p L + 2β l 2 0 L 2 e - L-Lp l 0 -1 + β (L -L p ) L p L 2 + l 0 L β 1 + e - L-Lp l 0 -2 L p L + α 1 -e - L-Lp l 0 . (54)
We observe that again F T (L) = µ k F N .

F 6 Fig. 2

 62 Fig. 2 Time evolution of F T (a) and x f (b) in the model by Maegawa et al. using N = 1 0 (as in [1 7]) and θ = 0. Before macroscopic stick-slip (reached when t ≈ 2.8 s), the loading curve is punctuated by partial relaxations associated with precursors to sliding, i.e. micro-slip fronts spanning a length smaller than L.

Fig. 4

 4 Fig. 4 Time evolution of F T (a) and x f (b) in the model by Maegawa et al. using N = 1 00 and θ = 0. Comparison with Fig. 2 shows a drastic increase in the time-frequency of micro-slip events.

Fig. 5

 5 Fig. 5 Number of different kinds of events as a function of N, using θ = 0. (a) Number of global events (solid line) for t ∈ [5 s, 20 s] and total number of events (dashed line). (b) Number of precursors longer than any previous one (solid line) and total number of precursors (dashed line). All curves are stongly increasing functions of N, indicating unphysical resolution-dependence of the model's results.

Fig. 6 Fig. 7

 67 Fig.[START_REF] Johnson | Contact Mechanics[END_REF] Tangential force distribution τ n normalised by the normal force p n on each block, at three different times: at t = 0.5 s, t = 3 s, and at the arrest of the precursor of length L p /L = 0.7 using N = 1 00 and θ = 0. Unwanted large oscillations with a half-period equal to the lattice spacing appear.

7 t = 3 6 Fig. 8

 7368 Fig.[START_REF] Baumberger | [END_REF] The tangential force τ normalised by the normal force p on each block is plotted at three different times: at t = 0.5 s, t = 3 s, and at the arrest of the precursor of length L p /L = 0.7 from a simulation including the relative viscous damping and using N = 1 00, θ = 0 and η = √ 0.1 √ km. Comparison with Fig.6shows that lattice-controlled fluctuations have disappeared.

Fig. 9

 9 Fig.9 Sketch of the behaviour of the spring between block n and the track. (a) The static friction force on block n is exerted through a spring of stiffness k t attached to the track. (b) As the block is moving, the spring is stretched until it reaches its breaking strength µ s p n . (c) When the spring is broken, block n is subject to the kinetic friction force µ k p n . (d) As the block stops, the spring reattaches at x = u stick n such that the total force on the block is zero at the time of reattachment.

6 Fig. 10

 610 Fig.10 Tangential force τ normalised by the normal force p at initiation of the first precursor, when an interfacial stiffness is considered. The first block has just reached it threshold for slip (τ/p = µ s ). The tangential force decays exponentially according to Eq. (46). Results obtained using N = 1 00, θ = 0 and k t = 1 0 7 N/m.

F 6 Fig. 11

 611 Fig. 11 Time evolution of F T (a) and x f (b) in our improved model including both a relative viscous damping and an interfacial stiffness, using N = 1 00 and θ = 0. Comparison with Figs. 4 and 2 shows similar envelopes but very different numbers of events and amplitudes of the drops in F T /F N .

Fig. 12

 12 Fig. 12 Number of different kinds of events as a function of N, when the interfacial stiffness in taken into account. θ = 0, η = √ 0.1 √ km and l 0 = 5 mm. (a) Number of global events for t ∈ [5 s, 20 s] (solid line)and total number of events (dashed line). (b) Number of precursors longer than any previous one (solid line) and total number of precursors (dashed line). Comparison with Fig.5shows that introduction of an interfacial tangential stiffness suppresses the resolution dependence of the numbers of global events and of precursors longer than any previous one.

Fig. 13 (

 13 Fig. 13 (a) Three different antisymmetric initial tangential force distributions aimed at modelling the effect of friction-frustrated differential Poisson expansion during normal loading. •: homogeneous distribution. × and +: linear profiles of increasing slope (β = 0.225 and 0.45 in Eq. (52) respectively). (b) The length of precursors corresponding to the three initial tangential force profiles shown in (a). Results obtained using N = 1 00, θ = 0, η = √ 0.1 √ km and l 0 = 5 mm. Solid lines are the analytical predictions of the precursor lengths discussed in Sec. 4.

6 Fig. 14

 614 Fig.14 The tangential force profile at the arrest of the 1 4th precursor plotted in Fig.13 (red dots) and the predicted profile used in the analytical prediction of precursor lengths (black drawn line) using N = 1 00, θ = 0, η = √ 0.1 √ km, l 0 = 5 mm and β = 0.225 in Eq. (52).

µ

  dx, n → xN/L. (48)The tangential force after a precursor of length L p is then given by k p(x) dx+ L Lp α p(L p ) -τ 0 (x) e -x-Lp l 0 + τ 0 (x) dx .

  (51) and (52) into Eq. (50) yieldsF T = N L Lp 0 µ k p dx + L Lp α p -β p 2(x -L/2) L e -x-Lp l 0 + β p 2(x -L/2) L dx . (53)

Acknowledgements We thank J. L. Vinningland for discussions. We acknowledge funding from the European Union (Marie Curie Grant No. PIEF-GA-2009-237089). This article was supported by a Center of Excellence grant to PGP from the Norwegian Research Council.

A Relative viscous damping in a linear chain of blocks

If friction forces are ignored, the equation of motion for an infinite chain of blocks connected by springs is given by m ün = k(u n+1 -2u n + u n-1 ) + η( un+1 -2 un + un-1 ).

We then assume a solution of the form

where ζ κ ∈ C and κ ∈ R. Inserting Eq. ( 23) into Eq. ( 22) yields the relation

which can be simplified to

since

The complex parameter ζ κ is then given by

The system is critically damped when Eq. ( 25) only has one solution for ζ κ , which occurs when the square root is zero:

The oscillations that are to be reduced have a wavelength λ = 2a, i.e. a wave number κ = 2π/λ = π/a. Inserting this into Eq. ( 28) leads to

which is the value of the damping coefficient η for which waves of wavelength λ = 2a are critically damped. Since the absolute value of sin is always smaller than one, choosing η = √ km will cause all other waves to be under-damped.

B Tangential force profiles and characteristic length with a tangential stiffness of the interface

An analytical expression for the characteristic length l 0 can be found. In order to do so, the following assumptions are made: N ≫ 1, l 0 /L ≪ 1 and slow loading compared to the internal dynamics of the system, which enables a static analysis. The system is first placed in a static state with an initial shear force profile given by τ 0 n , and then loaded slowly from the left. The equilibrium of all non-edge blocks writes

We introduce a new variable u ′ n defined by

where u 0 n is the initial position of block n. Inserting Eq. ( 31) into Eq. ( 30) yields

where

The two terms τ 0 n and -k t u 0 nu stick n cancel in Eq. ( 32) since the initial state is static, and thus

The above equation can be rewritten to

where a = L/(N x -1) is the lattice spacing. Since N ≫ 1, the first term in Eq. ( 35) can be replaced by the second spatial derivative, and replacing u ′ n with u ′ (na) = u ′ (x) yields

which has the general solution

The shear force is given by

By replacing again finite differences with second order derivatives,

and the general expression for the shear force profile can be found by using Eq. ( 37), which yields τ(x) = ka 2 l 2 0 L 2 Ae x/l 0 + Be -x/l 0 + τ 0 (x).

(41

The system is loaded from the left, and at the beginning of an event the shear force on block 1 is equal to the static friction threshold µ s p 1 .