D Efimov

R Ushirobira

J A Moreno

W Perruquetti

Homogeneous Lyapunov functions: from converse design to numerical implementation

The problem of the synthesis of a homogeneous Lyapunov function for an asymptotically stable homogeneous system is studied. First, for systems with nonnegative degree of homogeneity, several expressions of homogeneous Lyapunov functions are derived, which depend explicitly on the supremum or the integral (over finite or infinite intervals of time) of the system solutions. Second, a numeric procedure is proposed, which ensures the construction of a homogeneous Lyapunov function. The analytical results are illustrated by simulations.

I. INTRODUCTION

Stability analysis is one of the central problems in the domains of control and dynamical systems [START_REF] Bacciotti | Liapunov Functions and Stability in Control Theory[END_REF], [START_REF] Bhatia | Dynamical Systems: Stability Theory and Applications[END_REF], [START_REF] Chetaev | The Stability of Motion[END_REF], [START_REF] Fradkov | Pogromsky. Introduction to oscillations and chaos[END_REF], [START_REF] Hahn | Stability of Motion[END_REF], [START_REF] Isidori | Nonlinear control systems: An Introduction[END_REF], [START_REF] Krasovskii | Stability of Motion: Applications of Lyapunov's Second Method to Differential Systems and Equations With Delay[END_REF], [START_REF] Lyapunov | The general problem of the stability of motion[END_REF], [START_REF] Sepulchre | Constructive nonlinear control[END_REF]. The main approach to check the stability of a nonlinear dynamical system consists in applying the Lyapunov function method, that gives necessary and sufficient conditions for stability. The main issue with the application of this method is that there exists no procedure to design such a Lyapunov function for a generic dynamical system. There exist analytical methods to construct a Lyapunov function for linear or Lurie systems [START_REF] Khalil | Nonlinear Systems[END_REF], and several numerical approaches have been proposed for linear and partially linear (linear parameter-varying) systems [START_REF] Briat | Linear Parameter-Varying and Time-Delay Systems[END_REF], see also the survey [START_REF] Giesl | Review on computational methods for Lyapunov functions[END_REF]. The sum of squares approach also can be effectively applied for some classes of systems (like polynomial ones) [START_REF] Peet | A converse sum of squares lyapunov result with a degree bound[END_REF].

Homogeneous systems form a subclass of nonlinear systems [START_REF] Bacciotti | Liapunov Functions and Stability in Control Theory[END_REF], [START_REF] Bernuau | On homogeneity and its application in sliding mode control[END_REF], that include linear systems as a special case. These systems inherit many useful properties in a generic setting from the linear case: scalability of solutions, estimation of the rate of convergence and robustness [START_REF] Bernuau | Verification of ISS, iISS and IOSS properties applying weighted homogeneity[END_REF], [START_REF] Bernuau | On homogeneity and its application in sliding mode control[END_REF]. A remarkable fact about homogeneous systems is that they possess homogeneous Lyapunov functions [START_REF] Zubov | Methods of A.M. Lyapunov and Their Application[END_REF], [START_REF] Rosier | Homogeneous lyapunov function for homogeneous continuous vector field[END_REF]. However, no numeric routine has been proposed so far to construct such a Lyapunov function using homogeneity properties (the works [START_REF] Zubov | Methods of A.M. Lyapunov and Their Application[END_REF], [START_REF] Rosier | Homogeneous lyapunov function for homogeneous continuous vector field[END_REF] assumed that a Lyapunov function is already given, and a special transformation is then presented providing homogeneity property), and the present paper is going to fill this gap.

The aim of this work is to propose several homogeneous Lyapunov functions based on converse Lyapunov function methods, and also to develop numerical tools for calculating an explicit homogeneous Lyapunov function for a given nonlinear homogeneous system. The proposed results are based on an important feature of homogeneous Lyapunov functions: their construction and analysis can be performed on a sphere only, making calculations simpler. The numerical algorithms obtained here for constructing homogeneous Lyapunov functions can be used for checking the stability of homogeneous systems (since a system on a sphere can be approximated by a homogeneous one [START_REF] Efimov | Conditions of oscillations and multi-homogeneity[END_REF], this approach can be applied to a wide class of nonlinear systems locally). Their extension can also be implemented for control design using universal formulas [START_REF] Sontag | A "universal" construction of Arstein's theorem on nonlinear stabilization[END_REF], which are implicitly dependent on (control) Lyapunov functions.

The outline of this paper is as follows. The notation and preliminary results are introduced in sections II and III. The analysis of robustness of homogeneous Lyapunov functions is carried out in Section IV. Several explicit analytic expressions for homogeneous Lyapunov functions are proposed in Section V. Numerical routines using those expressions to design a homogeneous Lyapunov function are presented in Section VI. The proposed approach for constructing homogeneous Lyapunov functions is illustrated by numerical examples in Section VII. Concluding remarks and discussion appear in Section VIII.

II. NOTATION

• R + = {x ∈ R : x ≥ 0}, where R is the set of real numbers, N is the set of natural numbers.

• | • | denotes the absolute value in R, . denotes the Euclidean norm on R n (for n ∈ N), x A = inf ξ∈A x -ξ is the distance from a point x ∈ R n to a set A ⊂ R n . • A continuous function α : R + → R + belongs to the class K if α(0) = 0 and the function is strictly increasing.

The function α : R + → R + belongs to the class K ∞ if α ∈ K and it increases to infinity. A continuous function

The preliminary version of this work has been partially presented in [START_REF] Efimov | On the numerical construction of homogeneous lyapunov functions[END_REF] without proofs and with a different numeric procedure. D. Efimov, R. Ushirobira and W. Perruquetti are with Inria, Univ. Lille, CNRS, UMR 9189 -CRIStAL, Centrale Lille, F-59000 Lille, France. J. A. Moreno is with Instituto de Ingeniería, UNAM, Coyoacán, 04510 México DF, Mexico. This work was partially supported by HoTSMoCE Inria Associate Team and by ANR 15 CE23 0007 (Project Finite4SoS).

β : R + × R + → R + belongs to the class KL if β(•, t) ∈ K ∞ for each fixed t ∈ R + and lim t→+∞ β(s, t) = 0 for each fixed s ∈ R + .

• A sequence of integers 1, 2, ..., n is denoted by 1, n.

III. PRELIMINARIES

Consider the following nonlinear system:

ẋ(t) = f (x(t)), t ≥ 0, (1)
where x(t) ∈ R n is the state and f : R n → R n ensures forward existence and uniqueness of the system solutions at least locally in time, f (0) = 0. For an initial condition x 0 ∈ R n denote the corresponding solution by X(t, x 0 ) for any t ≥ 0 for which the solution exists (the solutions are understood in the Carathéodory sense). A set A ⊂ R n is called forward invariant for (1) if x 0 ∈ A implies that X(t, x 0 ) ∈ A for all t ≥ 0.

A. Stability definitions

Following [START_REF] Bacciotti | Liapunov Functions and Stability in Control Theory[END_REF], [START_REF] Khalil | Nonlinear Systems[END_REF], [START_REF] Lin | A smooth converse Lyapunov theorem for robust stability[END_REF], [START_REF] Polyakov | Nonlinear feedback design for fixed-time stabilization of linear control systems[END_REF], let Ω be an open neighborhood of a forward invariant set A ⊂ R n of (1).

Definition 1. For the system (1), the set A is said to be: (a) Lyapunov stable if for any x 0 ∈ Ω the solution X(t, x 0 ) is defined for all t ≥ 0, and for any > 0 there is δ > 0 such that for any x 0 ∈ Ω, if x 0 A ≤ δ then X(t, x 0 ) A ≤ for all t ≥ 0;

(b) asymptotically stable if it is Lyapunov stable and for any κ > 0 and > 0 there exists T (κ, ) ≥ 0 such that for any x 0 ∈ Ω, if x 0 A ≤ κ then X(t, x 0 ) A ≤ for all t ≥ T (κ, );

(c) finite-time stable if it is Lyapunov stable and finite-time converging from Ω, i.e. for any x 0 ∈ Ω there exists 0 ≤ T < +∞ such that X(t, x 0 ) ∈ A for all t ≥ T . The function

T A (x 0 ) = inf{T ≥ 0 : X(t, x 0 ) ∈ A ∀t ≥ T } is called the settling time of the system (1); (d) fixed-time stable if it is finite-time stable and sup x0∈Ω T A (x 0 ) < +∞.
The set Ω is called the domain of stability/attraction.

If Ω = R n , then the corresponding properties are called global Lyapunov/asymptotic/finite-time/fixed-time stability of (1) at A.

For a V : R → R denote the upper Dini derivative:

V + (t) = lim sup h→0 + V (t + h) -V (t) h , ∀ t ∈ R + .
If V is locally Lipschitz continuous then V + is finite, and if V is differentiable then V + is the usual derivative of V . For V : R n → R the generalized directional derivative at x ∈ R n in the direction d ∈ R n is defined by:

D + V (x)d = lim sup y→x h→0 + V (y + hd) -V (y) h .

B. Weighted homogeneity

Following [START_REF] Bacciotti | Liapunov Functions and Stability in Control Theory[END_REF], [START_REF] Zubov | Methods of A.M. Lyapunov and Their Application[END_REF], for strictly positive real numbers r i (i = 1, n) called weights and λ > 0, one can define:

• the vector of weights r = (r 1 , . . . , r n ) T , r max = max 1≤j≤n r j and r min = min 1≤j≤n r j ;

• the dilation matrix function Λ r (λ) = diag (λ ri ) n i=1 (note that ∀x ∈ R n and ∀λ > 0 we have Λ r (λ)x = (λ r1 x 1 , . . . , λ ri x i , . . . , λ rn x n ) T );

• the r-homogeneous norm of x ∈ R n is x r = ( n i=1 |x i | ρ r i ) 1
ρ for ρ ≥ r max (it is not a norm in common sense since the triangle inequality does not hold);

• for s > 0 the sphere and the ball in the homogeneous norm S r (s) = {x ∈ R n : x r = s} and B r (s) = {x ∈ R n :

x r ≤ s}, respectively.

Definition 2. A function g : R n → R is r-homogeneous of degree µ ∈ R if ∀x ∈ R n and ∀λ > 0 we have:

λ -µ g(Λ r (λ)x) = g(x).
A vector field f : R n → R n is r-homogeneous of degree ν ∈ R, with ν ≥ -r min if ∀x ∈ R n and ∀λ > 0 we have:

λ -ν Λ -1 r (λ)f (Λ r (λ)x) = f (x)
, which is equivalent to the i-th component of f being a r-homogeneous function of degree r i + ν.

The system (1) is r-homogeneous of degree ν if the vector field f is r-homogeneous of degree ν.

For r-homogeneous system (1) of degree ν, the solutions also admit a kind of homogeneity as functions of time [START_REF] Bacciotti | Liapunov Functions and Stability in Control Theory[END_REF], i.e. ∀x 0 ∈ R n and ∀t ∈ R:

X(t, Λ r (λ)x 0 ) = Λ r (λ)X(λ ν t, x 0 ) ∀λ > 0.
By its definition, • r is an r-homogeneous function of degree 1, and there exist σ, σ ∈ K ∞ such that

σ( x r ) ≤ x ≤ σ( x r ) ∀x ∈ R n .
Consequently, due to this equivalence of the norms • and • r , the homogeneous norm can be used in Definition 1. A variant of selection of σ, σ for a particular case considered in this work is given below:

Proposition 1. If the restrictions r max ≤ 1 and ρ ≥ 1 are satisfied, then

σ -1 (s) = n 1 ρ s 1 rmax s ≤ 1 s 1 r min s > 1 , σ(s) = n ρ-1 2 s rmin s ≤ 1 s rmax s > 1 . Proof: Under introduced restrictions ρ r i ≥ 1 ∀i = 1, n.
Using Jensen's inequality for ρ ≥ 1 and any a i ∈ R + for all i = 1, n, since the power function is concave, we obtain:

n i=1 a i 1 ρ = n 1 ρ n i=1 a i n 1 ρ ≥ n 1 ρ -1 n i=1 a 1 ρ i ,
from which the following inequality is satisfied:

x r = n i=1 |x i | ρ r i 1 ρ ≥ n 1 ρ -1 n i=1 |x i | 1 r i ≥ n 1 ρ -1 |x i | 1 r i .
Therefore,

x = n i=1 |x i | 2 ≤ n ρ-1 n i=1 x 2ri r ≤ n ρ-1 2 x rmin r x r ≤ 1 x rmax r x r > 1 ,
and we derived the desired expression of σ. Similarly

x r = n i=1 |x i | ρ r i 1 ρ ≤ n i=1 x ρ r i 1 ρ ≤ n 1 ρ x 1 rmax x ≤ 1 x 1 r min x > 1 , (2)
and inverting the above inequality we have an expression of σ.

If system (1) is r-homogeneous with degree ν, then a direct computation shows that it is also r-homogeneous for r = r -1 max (r 1 , . . . , r n ) T with degree ν = ν rmax . Therefore, the conditions of Proposition 1 can be always satisfied for any homogeneous system with r max = ρ = 1.

Corollary 1. Let r max ≤ 1 and ρ ≥ 1, then • r is locally Lipschitz continuous.

Proof: The result is a direct consequence of (2) and the fact that the power function x → x α for α ≥ 1 and the norm function • are locally Lipschitz continuous (a composition of Lipschitz continuous functions inherits the same property).

Further in this work we will always assume without loosing generality that the conditions r max ≤ 1 and ρ ≥ 1 are satisfied in order to use the established above continuity property of the homogeneous norm • r .

C. Time of convergence in homogeneous systems

An important advantage of homogeneous systems is that their rate of convergence can be evaluated qualitatively based on their degree [START_REF] Nakamura | Smooth lyapunov functions for homogeneous differential inclusions[END_REF]: 1) is r-homogeneous of degree ν and asymptotically stable at the origin, then it is (i) globally finite-time stable at the origin if ν < 0;

Theorem 1. If (
(ii) globally exponentially stable at the origin if ν = 0;

(iii) globally fixed-time stable with respect to the unit ball B r (1) if ν > 0.

Following that, we can further quantify the system convergence using the definition below.

Definition 3. For q > 1 define T q : R n → R + such that

X(T q (x 0 ), x 0 ) r = q -1 x 0 r , i.e. it is the function of contraction in q times.
From this definition, it is easy to show that T q is r-homogeneous with degree -ν (the degree of homogeneity for (1)):

T q (Λ r (λ)x 0 ) = λ -ν T q (x 0 ), ∀x 0 ∈ R n , ∀λ > 0,
and T q (x 0 ) is uniquely defined for each x 0 , but it can be set-valued. For a r-homogeneous system (1) of degree ν, if it is asymptotically stable at the origin, then there exists a continuously differentiable and r-homogeneous function V : R n → R + of degree µ > -ν such that for all x ∈ R n :

a x µ r ≤ V (x) ≤ b x µ r , (3)
D + V (x)f (x) ≤ -cV 1+ ν µ (x), (4)
for some 0 < a ≤ b and c > 0 [START_REF] Zubov | Methods of A.M. Lyapunov and Their Application[END_REF], [START_REF] Rosier | Homogeneous lyapunov function for homogeneous continuous vector field[END_REF]. Then for any

x 0 ∈ R n , X(t, x 0 ) r ≤ β( x 0 r , t) ∀t ≥ 0, (5)
where

β(s, t) =                    ( b a ) 1 µ s 1+a ν µ ν µ cs ν t 1 ν ν > 0      ( b a ) -ν µ s -ν + ν µ ct -1 ν t < - ( b a ) -ν µ s -ν ν µ c 0 t ≥ - ( b a ) -ν µ s -ν ν µ c ν < 0 ( b a ) 1 µ e -ct µ s ν = 0
is a function from class KL, and it represents a generic parameterization of the upper bound β( x 0 r , t) of X(t, x 0 ) r for all t ≥ 0 in homogeneous systems.

Definition 4. For q > 1, a continuous function T q : R n → R + such that X(t, x 0 ) r ≤ q -1 x 0 r for all t ≥ T q (x 0 ) for any x 0 ∈ R n is called an upper bound function of contraction in q times.

Taking into account the expression of β given in (5), we can conclude that a possible selection is

T q (x 0 ) = T q ( x 0 r ) with for s > 0, T q (s) =            q ν ( b a ) ν µ -1 a ν µ ν µ c s -ν ν > 0 ( b a ) -ν µ -q ν -ν µ c s -ν ν < 0 µ c ln q( b a ) 1 µ ν = 0 , (6)
which is uniform on the spheres and skipping some ambiguity, the same symbol T q is used for compactness of notation. Then T q (x 0 ) can also be considered r-homogeneous with degree -ν:

T q (Λ r (λ)x 0 ) = λ -ν T q (x 0 ), ∀ x 0 ∈ R n , ∀ λ > 0.

IV. ON HOMOGENEOUS PERTURBATIONS OF HOMOGENEOUS LYAPUNOV FUNCTIONS

This section is devoted to the analysis of the influence of homogeneous additive perturbations on homogeneous Lyapunov functions, and the establishment of conditions providing that such a perturbation leads to another homogeneous Lyapunov function.

To this end, let us assume that for a r-homogeneous system (1) of degree ν, there exists a locally Lipschitz continuous and r-homogeneous function V : R n → R + of degree µ > -ν, then the analysis and the verification of all properties of V can be performed on the sphere S r (1) only. Indeed, for any x ∈ R n there is y ∈ S r (1) such that x = Λ r ( x r )y, then

a x µ r ≤ V (x) = x µ r V (y) ≤ b x µ r ,
where a = inf y∈Sr(1) V (y) and b = sup y∈Sr(1) V (y) (same as in ( 3)), and

D + V (x)f (x) = x ν+µ r D + V (y)f (y) ≤ -d x ν+µ r (7)
where

d = -inf y∈Sr(1) D + V (y)f (y) (c = b -1-ν µ d in (4)
). Now, for such a function V let us show that there is a family of locally Lipschitz continuous and r-homogeneous Lyapunov functions V : R n → R + of degree µ sufficiently close to V . To this end, denote ε(x) = V (x) -V (x) as the difference between these functions. Clearly ε is a r-homogeneous function of degree µ. We wish to formulate the restrictions on ε so that V (x) = V (x) + ε(x) would be another Lyapunov function for (1) inheriting the properties of V .

Proposition 2. For a locally Lipschitz continuous and r-homogeneous function V : R n → R + of degree µ > -ν, assume that the estimates (3), ( 7) are satisfied for constants 0 < a ≤ b and d > 0. Let ε : R n → R be a locally Lipschitz continuous and r-homogeneous function of degree µ such that Proof: Obviously, V is locally Lipschitz continuous and r-homogeneous as a sum of two functions possessing these properties. Next again, for any x ∈ R n there is y ∈ S r (1) such that x = Λ r ( x r )y. By definition and by homogeneity:

-a < ε, d > ε, where ε = inf y∈Sr(1) ε(y) and ε = sup y∈Sr(1) D + ε(y)f (y). Then V (x) = V (x) + ε(x)
V (x) = x µ r V (y) = x µ r (V (y) + ε(y))
and there exist a = inf y∈Sr(1) V (y) and b = sup y∈Sr(1) V (y). If 0 < a + ε then the relations 0 < a ≤ b are satisfied, and for all

x ∈ R n a x µ r ≤ V (x) ≤ b x µ r .
Similarly,

D + V (x)f (x) = x ν+µ r D + V (y)f (y) = x ν+µ r (D + V (y)f (y) + D + ε(y)f (y))
and there exists d = -inf y∈Sr(1)

D + V (y)f (y). If -inf y∈Sr(1) D + V (y)f (y) = d > ε = sup y∈Sr(1) D + ε(y)f (y), then d > 0.
Remark 1. To demonstrate that the restriction d > ε is always feasible for some ε, an upper estimate of ε can be used:

ε = sup y∈Sr(1) ∂ + ε(y) sup y∈Sr(1) f (y) ,
where ∂ + ε(y) is the generalized gradient defined in [START_REF] Clarke | Nonsmooth Analysis and Control Theory[END_REF]. Then for any f in (1), by minimizing the values of ε(y) and ∂ + ε(y) , it is possible to ensure that all conditions of Proposition 2 are satisfied.

V. CONSTRUCTING HOMOGENEOUS LYAPUNOV FUNCTIONS

In this section several explicit formulas for locally Lipschitz continuous and homogeneous functions are proposed completing the results of [START_REF] Zubov | Methods of A.M. Lyapunov and Their Application[END_REF], [START_REF] Rosier | Homogeneous lyapunov function for homogeneous continuous vector field[END_REF], where only the existence of smooth homogeneous Lyapunov functions has been substantiated without an expression that can be used for numerical synthesis. For this purpose, in this section we will assume that a homogeneous system (1) is asymptotically stable and the estimate (5) is given. Assumption 1. Let (1) be asymptotically stable at the origin with a locally Lipschitz continuous and r-homogeneous vector field f of degree ν ≥ 0.

A possible expression for the function β ∈ KL is given just below [START_REF] Briat | Linear Parameter-Varying and Time-Delay Systems[END_REF] and the parameters of the function β, i.e. a, b, c and µ, are assumed to be fixed (as usual for the converse results). It is worth to stress that for the analysis in the sequel, these parameters are not related with a Lyapunov function of (1), they come from a possible parametrization of β given after (5). Remark 2. Recall that under Assumption 1, if the system (1) is locally Lipschitz continuous, then for any compact set of initial conditions E ⊂ R n and any time 0 T < +∞, there exists K E,T ∈ R + such that [START_REF] Khalil | Nonlinear Systems[END_REF] (assuming that the solutions originated in E are defined on [0, T ])

X(t, x 1 ) -X(t, x 2 ) K E,T x 1 -x 2
for all 0 t T and all x 1 , x 2 ∈ E. Using the Lipschitz continuity of • r established in Corollary 1 and the fact that for an asymptotically stable system X(t, x 0 ) ∈ E for all t ≥ 0 and x 0 ∈ E for some compact

E ⊂ R n ( E = B r (( b a ) 1
µ
) for E = B r (1) according to (5)), then there exists

L E,T ∈ R + such that | X(t, x 1 ) r -X(t, x 2 ) r | L E,T x 1 -x 2
for all 0 t T and any x 1 , x 2 ∈ E. These constants K E,T and L E,T will be used below in the proofs.

A. Design based on supremum of trajectories

First, the converse Lyapunov function method initiated in the works of [START_REF] Persidskii | On a theorem of Liapunov[END_REF], [START_REF] Yoshizawa | On the stability of solutions of a system of differential equations[END_REF] is used, which is developed below for the needs of homogeneous systems.

Lemma 1. Let Assumption 1 be satisfied. Consider

V (x) = sup t 0 {v(X(t, x)) k(v(x) ν t)}, v(x) = sup t 0 X(t, x) r ,
where k : R + → R + is a continuously differentiable function satisfying 0 < κ 1 k(t) κ 2 < +∞ with a monotonically decreasing function k(t) > 0 such that k(t)t ≤ κ 3 < +∞ for all t 0. Then V is a locally Lipschitz continuous and r-homogeneous of degree 1 function satisfying

κ 1 x r V (x) for all x ∈ R n , D + V (x)f (x) < 0 for a.a. x ∈ R n \ {0}.
An example of the function k(t) needed in this lemma is

k(t) = (κ 1 + κ 2 t)(1 + t) -1 , k(t) = (κ 2 -κ 1 )(1 + t) -2 , κ 2 > κ 1 > 0. Proof: For any x 0 ∈ R n , define v(x 0 ) = sup t 0 X(t, x 0 ) r ,
by construction x 0 r v(x 0 ) and v(0) = 0. From the attractivity of the origin and the continuity of X(t, •), for any x 0 ∈ R n there exists a finite time T x0 ∈ R + such that v(x 0 ) = sup 0 t Tx 0 X(t, x 0 ) r . For a homogeneous system (1), a useful and simple choice is T x0 = T q (x 0 ) for some properly selected q > 1 (the function of contraction in q times from Definition 3). To analyze continuity of the function v, consider

|v(x 1 ) -v(x 2 )| = | sup t 0 X(t, x 1 ) r -sup t 0 X(t, x 2 ) r | = | sup 0 t Tx 1 X(t, x 1 ) r -sup 0 t Tx 2 |X(t, x 2 ) r | sup 0 t T | X(t, x 1 ) r -X(t, x 2 ) r |,
where T = max{T x1 , T x2 } and x 1 , x 2 ∈ R n . For any compact E ⊂ R n there exists T E = sup x0∈E T x0 with the property T E < +∞. Keeping this in mind, and Lipschitz continuity of the system (1) (Remark 2), we see that

|v(x 1 ) -v(x 2 )| sup 0 t T E | X(t, x 1 ) r -X(t, x 2 ) r | L E,T E x 1 -x 2
for all x 1 , x 2 ∈ E, and the function v is Lipschitz continuous on the set E, and locally Lipschitz continuous in R n as needed. Moreover, the function v is not increasing on any trajectory of the system (1), indeed for any

x 0 ∈ R n : v(X(t, x 0 )) = sup τ 0 X(τ, X(t, x 0 )) r = sup τ t X(τ, x 0 ) r sup τ 0 X(τ, x 0 ) r = v(x 0 ).
In addition, v is r-homogeneous of degree 1 since for all x 0 ∈ R n and λ > 0:

v(Λ r (λ)x 0 ) = sup t 0 X(t, Λ r (λ)x 0 ) r = sup t 0 Λ r (λ)X(λ ν t, x 0 ) r = λ sup t 0 X(λ ν t, x 0 ) r = λ sup τ 0 X(τ, x 0 ) r = λv(x 0 ),
where the change of variables τ = λ ν t was used. Now, define a new function for all x 0 ∈ R n :

V (x 0 ) = sup t 0 {v(X(t, x 0 )) k(v(x 0 ) ν t)},
where k : R + → R + is a continuously differentiable function satisfying the conditions given in the formulation of lemma.

The function V has a lower bound κ 1 x 0 r V (x 0 ) and V (0) = 0. In addition, for all x 0 ∈ R n and λ > 0:

V (Λ r (λ)x 0 ) = sup t 0 {v(X(t, Λ r (λ)x 0 )) k(v(Λ r (λ)x 0 ) ν t)} = sup t 0 {v(Λ r (λ)X(λ ν t, x 0 )) k(λ ν v(x 0 ) ν t)} = λ sup t 0 {v(X(λ ν t, x 0 )) k(λ ν v(x 0 ) ν t)} = λ sup τ 0 {v(X(τ, x 0 )) k(v(x 0 ) ν τ )} = λV (x 0 ),
where the substitution τ = λ ν t has been used again, and V is r-homogeneous of degree 1. Next, for any 0 < ι 1 ≤ ι 2 < +∞ and any

x 1 , x 2 ∈ S = {x ∈ R n : ι 1 ≤ v(x) ≤ ι 2 } consider |V (x 1 ) -V (x 2 )| = | sup t 0 {v(X(t, x 1 ))k(v(x 1 ) ν t)} -sup t 0 {v(X(t, x 2 )) k(v(x 2 ) ν t)}| = | sup t 0 {v(X(t, x 1 ))k(v(x 1 ) ν t)} -sup t 0 {v(X(t, x 2 )) k(v(x 1 ) ν t) -sup t 0 {v(X(t, x 2 )) k(v(x 2 ) ν t) + sup t 0 {v(X(t, x 2 )) k(v(x 1 ) ν t)}| ≤ | sup t 0 {v(X(t, x 1 )) k(v(x 1 ) ν t)} -sup t 0 {v(X(t, x 2 )) k(v(x 1 ) ν t)}| +| sup t 0 {v(X(t, x 2 )) k(v(x 1 ) ν t)} -sup t 0 {v(X(t, x 2 )) k(v(x 2 ) ν t)| sup t 0 |k(v(x 1 ) ν t)[v(X(t, x 1 )) -v(X(t, x 2 ))]| + sup t 0 |v(X(t, x 2 ))[k(v(x 1 ) ν t) -k(v(x 2 ) ν t)]| κ 2 sup t 0 |v(X(t, x 1 )) -v(X(t, x 2 ))| + ι 2 sup t 0 |k(v(x 1 ) ν t) -k(v(x 2 ) ν t)|.
For any x 1 , x 2 ∈ S, there is 0 < T < +∞ such that

sup t 0 |v(X(t, x 1 )) -v(X(t, x 2 ))| = sup 0 t T |v(X(t, x 1 )) -v(X(t, x 2 ))|.
This property follows from the fact that v is not increasing and that v(X(+∞, x 1 )) = v(X(+∞, x 2 )) = 0 for an asymptotically stable system (1). Next, using the Mean value theorem

sup t 0 |k(v(x 1 ) ν t) -k(v(x 2 ) ν t)| = sup t 0 | k(θt)t(v(x 1 ) ν -v(x 2 ) ν )|, where θ ∈ [v(x 1 ) ν , v(x 2 ) ν ] ⊂ [ι ν 1 , ι ν 2 ], then sup t 0 |k(v(x 1 ) ν t) -k(v(x 2 ) ν t)| ≤ κ 3 ι -ν 1 |v(x 1 ) ν -v(x 2 ) ν | ≤ κ 4 x 1 -x 2
and the existence of κ 4 > 0 follows from the fact that v is locally Lipschitz, ν ≥ 0 and min{v(x 1 ), v(x 2 )} ≥ ι 1 > 0. Therefore,

|V (x 1 ) -V (x 2 )| κ 2 sup 0 t T |v(X(t, x 1 )) -v(X(t, x 2 ))| + κ 4 x 1 -x 2 κ 2 L S,T X(t, x 1 ) -X(t, x 2 ) + κ 4 x 1 -x 2 (κ 2 L S,T K S,T + κ 4 ) x 1 -x 2
for all x 1 , x 2 ∈ S. Therefore, the function V is locally Lipschitz continuous on S for any 0 < ι 1 ≤ ι 2 < +∞, and by homogeneity it inherits this property for all R n . In addition, V is strictly decreasing for any x 0 ∈ R n \ {0}:

V (X(t, x 0 )) = sup τ 0 {v(X[τ, X(t, x 0 )]) k(v(X(t, x 0 )) ν τ )} = sup τ t {v(X[τ, x 0 ]) k[v(X(t, x 0 )) ν (τ -t)]} < sup τ 0 {v(X[τ, x 0 ]) k(v(x 0 ) ν τ )} = V (x 0 ),
where on the last step the facts, that the function v is not increasing and k is strictly increasing, have been used. Denote

V (t) = V (X(t, x 0 )) for any x 0 ∈ R n and t ∈ R + , then V + (t) = D + V (X(t, x 0 ))f (X(t, x 0 )
) almost everywhere and V + (t) := lim sup

h→0 + h -1 [V (X(h, x 0 )) -V (x 0 )],
then V + (t) < 0 (since it has been established above that V (X(t, x 0 )) < V (x 0 )) and, hence,

D + V (x 0 )f (x 0 ) < 0 for almost all x 0 ∈ R n \ {0}.
The drawback of the above construction of v and V is that for their calculation we have to know the solutions X(t, x 0 ) for all t ≥ 0, which is an obstruction for application. Using homogeneity once more, this shortage can be avoided: Lemma 2. Let Assumption 1 be satisfied. Then

V (x) = sup 0 t T q (x) {v(X(t, x)) k(v(x) ν t)}, v(x) = sup 0 t T q (x) X(t, x) r ,
where

k : R + → R + is a continuously differentiable function satisfying 0 < κ 1 k(t) κ 2 < ( b a ) -1
µ qκ 1 with a monotonically decreasing function k(t) > 0 such that k(t)t ≤ κ 3 < +∞ for all t 0, is a locally Lipschitz continuous and r-homogeneous of degree 1 function such that

κ 1 x r V (x) for all x ∈ R n , D + V (x)f (x) < 0 for a.a. x ∈ R n \ {0}. Proof: For any x 0 ∈ R n define v(x 0 ) = sup 0 t T q (x0) X(t, x 0 ) r ,
where T q (x 0 ) is the upper bound function of contraction in q > 1 times from Definition 4. Note that v(x 0 ) = sup

0 t T q (x0) X(t, x 0 ) r = sup t≥0 X(t, x 0 ) r
and the proof of all properties of the function v is the same as in Lemma 1. Now, define a new function for all x 0 ∈ R n :

V (x 0 ) = sup 0 t T q (x0) {v(X(t, x 0 )) k(v(x 0 ) ν t)},
where k : R + → R + is a continuously differentiable function satisfying the conditions of the lemma. The function V has a lower bound κ 1 x 0 r V (x 0 ) and V (0) = 0. In addition, for all x 0 ∈ R n and λ > 0:

V (Λ r (λ)x 0 ) = sup 0 t T q (Λr(λ)x0) {v(X(t, Λ r (λ)x 0 )) k(v(Λ r (λ)x 0 ) ν t)} = sup 0 t λ -ν T q (x0) {v(Λ r (λ)X(λ ν t, x 0 )) k(λ ν v(x 0 ) ν t)} = λ sup 0 t λ -ν T q (x0) {v(X(λ ν t, x 0 )) k(λ ν v(x 0 ) ν t)} = λ sup 0 τ T q (x0) {v(X(τ, x 0 )) k(v(x 0 ) ν τ )} = λV (x 0 ),
where the substitution τ = λ ν t has been used again, and V is r-homogeneous of degree 1. Next, for any 0 < ι 1 ≤ ι 2 < +∞ and any

x 1 , x 2 ∈ S = {x ∈ R n : ι 1 ≤ v(x) ≤ ι 2 } consider |V (x 1 ) -V (x 2 )| = sup 0 t T q (x1) {v(X(t, x 1 )) k(v(x 1 ) ν t)} - sup 0 t T q (x2) {v(X(t, x 2 )) k(v(x 2 ) ν t)} ≤ sup 0 t T q (x1) {v(X(t, x 1 ))k(v(x 1 ) ν t)} - sup 0 t T q (x2) {v(X(t, x 2 ))k(v(x 1 ) ν t)} + 0 t T q (x2) {v(X(t, x 2 )) k(v(x 1 ) ν t)} - sup 0 t T q (x2) {v(X(t, x 2 )) k(v(x 2 ) ν t)} sup 0 t T |k(v(x 1 ) ν t) [v(X(t, x 1 )) -v(X(t, x 2 ))]| + sup 0 t T q (x2) |v(X(t, x 2 )){k(v(x 1 ) ν t) -k(v(x 2 ) ν t)}| κ 2 sup 0 t T |v(X(t, x 1 )) -v(X(t, x 2 ))| + ι 2 sup 0 t T |k(v(x 1 ) ν t) -k(v(x 2 ) ν t)| ,
where as before T = sup x0∈S T q (x 0 ). Recalling the properties of the function k(t) and using the Mean value theorem we derive: sup

0 t T |k(v(x 1 ) ν t) -k(v(x 2 ) ν t)| = sup 0 t T k(θt)t(v(x 1 ) ν -v(x 2 ) ν ) for θ ⊂ [ι ν 1 , ι ν 2 ], then sup 0 t T k(θt)t(v(x 1 ) ν -v(x 2 ) ν ) ≤ κ 3 ι -ν 1 |v(x 1 ) ν -v(x 2 ) ν |
and since ι 1 > 0 there is κ 4 > 0 such that

κ 3 ι -ν 1 |v(x 1 ) ν -v(x 2 ) ν | ≤ κ 4 x 1 -x 2 .
Consequently, for a locally Lipschitz and stable system (1) we obtain (Remark 2)

|V (x 1 ) -V (x 2 )| κ 2 sup 0 t T |v(X(t, x 1 )) -v(X(t, x 2 ))| + κ 4 x 1 -x 2 κ 2 L S,T X(t, x 1 ) -X(t, x 2 ) + κ 4 x 1 -x 2 (κ 2 L S,T K S,T + κ 4 ) x 1 -x 2
for all x 1 , x 2 ∈ S. Therefore, the function V is locally Lipschitz continuous on S for any 0 < ι 1 ≤ ι 2 < +∞, and by homogeneity it inherits this property for all R n . Recall that v(X(t, x 0 )) is not increasing and k(t) is strictly increasing in time, then V is strictly decreasing for any x 0 ∈ R n \ {0}, and to demonstrate this fact note that

V (X(t, x 0 )) = sup 0 τ T q (X(t,x0)) {v(X[τ, X(t, x 0 )]) k(v(X(t, x 0 )) ν τ )} = sup 0 τ T q (X(t,x0)) {v(X[τ + t, x 0 ]) k(v(X(t, x 0 )) ν τ )} = sup t s t+T q (X(t,x0)) {v(X[s, x 0 ]) k[v(X(t, x 0 )) ν (s -t)]}. If t + T q (X(t, x 0 )) ≤ T q (x 0 ), then obviously V (X(t, x 0 )) ≤ sup t s T q (x0) {v(X[s, x 0 ]) k[v(X(t, x 0 )) ν (s -t)]} < sup 0 s T q (x0) {v(X[s, x 0 ]) k[v(X(t, x 0 )) ν s]} ≤ sup 0 s T q (x0) {v(X[s, x 0 ]) k[v(x 0 ) ν s]} = V (x 0 ).
If t + T q (X(t, x 0 )) > T q (x 0 ), then X(t, x 0 ) r < x 0 r and (if t > T q (x 0 ) then the first supremum above disappears and the equality sign has to be replaced with the sign less or equal)

V (X(t, x 0 )) = max sup t s T q (x0) {v(X[s, x 0 ]) k[v(X(t, x 0 )) ν (s -t)]}, sup T q (x0) s t+T q (X(t,x0)) {v(X[s, x 0 ]) k[v(X(t, x 0 )) ν (s -t)]} ≤ max sup t s T q (x0) {v(X[s, x 0 ]) k[v(X(t, x 0 )) ν (s -t)]}, κ 2 v(X[T q (x 0 ), x 0 ]) ≤ max sup t s T q (x0) {v(X[s, x 0 ]) k[v(X(t, x 0 )) ν (s -t)]}, κ 2 ( b a ) 1 µ q -1 x 0 r . Let κ 2 ( b a ) 1 µ q -1 < κ 1 , then V (X(t, x 0 )) < max sup t s T q (x0) {v(X[s, x 0 ]) k[v(X(t, x 0 )) ν (s -t)]}, V (x 0 )
< max sup

0 s T q (x0) {v(X[s, x 0 ]) k[v(x 0 ) ν s]} , V (x 0 ) = V (x 0 )
and we derived the required property that

V (X(t, x 0 )) < V (x 0 )
for all t > 0. Remark 3. Note that the functions v(x) given in lemmas 1 and 2 are also locally Lipschitz continuous and r-homogeneous of degree 1 such that

x 0 r v(x) ( b a ) 1 µ x 0 r for all x ∈ R n , D + v(x)f (x) 0 for a.a. x ∈ R n \ {0}.
Thus, v(x) are non-strict Lyapunov functions for [START_REF] Bacciotti | Liapunov Functions and Stability in Control Theory[END_REF]. It is worth to note that it is difficult to make a numeric derivation of v since any small computational error may lead to a function v with a positive derivative along trajectories due to its non-strictness.

B. Design based on integral of trajectories

Following [START_REF] Massera | On Liapounoff's conditions of stability[END_REF], [START_REF] Kurzweil | On the reversibility of the first theorem of Lyapunov concerning the stability of motion[END_REF], an alternative integral construction of a strict Lyapunov function can be used: Lemma 3. Let Assumption 1 be satisfied. Then for µ > max{1, ν},

W (x) = +∞ 0 X(t, x) µ r dt
is a locally Lipschitz continuous and r-homogeneous of degree µ -ν function such that

T min q q -µ x µ-ν r W (x), ∀ x ∈ R n , D + W (x)f (x) < 0 for a.a. x ∈ R n \ {0},
where T min q = inf x0∈Sr(1) T q (x 0 ).

Proof: For any x 0 ∈ R n consider

W (x 0 ) = +∞ 0 X(t, x 0 ) µ r dt,
which is a well-defined integral for any x 0 ∈ R n under Assumption 1 for the function β ∈ KL given below (5) (a direct computation confirms this fact). The obtained function W is r-homogeneous of degree µ -ν since for all x 0 ∈ R n and λ > 0:

W (Λ r (λ)x 0 ) = +∞ 0 X(t, Λ r (λ)x 0 ) µ r dt = λ µ +∞ 0 X(λ ν t, x 0 ) µ r dt = λ µ-ν +∞ 0 X(τ, x 0 ) µ r dτ = λ µ-ν W (x 0 ),
where the change of variables τ = λ ν t was used. By construction W (0) = 0 and for any x 0 ∈ S r (1):

W (x 0 ) = +∞ 0 X(t, x 0 ) µ r dt ≥ Tq(x0) 0 q -µ dt = T q (x 0 )q -µ ,
hence by homogeneity T min q q -µ x µ-ν r W (x) for all x ∈ R n . Since µ > 1, X(t, x 0 ) is locally Lipschitz continuous in x 0 for an asymptotically stable system (1) with locally Lipschitz f , and the Lipschitz continuity of • r is established in Corollary 1, then W is locally Lipschitz continuous in R n as needed.

Moreover, the function W is strictly decreasing on any trajectory of the system (1), indeed for any x 0 ∈ R n :

W (X(t, x 0 )) = +∞ 0 X(τ, X(t, x 0 )) µ r dτ = +∞ 0 X(τ + t, x 0 ) µ r dτ = +∞ t X(s, x 0 ) µ r ds < +∞ 0 X(s, x 0 ) µ r ds = W (x 0 ),
where the strict sign is obtained since for any t > 0 and x 0 = 0 there are time instants s ∈ (0, t) such that X(s, x 0 ) r = 0.

For the case of conventional homogeneity, Lemma 3 was proven in [START_REF] Zubov | Methods of A.M. Lyapunov and Their Application[END_REF] (Theorem 36). Finally, a new more practical construction on a fixed time interval might be preferable: Theorem 2. Let Assumption 1 be satisfied. Then there exists q > 1 such that

U (x) = T q ( x r ) 0 X(t, x) µ r dt,
where T q ( x r ) is the upper bound function of contraction in q times (an example is given in (6)), for µ > max{1, ν + ν 2 } is a locally Lipschitz continuous and r-homogeneous of degree µ -ν function such that

T min q q -µ x µ-ν r U (x) T q (1) b a x µ-ν r , ∀ x ∈ R n , D + U (x)f (x) < 0 for a.a. x ∈ R n \ {0},
where T min q = inf x∈Sr(1) T q (x).

Proof: For any

x 0 ∈ R n define U (x 0 ) = T q ( x0 r ) 0 X(t, x 0 ) µ r dt, for U (ξ) = N j=1 θ j p (g(ξ, ξ j )) . (8)
Indeed, the vector θ is the solution of the equation

Υ N = Π N θ, (9)
where

Υ N = [U h 1 . . . U h N ]
T ∈ R N is the vector composed by the corresponding values of U h j and the matrix Π N is nonsingular (symmetric and positive definite) thanks to the properties of the function p.

Assumption 2. The function p (g(ξ, ξ 0 )) is Lipschitz continuous with respect to ξ for any ξ, ξ 0 ∈ S r (1).

This restriction is satisfied for many examples of zonal strictly positive definite functions p (see [START_REF] Cheney | Approximation and Interpolation on Spheres[END_REF], [START_REF] Hubbert | Radial basis function interpolation on the sphere[END_REF] and above). Assumption 3. There exist α > 0 and > 0 such that for any N > 0

Π -1 N 2 ≤ N 1+α
, where • 2 is the induced matrix norm.

Thus, for a sufficiently big value of N (for a sufficiently dense grid on S r (1)) the derivative of ( 8) can be made sufficiently small, which implies that the result of Proposition 2 can be applied and the conditions on ε and ε are satisfied: Theorem 3. Let assumptions 1, 2 and 3 be satisfied. Then there exist q > 1, N > 0 and h > 0 such that for any τ > 0 and µ > max{1, ν + ν 2 }, a locally Lipschitz continuous and r-homogeneous of degree µ -ν Lyapunov function for the system (1) can be found on S r (1) in the form [START_REF] Chetaev | The Stability of Motion[END_REF], where the vector of parameters θ is the solution of (9).

Proof: Since all conditions of Corollary 2 are satisfied, there exist q > 1, N > 0 and h > 0 such that for any τ > 0 and µ > max{1, ν + ν 2 }, a locally Lipschitz continuous and r-homogeneous of degree µ -ν Lyapunov function U : R n → R + for the system (1) can be calculated at the points ξ j for all j = 1, N :

U h j = U (ξ j ).
From the proof of Theorem 2, the Lipschitz constant of U on the sphere S r (1) is upper bounded by the constant L U = T q (1)L Sr(1),T q (1) , which depends only on the properties of the system solutions and it is independent on N . Due to this property, for any N > 0 the functions U obtained in Corollary 2 are uniformly upper bounded on S r (1) by some constant U max ∈ R + (this bound is related with solutions and not with N ), i.e. U h j ≤ U max for any fixed q, h, µ, τ and any N . Since S r (1) is a compact set, then the vector Υ N and the matrix Π N are elementwise bounded uniformly for any N > 0. By Assumption 3, the norm of the inverse of Π N has an upper bound N 1+α for some > 0, then the solution θ of ( 9) is also elementwise upper bounded and decreasing with N :

θ max ≤ θ ≤ Π -1 N 2 Υ N ≤ N Π -1 N 2 Υ N max ≤ N -α Υ N max ≤ N -α U max ,
where • max is the max norm of a vector (elementwise maximum). Define

ε(x) = U (x) - N j=1 θ j p (g(x, ξ j ))
as the approximation error by (8) of the Lyapunov function U . Therefore, due to Assumption 2 there is > 0 such that

L ε = max{L U , N -α U max } is a Lipschitz constant of ε on S r (
|D + N j=1 θ j p (g(y, ξ j )) f (y)| ≤ θ max sup y∈Sr (1)
f (y) ,

and ε ≤ N -α U max sup y∈Sr (1) f (y)
is a strictly decreasing function of N . Therefore, the values -ε and ε are decreasing with N for the approximation (8), and there is a sufficiently big N > 0 such that the conditions of Proposition 2 are satisfied. Note that due to the introduced assumptions and to properties of the function U in (8), the directional derivative D + U (y)f (y) is well defined for all y ∈ S r [START_REF] Bacciotti | Liapunov Functions and Stability in Control Theory[END_REF]. Thanks to the r-homogeneity of f with degree ν and to the r-homogeneity of U of degree µ -ν, D + U (x)f (x) = x µ r D + U (y)f (y) for any x ∈ R n , where y ∈ S r (1) is such that x = Λ r ( x r )y, and D + U (x)f (x) is well-defined for all x ∈ R n .

VII. EXAMPLES

For simplicity of presentation, the case n = 2 is investigated in this section, and the following expressions and values of parameters are used for computation of a homogeneous Lyapunov function: µ = max{1, ν + ν 2 } + 0.15, q = 2, τ = 1.25, h = 0.01, p(s) = |s|.

In the examples below, these values have been fixed and the value of N was increased (similarly h -1 and q may be gradually augmented) until convergence of the procedure. for x ∈ R 2 , then ν = 0.5 for r = [0. 5 1]. The results of simulation are shown in Fig. 2 for N = 100 (trajectories in Fig. 2,a, the values of U in Fig. 2,b and the derivative in Fig. 2,c). The results of these calculations also confirm our theoretical findings.

VIII. CONCLUSION

The problem of construction of a homogeneous Lyapunov function for an asymptotically stable homogeneous system is revisited. First, for systems with nonnegative degree of homogeneity, several expressions of homogeneous Lyapunov functions are proposed, which depend explicitly on the system solutions (on finite or infinite intervals of time). Second, a procedure is presented, which under some technical assumptions ensures that a homogeneous Lyapunov function can be numerically constructed. The results are illustrated by simulations for linear and nonlinear cases. The proposed numeric design of Lyapunov functions can be used for verification of stability of homogeneous systems. Future directions of research include extensions of our approach for the cases of negative degree and discontinuous dynamics. Guidelines for parameter tuning have also to be developed.

  is a locally Lipschitz continuous and r-homogeneous Lyapunov function for[START_REF] Bacciotti | Liapunov Functions and Stability in Control Theory[END_REF].

θ

  [START_REF] Bacciotti | Liapunov Functions and Stability in Control Theory[END_REF]. Since ε(ξ j ) = 0 due to[START_REF] Clarke | Nonsmooth Analysis and Control Theory[END_REF], then the estimate sup y∈Sr(1) |ε(y)| ≤ sup y∈Sr(1) |ε(y) -ε(ξ j )| ≤ sup y∈Sr(1) L ε |y -ξ j | is valid for any j = 1, N . To obtain a tighter upper estimate, for every y ∈ S r (1), we can consider j * (y) ∈ 1, N corresponding to inf |y -ξ j |. So for all y ∈ S r (1), there is κ > 0 such that |y -ξ j * (y) | ≤ κ N since ξ j , j = 1, N form a uniform grid on the sphere S r (1). Therefore, sup y∈Sr(1) |ε(y)| ≤ L ε κ N is decreasing with N . Thus, ε = inf y∈Sr(1) ε(y) ≥ -sup y∈Sr(1) |ε(y)| ≥ -L ε κ N admits a lower bound strictly increasing in N . By definition j p (g(y, ξ j )) f (y)| since D + U (y)f (y) < 0 for all y ∈ S r (1). As it is shown in Proposition 1.1.a of [9]: sup y∈Sr(1)

A. Nonlinear system 1 Consider ( 1 5 -x 1 1 2 3

 11512 ) withf (x) = x 2 -x 1 for x ∈ R 2 , then ν = 2 3 for r = [ 1 3 1].The results of simulation are shown in Fig.1for N = 50. The state trajectories are presented in Fig.1,a. The obtained values of U are plotted in Fig.1,b (the bold points represent U h j ), the derivative of U is given in Fig.1,c. As we can conclude from these results the obtained values correspond to a Lyapunov function, and it is homogeneous by construction.

Figure 1 .Figure 2 .

 12 Figure 1. Results of calculations for nonlinear system 1: a. State trajectories; b. The level of Lyapunov function U (y) for y ∈ Sr(1); c. The values of the derivative D + U (y)f (y) for y ∈ Sr(1)

where T q ( x 0 r ) for q > 1 is from Definition 4, which is defined as a function of x 0 r . The function U is r-homogeneous of degree µ -ν since for all x 0 ∈ R n and λ > 0:

where the change of variables τ = λ ν t was used. By construction U (0) = 0 and for any x 0 ∈ S r (1):

q -µ dt = T q (x 0 )q -µ , and similarly using [START_REF] Briat | Linear Parameter-Varying and Time-Delay Systems[END_REF]:

Hence, by homogeneity

for all x ∈ R n . Further in the proof we will also need the following estimate:

where 1 < q < q (then T q (1) < T q (1)). Denote

If b > a (we can impose such a restriction without loosing generality), then κ > 0 and by homogeneity for all x ∈ R n :

To analyze the continuity of the function U , consider x 1 , x 2 ∈ S r (1) and

Due to Lipschitz continuity of the system (1) there exists L Sr(1),T q (1) ∈ R + (Remark 2) such that

for all 0 t T q (1) and any x 1 , x 2 ∈ S r (1). Consequently,

T q (1)L Sr(1),T q (1) x 1 -x 2 for all x 1 , x 2 ∈ S r (1), and the function U is Lipschitz continuous on the unit sphere S r (1), and by homogeneity it is locally Lipschitz continuous in R n as needed. Moreover, the function U is decreasing on any trajectory of the system (1), indeed for any x 0 ∈ R n :

where the strict sign appears since X(s, x 0 ) = 0 for s ∈ [0, t]. Consider the case < t < T γq ( x 0 r ) where q -1 < γ < 1 is a constant to be determined later (if for all x 0 ∈ R n and t ∈ (0, T γq ( x 0 r )) the function U is strictly decreasing, then it is enough for our purpose), then for t < T γq ( x 0 r ) we have

Hence, in this case again U (X(t, x 0 )) < U (x 0 ). Finally, let us analyze the case

x 0 r has to be considered only, and due to homogeneity of T q , the condition T q ( X(t, x 0 ) r ) > T q ( x 0 r ) implies that the analysis can be even further restricted to the case κ T q (1) b a x 0 r ≤ X(t, x 0 ) r < x 0 r and T q ( x 0 r ) < T q ( X(t, x 0 ) r ) ≤ T q ( x 0 r ), where

. Taking into account all these restrictions we have:

Let us consider x 0 ∈ S r (1), the solution X(t, x 0 ) is a continuously differentiable function for almost all t ≥ 0, then for all t ∈ (0, T γq (1)) and all x 0 ∈ S r (1) there exists L > 0 such that | x 0 ν r -X(t, x 0 ) ν r | ≤ Lt (recall that ν ≥ 0 and the Lipschitz continuity of • r is established in Corollary 1). Consequently, for x 0 ∈ S r (1)

(1)t.

Recall that

, then there exists q > 1 and γ > q -1 such that

Then the property U (X(t, x 0 )) < U (x 0 ) is proven for any x 0 ∈ S r (1), and by homogeneity for any x 0 ∈ R n and t ∈ (0, T γq ( x 0 r )) we have shown that U (X(t, x 0 )) < U (x 0 ).

Remark 4. As we can see in the proof above, the analysis of the properties of U is performed on a compact set separated from the origin (between S r (q -1 ) and S r (1)) and extended globally using the homogeneity arguments. Then the case of ν < 0 or even discontinuous dynamics in ( 1) can be treated similarly if the continuity is lost at the origin only.

VI. NUMERIC DESIGN Let us consider how the Lyapunov function proposed in Theorem 2 can be constructed numerically. By definition, for an r-homogeneous of degree µ -ν Lyapunov function U : R n → R + , for any x ∈ R n there is a unique y ∈ S r (1) such that x = Λ r ( x r )y and

Consequently, it is enough to approximate the values of U on S r (1), and next to reconstruct its values globally via homogeneity arguments. This section has two parts: first, pointwise derivation of values of a homogeneous Lyapunov function is discussed; second, an interpolation of its values on S r (1) is presented under mild restrictions leading to a simple numeric procedure. For a given fixed discretization step h > 0, let X h (t i , x 0 ) denote an approximation at instants t i = ih, i ≥ 0 of the solution X(t, x 0 ) of (1) using Euler method [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] (any other discretization method can also be used).

A. Pointwise calculation

For N > 0, let ξ j ∈ S r (1) with j = 1, N form a uniform grid on the unit sphere S r (1), and let us introduce into consideration the following variables:

where τ > 1 is a tuning parameter.

Corollary 2. Let Assumption 1 be satisfied. Then there exist q > 1, N > 0 and h > 0 such that for any τ > 0 and µ > max{1, ν + ν 2 }: U h j = U (ξ j ) ∀j = 1, N , where U : R n → R + is a locally Lipschitz continuous and r-homogeneous of degree µ -ν Lyapunov function for the system (1).

Proof: Since all conditions of Theorem 2 are satisfied, then for a sufficiently big q > 1 a Lyapunov function can be selected in the form given in the theorem:

By the properties of the (explicit) Euler method [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] (see also Theorem 7 in [START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF] for ν > 0), for N sufficiently big and h sufficiently small, the following inequalities are satisfied for any > 0:

By Proposition 2, there exists a family of r-homogeneous of degree µ-ν Lyapunov functions U (x) that belong to -vicinity of U (x) for some > 0. Since the value of can be made arbitrary small by increasing N and decreasing h, the result follows.

B. Interpolation

The result of Corollary 2 shows how pointwise values of a homogeneous Lyapunov function on S r (1) can be derived. Next, based on this result let us make an interpolation. To this end, note that there is always a homogeneous norm such that S r (1) = S n = {x ∈ R n : x = 1} (for example, an implicit definition of a canonical homogeneous norm from [START_REF] Polyakov | On finite-time stabilization of evolution equations: A homogeneous approach[END_REF] can be used), then for the interpolation we will focus on this case only, and denote the geodesic distance on S n as g(x, y) = arccos(x T y) for any x, y ∈ S n . Following the theory of radial basis function interpolation on the sphere [START_REF] Cheney | Approximation and Interpolation on Spheres[END_REF], [START_REF] Hubbert | Radial basis function interpolation on the sphere[END_REF], a continuous function p : [0, π] → R is (zonal) strictly positive definite on the sphere S n if for all distinct points ξ j ∈ S n with j = 1, N for all N > 0, the matrix

is positive definite, that is λ min (Π N ) > 0, where λ min (Π N ) is the minimum eigenvalue of Π N (some examples of such functions can be found in [START_REF] Cheney | Approximation and Interpolation on Spheres[END_REF], [START_REF] Hubbert | Radial basis function interpolation on the sphere[END_REF], e.g., p(s) = e cos(s) or p(s) = (2 -cos(s)) -1 ). Then, selecting a zonal strictly positive definite function p, there exists always θ = [θ 1 . . . θ N ] T ∈ R N such that U h j = U (ξ j ) ∀j = 1, N


Laurent Mezeix

Christophe Bouvet

Serge Crézé

Dominique Poquillon

Fibers networks as a new type of core material. Processing and mechanical properties

Fibers networks materials have been elaborated from different fibers: metallic fibers, glass fibers and carbon fibers. Cross-links have been achieved using epoxy spraying. The scope of this paper is to analyze the mechanical behavior of these materials and to compare it with available models. The first part of this paper deals with elaboration of fibers network materials. In the second part, compression tests are performed. The specific mechanical behavior obtained is discussed.

INTRODUCTION

Entangled materials can be made from natural materials (wool, cotton…) as well as artificial ones (steel wool, glass wool…). Bonded metal fibers network materials offer advantages [1][2][START_REF] Delannay | Proceedings of MetFoam'99[END_REF][START_REF] Yamada | [END_REF][5][6][7] for use like heat exchanger [8] or insulation [9]. Indeed, they present a low relative density, high porosity and simplicity of production by cost-effective routes with considerable versatility concerning metal composition and network architecture.

On the other hand, sandwich panel consists of two thin skins separated by a thick core. Core material is usually in the form of honeycomb, foam or balsa. Recently, a novel type of sandwich has been developed with bonded metallic fibers as core material [10][11][12][13][START_REF] Dean | Proceedings of ICCS[END_REF][START_REF] Zhou | [END_REF][16]. This material presents an attractive combination of properties like high specific stiffness, good damping capacity and energy absorption. Metal fibers are bonded with a polymeric adhesive [START_REF] Zhou | [END_REF] or fabricated in a mat-like form and consolidated by solid state sintering [12,16].

Entangled cross-linked carbon fibers have been also studied for use as core material by Laurent Mezeix [START_REF] Mezeix | Proceedings of ICCS8[END_REF]. Indeed entangled cross-linked carbon fibers present many advantages for application as core material: open porosity, possibility to mix different fibers to optimize properties or possibility to reeve electric or control cables on core material. Only a few studies are devoted to the mechanical behavior of material made from entangled crosslinked fibers [16,[START_REF] Clyne | [END_REF].

In the present paper, fibers have been chosen in function of their nature, applications and costs. Due to their high performance, carbon fibers are intensely used in many applications: aeronautics, sport equipments or high performance vehicles. Because of their low cost, glass fibers are widely used in many applications like thermal insulation for pipe, building. Bonded metal fiber network are used as core material [10][11][12][13][START_REF] Dean | Proceedings of ICCS[END_REF][START_REF] Zhou | [END_REF][16]. The purpose of the present study is to analyze the mechanical behavior and to model the initial stiffness made from three types of cross-linked fibers.

EXPERIMENT

Material elaboration

Carbon fibers (200 Tex) consist in a yarn of stranded carbon filaments. Fibers were provided by Toho Tenax. The Young modulus of the carbon fiber is 240 GPa. Fibers diameter is 7 µm and the initial epoxy coating represents 1 wt%. Assuming than the coating is uniform, then, its thickness is about 30 nm. Stainless steel fibers were provided by UGITECH. Fibers diameter is 12 µm, the Young modulus is 180 GPa. Glass fibers are obtained from yarns (600 Tex) that were provided by the company PPG Fiber glass. Fibers diameter is 12 µm and their Young modulus is about 73 GPa.

For aeronautical applications, many sandwich-panel skins are made using carbon/epoxy prepreg. That the reason why epoxy resin has been chosen for cross-linking fibers. Epoxy resin was provided by the company SICOMIN. The provided resin has a low viscosity (285 mPa.s) and polymerization duration is 2 hours at 80°C. For all the tests carried out during this work, specimens are carefully weighted using SARTORIUS balance (±10µg). Resin is heated up to 35 °C to decrease viscosity and thus allows a better vaporization. A paint spray gun (Fiac UK air compressors) is used to spray epoxy. Materials were observed using a Scanning Electron Microscope (LEO435VP) operating at 15 kV. Fibers are cut so the fibers length equals 40 mm. For core application, one key parameter is the density. A previous study has shown that the yarns size needs to be decreased by separating filaments to obtain the lowest final density possible [START_REF] Mezeix | Material Science Master's degree[END_REF]. In the present work, separation of yarns was obtained thanks to a blower room. Then, the entangled fibers are cross-linked. In that case, epoxy is sprayed using paint spray gun during the final minutes of the entanglement. We have chosen to perform tests on fiber networks materials with fibers density of 150 kg/m 3 . As the volume of the mould is known and the fibers mass is carefully weighted, the volumetric density is controlled. It is important to notice as the density of bulk carbon in the carbon fibers (1760 kg/m 3 ) is lower than the glass fibers one (2530 kg/m 3 ), which is lower than the stainless steel one (7860 kg/m 3 ), the relative volumetric density of the tests materials differs. Figure 1 shows typical SEM observations on networks made of carbon, glass and stainless steel fibers. We can notice clearly joints between filaments.

Compression tests

The quasi-static compressive response of entangled fibers is measured using a screwdriven test machine MTS with a 5 kN load cell. The sample is introduced between the punches and the compression test is then carried out (no lateral confinement). The sample size is 30x30x30 mm 3 . The punch velocity is v 0 = 1.8 mm.min -1 corresponding to a nominal strain rate of ε& = 10 -3 s -1 . To analyze the experimental results, we used the usual definition for the true strain and true stress.

RESULTS AND DISCUSSION

Cross-linked fibers samples have been tested in compression. Density used for each type of fibers is 150 kg/m 3 . The initial fiber density of the material is determined as explained above. The volume fraction of fiber differs between the three types of materials. 150 kg/m 3 correspond to a volume fraction of 1.9% for the material made of stainless steel fibers, to 5.9% for the material made of glass fibers and to 8.5% for the material made of carbon fibers. Resin is sprayed during entanglement, which is weighted again after epoxy polymerization. So the additional mass of resin can be measured and corresponds to 50 kg/m 3 .

Figure 2 shows comparison during compression tests between materials made with cross-linked carbon fibers, stainless steel fibers and glass fibers. The better stiffness is achieved with the glass fibers. Furthermore, the plateau level [START_REF] Gibson | Cellular solids: structure and properties[END_REF] of stress during compression is quite high. Densification is the dominant mechanism for strain above 70%. The initial stiffness of the networks glass fibers is 8.4 MPa, 1.1 MPa for carbon fibers and 0.6 MPa for stainless steel fibers. The average stress level of the plateau before densification is 0.8 MPa, 0.4 MPa and 0.25 MPa for respectively glass, carbon and stainless steel fibers. In order to get a better understanding of the macroscopic behavior of the cross-linked architecture, important microscopic information is the average distance d av between two cross links. An analysis of SEM pictures has been carried out and results are given in Table 2. As expected after the comparison on the compression behavior, the shortest distance d av is obtained for glass fibers. This is the reason why the initial stiffness of the material is better when compared to the deflection of the beam which depends on the cubic of the distance between two cross links.

This point was not expected as the relative density of the entanglement is higher for the carbon than for the glass. Different surface properties, different reactivity with the epoxy could explain that point.

There are many discussions about the determination of the average distance between fiber contacts. Authors have proposed models to quantify the number of contacts per fiber, in layered structures [START_REF] Batchelor | [END_REF]22] and in three-dimensional network fibers [23,24,25]. The average number of contacts per fiber <c>, in the case of 3D random network, is given by equation (1) as a function of L, the fiber length, of f, the volume fraction of fibers (f) and of D, the fiber diameter. Then, the distance, d av , between joints can be obtained (see table 1) by d av = L/<c>:

f D L c 2 = 〉 〈 (1)
If the initial density of the fibers entanglement is the same for all fibers (150 kg/m 3 ), the volume fraction of the three tested material is not the same. We can notice that the distance between joints obtained by equation ( 1) is closed from SEM observations for material made of glass and material made of stainless steel fibers. In the case of material made of carbon fibers, there is a large difference between the two values. Further investigations on the microscopic organization of the material elaborated in that study are still necessary. X-ray topography measurement would provide more valuable information [16]. T.W. Clyne et al. developed a simple analytical model based on the bending of inclined individual fiber segments [START_REF] Clyne | [END_REF]26,27], the Young's modulus of the network is given as a function E f , the Young's modulus of the fibers and of the geometrical parameters of the network (f, D and d av ) :

2 32 9       = D d f E E av f a (2)
The values obtained with this approach can be compared with the one given by Gibson and Ashby for open cell foams [START_REF] Gibson | Cellular solids: structure and properties[END_REF]. This second model is also based on beam deflections, but with different boundaries conditions (more constrained). The Young's modulus of the network is then given by: Young modulus has been calculated applying equations ( 2) and ( 3) and is given in Table 2. They are compared with the experimental data (initial stiffness of the material). We can notice that the Young modulus obtained using equation ( 3) is closed to the experimental values for carbon and glass fibers. This point could mean that the epoxy drops at the cross links limit the deflection of fibers. In the case of stainless steel fibers, the difference between experiment and equation (3) could be explained by the plastification of fibers and the poor reticulation of the network.

4 4 3       = D d E E av f a π (3)
Tomography data before the compression test and at the beginning of the plateau would be very useful to determine the fibers orientation, the isotropy of the initial entanglement, the average number of contacts per fibers and could help to get a better understanding of the fiber slippage and fiber orientation changes induced by compression.

These observations are now possible and have just been achieved in another study [16]. We plan to get comparable data on our specimens.

CONCLUSIONS

Original materials have been elaborated using fibers networks bonded with epoxy resin. Three different families of fibers have been tested, glass fibers, carbon fibers and stainless steel fibers. The material obtained remains light (200 kg/m 3 ) as the process developed in this study optimize the quantity of epoxy used. The best stiffness is obtained for glass fibers mainly because the shortest distance between cross-links compared to the carbon case. The initial stiffness of the cross-links architecture seems to follow the model proposed by Ashby.

Figure 1 .

 1 SEM observation on materials made with (a) cross-linked carbon fibers, (b) crosslinked glass fibers and (c) cross-linked stainless steel fibers.

Figure 2 .

 2 Figure 2. Comparison of compression stress/strain curves of architectures made with entangled fibers cross-linked with epoxy resin (initial density = 150 kg/m 3 ).

Table 1 :

 1 Comparison of distances between cross-links.

			Average distance	Average distance
	Material	<c>, number of contacts per fiber	between joints : calculated value	between joints observed by SEM
			d av = L/<c> [mm]	[mm]
	Network of stainless steel fibers	127	0.31	0.30
	Network of carbon fibers	974	0.04	0.20
	Network of glass fibers	395	0.10	0.15

Table 2 :

 2 Comparison of Young modulus obtained during compression tests with the values calculated using equations (2) and (3).

	E: Young modulus [MPa]


Stéphane Paquelet

Vincent Savaux

On the Symmetry of FIR Filter with Linear Phase

Keywords: FIR Filter, Group Delay, Linear Phase

This paper deals with signal processing theory related to finite impulse response (FIR) filters with linear phase. The aim is to show a new and complete proof of the equivalence between linear phase and symmetry or antisymmetry of the real coefficients of the filter. Despite numerous proofs are available in the literature, they are usually incomplete, even though the result is commonly used by the signal processing community. We hereby address a pending issue in digital signal processing: we first prove the uniqueness of the group delay for any decomposition amplitude-phase of the frequency response. Based on this first step, we then derive a complete proof of the equivalence: a FIR filter has (anti)symmetric coefficients if and only if the phase is linear. It must be emphasized that this brief paper deals with theoretical aspects of FIR filters.

I. Introduction

Analysis and design of digital finite impulse response (FIR) filters with linear phase are basics in the field of signal processing [START_REF] Proakis | Digital Signal Processing: Principles, Algorithms, and Applications[END_REF]- [START_REF] Saramäki | Handbook for Digital Signal Processing[END_REF]. Such filters are widely used in electronic devices for their inherent stability. Furthermore, numerous techniques for FIR filter design and applications have been proposed, from early 70's [START_REF] Rabiner | FIR Digital Filter Design Techniques Using Weighted Chebyshev Approximation[END_REF] to the present day [START_REF] Pak | Distributed Hybrid Particle/FIR Filtering for Mitigating NLOS Effects in TOA Based Localization Using Wireless Sensor Networks[END_REF], [START_REF] Ahn | Deadbeat Dissipative FIR Filtering[END_REF]. However, unlike papers like [START_REF] Rabiner | FIR Digital Filter Design Techniques Using Weighted Chebyshev Approximation[END_REF], [START_REF] Pak | Distributed Hybrid Particle/FIR Filtering for Mitigating NLOS Effects in TOA Based Localization Using Wireless Sensor Networks[END_REF], [START_REF] Ahn | Deadbeat Dissipative FIR Filtering[END_REF] which deal with design and applications of FIR filters, the present paper is dedicated to theoretical aspects related to FIR filters.

Manuscript submitted in June 2017

Stéphane Paquelet and Vincent Savaux are with IRT b<>com, Rennes, FR * Corresponding author: e-mail: vincent.savaux@b-com.com Phone: +33 256358216.

It is commonly accepted that a FIR filter has a linear phase if and only if the taps are either symmetric or antisymmetric. However, to the best of our knowledge, although this result is commonly accepted, the proof leading to this fundamental property is not satisfactorily demonstrated in the literature [START_REF] Proakis | Digital Signal Processing: Principles, Algorithms, and Applications[END_REF], [START_REF] Bilmes | Lecture 15: Generalized Linear Phases and FIR Types[END_REF], [START_REF] Diniz | Digital Signal Processing, System Analysis and Design[END_REF], in which mathematical shortcuts are used. In particular, showing that a (anti)symmetric filter has a linear phase is straightforward (see [START_REF] Bilmes | Lecture 15: Generalized Linear Phases and FIR Types[END_REF] for instance), whereas the inverse implication requires more difficult developments. Otherwise, the authors of [START_REF] Diniz | Digital Signal Processing, System Analysis and Design[END_REF] made an attempt to prove the double implication, but they inherently presupposed some features regarding the group delay and the differentiability of the frequency response of FIR filters.

In this paper, we revisit the proof of the equivalence between linear phase and (anti)symmetry of FIR filters. To this end, we first prove the existence of the decomposition module/phase of the frequency response of the filter, where both the module and phase functions are at least of class C 1 . Furthermore, we show that the group delay is unique for any decomposition. To the best of our knowledge, both proofs are original contributions of this paper, and they are a mandatory prerequisite for the proof of the equivalence between linear phase and (anti)symmetry, since it justifies the existence of both the phase/module decomposition and the group delay.

Once these properties are justified, we can prove the equivalence between linear phase and (anti)symmetry of FIR filters with proper assumptions. Showing that a FIR filter has a linear phase if it is (anti)symmetric is not tackled in this paper, since it has been undertaken in [START_REF] Proakis | Digital Signal Processing: Principles, Algorithms, and Applications[END_REF], [START_REF] Bilmes | Lecture 15: Generalized Linear Phases and FIR Types[END_REF].

However, we focus on the inverse implication, i.e. we prove that a FIR filter is (anti)symmetric if the phase is linear. This is another original contribution of this paper. Moreover, we consider a very general approach in the sense that the filter coefficients are considered as elements of a mathematical sequence, and no assumption about FIR filters features is made. From this, the proof is based on the analysis of the linear independence of sin functions. To summarize, this paper closes the theoretical concepts related to linearity of the phase of FIR filters, (anti)symmetry of the taps, and uniqueness of the group delay.

The paper is organized as follows: in Section II, we prove the existence of a decomposition module/phase of the frequency response of any FIR filter, and the uniqueness of the group delay.

In Section III, we show the following implication: if a FIR filter has a linear phase, then the coefficients are symmetric or antisymmetric, which completes the existing results concenring the inverse implication. Finally, Section IV concludes this paper.

II. Background and Prerequisites

A. Standard Mathematical Notations

The following notations are used throughout the paper: the class C k (A) of a function f means that the k-th derivative of f exists and is continuous on the set A; the braces {.} define discrete sets, and |f | is the modulus of f . The mathematical equivalence between two assertions (A) ⇔ (B) means that we have the double implication (A) ⇒ (B) and (A) ⇐ (B).

B. Model

We consider a finite impulse response (FIR) filter composed of M real coefficients {h m , 0 ≤ m ≤ M -1}. The filter frequency response is obtained by means of the discrete-time Fourier transform (DTFT) as follows:

H(ν) = M -1 m=0 h m e -2jπmν , (1)
where

ν ∈ I = [-1 2 , 1 2 ]
is the normalized frequency. In order to highlight the amplitude response of the FIR (not to be confused with the magnitude |H| [START_REF] Rabiner | FIR Digital Filter Design Techniques Using Weighted Chebyshev Approximation[END_REF], [START_REF] Saramäki | Handbook for Digital Signal Processing[END_REF]), it is more relevant to consider, if it exists, the following decomposition of H(ν):

H(ν) = H r (ν)e jφ(ν) , (2)
where H r and φ must be (at least) C 1 (R) functions. Note that the differentiability class of H r is rarely discussed in the literature, even it is usually stated that it is sufficient for H r to be C 0 (R).

Now, it must be emphasized that the existence of the group delay, defined as α = -

1 2π dφ(ν) dν ,
presupposes the existence of a decomposition (2) with φ a C 1 (R) function, such that the derivative of φ is unique for any possible decomposition. As a consequence, we deduce that H r is necessarily ν) , where H(ν) is C ∞ by construction [START_REF] Proakis | Digital Signal Processing: Principles, Algorithms, and Applications[END_REF]. Such assertions are proved hereafter.

C 1 (R) since H r (ν) = H(ν)e -jφ(

C. Existence of the Decomposition (2)

Theorem 1. ∀ν ∈ I = [-1 2 , 1 2 ], ∃H r and φ of class C 1 (R), such that H(ν) defined in (1) can be decomposed as H(ν) = H r (ν)e jφ(ν) .
Proof. Since any polynomial can be decomposed as a product of monomials, one can show that every monomial is C 1 (R) to prove that the polynomial is C 1 (R). For the sake of clarity, we rewrite

(1) by using the z transform notation, where z = e 2jπν . Thus, if we suppose that h M -1 = 0, (

is factorized as:

H(z) = M -1 m=0 h m z -m = h M -1 M -1 m=0 (z -1 -z -1 m ), (3)
where z -1 m is the root of the monomial (z

-1 -z -1 m ). Let us denote z m = ρ m e jθm , where (ρ m , θ m ) ∈ R 2 .
We aim to show that, for any m ∈ {0, 1, .., M -1}, the monomials can be expressed as

z -1 -z -1 m = H m (ν) = H m,r (ν)e jφm(ν) , (4)
where H m,r and φ m are C 1 (R) functions, which will be proved afterward.

Two cases have to be dealt with, according to ρ m = 1 or ρ m = 1:

• if ρ m = 1, some straightforward trigonometric manipulations allow us to rewrite the monomial as follows:

z -1 -z -1 m = e -2jπν -e -jθm (5) = 2 sin πν - θ m 2 e -j(πν+ θm 2 + π 2 ) , (6)
where we identify The existence of the decomposition in (2) being proved, we can now prove that the group delay is unique for any decomposition module/phase of H. To this end, we show the uniqueness of the derivative of φ.

H m,r (ν) = 2 sin πν -θm 2 and φ m (ν) = -(πν + θm 2 + π 2 ). Note that both functions are of class C ∞ (R), hence C 1 (R). • If ρ m = 1, H m is of class C ∞ from R to C * ,

D. Uniqueness of the Derivative of φ

Theorem 2. The derivative of φ is unique for any decomposition of H defined in [START_REF] Rabiner | FIR Digital Filter Design Techniques Using Weighted Chebyshev Approximation[END_REF].

Proof. Let ν n and ν n+1 be two consecutive zeros of H. We consider the interval

J n =]ν n , ν n+1 [,
where H(ν) = 0. Furthermore, we denote by (H r,1 , φ 1 ) and (H r,2 , φ 2 ) two possible decompositions of H. Thus, for any ν

∈ J n , since |H r,1 (ν)| = |H r,2 (ν)|, we have H r,1 (ν) = κ(ν)H r,2 (ν)
, where κ(ν) = {±1}. Therefore, we deduce that:

e j(φ1(ν)-φ2(ν)) = ±1, (7)
and then φ1(ν)-φ2(ν) π = k n , where k n ∈ Z. We proved that φ 1 and φ 2 are continuous on R, hence they are on Jn ∋ ν n+1 , thus φ1(νn+1)-φ2(νn+1) π = k n , and on Jn+1 ∋ ν n+1 , thus φ1(νn+1)-φ2(νn+1) π = k n+1 . Then, for any n, we deduce that k n = k n+1 = k, and κ(ν) = (-1) k . To summarize, the above results show that every possible decomposition of H can be deduced from a unique one (H r , φ)

as ((-1) k H r , φ + kπ). Finally, this points out that the derivative of φ (which is proportional to the group delay if the phase is linear) is unique and independent of the decomposition, which concludes the proof.

Such as aforementioned, it is usual to consider H r (ν) = |H(ν)| in the literature, or otherwise the differentiability class of H r and φ is not discussed [START_REF] Rabiner | FIR Digital Filter Design Techniques Using Weighted Chebyshev Approximation[END_REF], [START_REF] Saramäki | Handbook for Digital Signal Processing[END_REF]. However, it can be noticed that if the C 0 (R) condition is imposed to H r instead of C 1 , as in the case of the modulus on the zeros of H, then φ is not C 1 (R) and the group delay is not defined anymore. Furthermore, it must be noticed that all the possible decompositions ((-1) k H r , φ + kπ), k ∈ Z, are of the same class.

Since at least one is of class C ∞ (see proof of Theorem 1 for instance), then all decompositions are C ∞ .

III. Symmetry of FIR Filter with Linear Phase

The aim of this section is to show that any FIR filter has a linear phase if and only if their real coefficients are either symmetric or antisymmetric. Mathematically, this can be expressed as

(A)          H(ν) = H r (ν)e jφ(ν) φ(ν) = β -2παν ⇔                  α = M -1 2 β ∈ {kπ, π 2 + kπ} h m = ε β h M -1-m (B), (8)
ε β =          1, if β = π 2 + kπ -1, if β = kπ , ( 9
)
k ∈ Z, (α, β) ∈ R 2 , ( 10
)
The proof of the implication (A) ⇐ (B) is straightforward and frequently proposed in the literature such as [START_REF] Proakis | Digital Signal Processing: Principles, Algorithms, and Applications[END_REF], [START_REF] Bilmes | Lecture 15: Generalized Linear Phases and FIR Types[END_REF], therefore it is not reported in this paper. However, we propose an original proof of the implication (A) ⇒ (B). This is undertaken by considering a very general context, where no a priori constraint on the values of the filter coefficients h m is imposed. In particular, it is mathematically possible to assume h 0 = 0 or h M -1 = 0, even though, in practice, it is unlikely that the extremal coefficients are zero.

A. Preliminaries

Let J be an interval J ⊂ I such as ∀ν ∈ J, M -1 m=0 h m cos(2πmν) = 0. The existence of J is guaranteed since the number of roots of a trigonometric polynomial is finite. If non-zero FIR filters with linear phase exists, for any ν ∈ J, the decomposition of the equation (1) = (2) into real and imaginary parts yields

         H r (ν) cos(β -2παν) = M -1 m=0 h m cos(2πmν) H r (ν) sin(β -2παν) = -M -1 m=0 h m sin(2πmν) . ( 11
)
Hence, for any ν ∈ J, we can write:

tan(β -2παν) = -M -1 m=0 h m sin(2πmν) M -1 m=0 h m cos(2πmν) . ( 12
)
Then, after some straightforward trigonometric manipulations, we obtain:

M -1 m=0 h m sin(2π (α -m) αm ν + β) = 0, ( 13
)
where α m is defined for the sake of readability. In the following, we analyze the linear independence of the elements of set {sin(2πα m ν + β), 0 ≤ m ≤ M -1, ν ∈ J} in order to deduce the symmetry/antisymmetry properties of the coefficients h m .

It is usually claimed in the literature that no additional result can be deduced from (13), except that "obvious" non-null solutions are obtained when β ∈ {kπ, π 2 + kπ}, k ∈ Z, and α = M -1 2 , and in that case the filter is either symmetric or antisymmetric. However, this is a shortcut that should be mathematically proved. We hereafter provide an original and complete proof of (A) ⇒ (B) from ( 13), from which we deduce the α value based on mathematical considerations.

B. Linear Independence and Notations

It must be emphasized that if the |α m | elements in (13) are all distinct, then the functions appears twice, denoted by E 2 . By construction,

ν ∈ J → sin(2πα m ν + β), 0 ≤ m ≤ M -1,
E 1 ∪ E 2 = {0, 1, .., M -1} and E 1 ∩ E 2 = {∅}:
it is a partition of {0, 1, .., M -1}.

Accordingly, the development of the sum in (13) leads to the following expression: In the latter case, m = α and h α ∈ R. Moreover, note that, if α is an integer, from (14) α ∈ {0, 1, .., M -1}, and α ∈ E 1 because |α -m| = 0 can appear only once.

m∈E1 h m sin(2πα m ν + β) (15) + m,m ′ ∈E2 (h m + h m ′ ) cos(2πα m ν) sin(β) (16)
+ (h m -h m ′ ) sin(2πα m ν) cos(β) = 0. (17)

C. h Symmetry and β Possible Values

h m + h m ′ = h m -h m ′ = 0, i.e.: h m = h m ′ = 0. ( 18
)
This result imposes specific values to β and coefficient pairing in order to obtain non-zero solutions, namely:

• β = 0 [π],

D. α Value

Let M in and M ax be the indexes defined such as follows:

M in = arg min m h m = 0, ( 19
)
M ax = arg max m h m = 0. ( 20
)
The above definition of M in and M ax is illustrated in Fig. 1, where M = 10. In this example, it can be observed that M in = 2 and M ax = 8. We now investigate the belonging of M in and M ax to E 1 and E 2 in order to derive the α value:

• Case 1: M in and M ax belong to different subsets. For instance, suppose that M in ∈ E 1 and M ax ∈ E 2 . By examining E 1 in Section III-C, we know that h M in = 0, then α M in = 0,

which yields α = M in and β = kπ, k ∈ Z. Since M ax ∈ E 2 , E 2 = {∅}, which contains at least two indexes m, m ′ necessarily larger than M in. Hence α = m+m ′ 2 > M in, which is in contradiction with M in = α.
The same reasoning can be carried out when M ax ∈ E 1 , which also leads to a contradiction.

As a consequence, we deduce that both M in and M ax belong to either E 1 or E 2 , such as hereafter investigated.

• Case 2: M in, M ax ∈ 2 : we deduce from (14

) that h M in = h 2α-M in and h M ax = h 2α-M ax ,
where 2α -M in and 2α -M ax belong to E 2 . Hence, we obtain M in ≤ 2α -M ax and 2α -M in ≤ M ax and therefore M in+M ax

2 ≤ α ≤ M in+M ax 2 , i.e.: α = M in + M ax 2 . ( 21
)
Combined with those presented in Section III-C, the result in (21) highlights that the coefficients of the filter are (anti)symmetric to α.

• Case 3: M in, M ax ∈ E 1 : the unique trivial solution is M in = M ax = α, β = kπ, k ∈ Z. This result corresponds to the case where the filter is composed of a unique non-zero coefficient

h m = h α ∈ R.

E. Summary

It has been shown that a constant group delay necessarily implies that:

• α = M in+M ax 2 , • β ∈ {kπ, π 2 + kπ}, k ∈ Z, and ε β =          1, if β = π 2 + kπ -1, if β = kπ .
• h m = ε β h 2α-m , which shows that the coefficients are (anti)symmetric to α.

It must be emphasized that in practice, the extremal coefficients are non-zero values, i.e.

M in = 0 and M ax = M -1. Eventually, even if it means discarding null coefficients at the edges, we recognize the usual expression α = M -1 2 such as defined in (10). Thus, we have shown that (A) ⇒ (B), which finally proves the equivalence (A) ⇔ (B).

IV. Conclusion

In this paper, we proposed an original proof of the equivalence between the linear phase of FIR filters and the (anti)symmetry of the coefficients. As a prerequisite, the uniqueness of the group delay for any decomposition of the filter frequency response has been proved as well. The proof of equivalence has been carried out by considering a very general model, in which the taps of the filters are elements of a mathematical sequence. From this paradigm, we analyzed the independence of sin functions in order to deduce some properties of the FIR filters, namely the (anti)symmetry property of the coefficients, as well as the values of the linear phase parameters α and β. It is worth noting that the proposed proof remains valid if FIR filters with complex coefficients are considered. In that case, however, the coefficients are conjugate (anti)symmetric.

  as H m (ν) = 0 for any ν ∈ I. Then, from lifting property of the map p : R → S 1 given by p(ν) = Hm(ν) |Hm(ν)| , where S 1 is the unit circle (see[START_REF] Hatcher | Algebraic Topology[END_REF] for more details), we know that the polar form of H m (ν) = H m,r (ν)e jφm(ν) exists, withH m,r = |H m |and φ m of the same class as H m . Hence H m,r and φ m are of class C 1 (R). This concludes the proof.

  are linearly independent, except if one of the sin function is null. In any case, the only solution is either the null filter, or a trivial single-coefficient filter. As a consequence, we deduce that a non-trivial solution exists only if there are at least two different indexes m and m ′ such that |α m | = |α m ′ |, or equivalently α m = -α m ′ , and hence: m + m ′ 2 = α. (14) This remark suggests splitting the index of the sum in (13) into two subsets: one gathering of indexes m for which a value |α m | appears only once, denoted by E 1 , and one for which |α m |

E 1

 1 Examination: By linear independence of the functions sin(2πα m ν + β), we deduce that if m ∈ E 1 , h m = 0, except in the case of a null sin function, i.e. α m = 0 and β = kπ, k ∈ Z.

Fig. 1 .

 1 Fig. 1. Filter coefficients, M in = 2 and M ax = 8.

E 2

 2 Examination: The linear independence of the families sin(2πα m ν) and cos(2πα m ν), in the second sum in (17), yields: if cos(β) sin(β) = 1 2 sin(2β) = 0, i.e. β = k π 2 , k ∈ Z, the expression (17) holds only if

  and in this case (17) holds if h m = h 2α-m , m ∈ E 2 , which corresponds to the symmetric case;• β = π 2 + kπ, k ∈ Z, and in this case (17) holds if h m = -h 2α-m , m ∈ E 2 , which correspondsto the antisymmetric case.

July 3, 2018 DRAFT


D Guillon

S Rivallant
email: samuel.rivallant@isae.fr

J J Barrau

C Petiot

P Thevenet

B Malherbe

EXPERIMENTAL AND NUMERICAL STUDY OF THE SPLAYING MODE CRUSH OF CFRP LAMINATES

Keywords: laminate, crushing, energy absorption, experimental/numerical correlation

Thanks to an innovative plate crushing test fixture, elementary crushing modes of CFRP laminates have been observed and analysed for static and dynamic loading. These analyses enabled to propose a 2D explicit finite element model at mesoscale for the simulation of splaying mode, with a good correlation between model and test data.

INTRODUCTION

Thanks to their high strength-to-weight ratio, composite materials are widely used in the field of transport and especially in aerospace applications, where weight savings are important. But specificity and complexity of composite fractures modes make difficult the prediction of the mechanical behaviour of composite structures subjected to crush loading.

Numerous studies [START_REF] Hull | A unified approach to progressive crushing of fibre-reinforced composite tubes[END_REF][START_REF] Mamalis | Crashworthy capability of composite material structures[END_REF] have been done in the last decades and the interest in crashworthiness is still present as show the recent works on the subject [START_REF] Feraboli | Efforts in the standardization of composite materials crashworthiness energy absorption[END_REF]. But there is still a lake in numerical simulation. Most of the models developed in the last few years [START_REF] Mcgregor | Simulation of progressive damage development in braided composite tubes under axial compression[END_REF][START_REF] Bisagni | Progressive crushing of fiberreinforced composite structutural components of a formula One racing car[END_REF] are based on global tests characterisation that make the model depend on global parameters, which do not permit to have predictive models. Some models [START_REF] Fleming | Finite element simulation of delamination with application to crashworthy design[END_REF][START_REF] Pinho | Numerical simulation of the crushing process of composite materials[END_REF] are based on material characteristics, but often need an a priori knowledge of the crush damage mode developed in the crush front. The challenge today in crashworthiness simulation is then to be able, from elementary material characterisation data, to predict both crush damage modes and energy absorption in any structures.

In this study, a model for splaying mode of crushing plates is presented. Simulation is based only on elementary material mechanical characteristics. Experimental tests made with an innovative plate crush test fixture enable to have precise data (load, displacement, global behaviour) for experiment/simulation correlation.

EXPERIMENTAL ACHIEVEMENT

This paragraph presents the test fixture and a summary of main results of the experimental study made on composite laminated plates subjected to crush loading. More experimental results can be found in [START_REF] Guillon | Initiation and propagation mechanisms of progressive crushing in carbon-epoxy laminated plates[END_REF].

Test fixture

The test fixture (figure 1) has been designed to improve plate crushing fixtures found in literature. It is made of four vertical adjustable guides which avoid buckling of the specimen, and two horizontal adjustable guides. These ones are localised just beneath the vertical guides so that a gap exists between horizontal guides and the metallic base plate where the laminate crashes. Such a concept has been developed at the same time by Feraboli [3]. This avoid tearing of the plate observed in most of the plate crushing fixtures [START_REF] Lavoie | Dynamic crush tests of energy absorbing laminated composites plates[END_REF][START_REF] Daniel | The crush behaviour of carbon fibre angle-ply reinforcement and the effect of interlaminar shear strength on energy absorption capability[END_REF], and ensure that the boundary conditions are the same along the whole width of the plate, just above the crushing front. Two main uprights carry the guides without interference with the crush front. Visualisation holes have been made to let the crush front visible. The laminate is fixed to a steel cylinder which is the interface with the static machine or the drop tower (dynamic tests). During crushing, plate and cylinder go down between the uprights and the vertical guides. At the bottom of the cylinder, a 120 kN Kistler piezoelectric sensor measures the crushing load. The small distance between the specimen and the sensor limits the mechanical filtering of the signal, allowing high precision in the dynamic crush load measurement. A high speed camera is used to have a real-time observation of the crushing front on the side of the plate. The filming speed is 20000 fps, with a 512x256 pixels resolution.

Specimen

Specimens are 160x60mm flat laminated plates. White graduations are drawn on the edge of the specimen each 5 mm. Each test is defined by three parameters:

-Trigger mechanism: three different triggers, shown in figure 2

Experimental results

The numerous experiments made enabled to observe the different crush modes of composite laminated plates (figure 4). Crush mechanisms are homogeneous through the width of the plate, thus the real-time pictures of the side of the crushed laminate allow observation and understanding of crush modes. Following is shortly described some of the main modes observed.

Undamaged splaying occurs when the laminate is flexible enough to bend without localized damages in plies during the crush. Small ratio between interface toughness and ply strength favour this crush mode. Energy absorbed is low: less than 10 kJ/kg, (see figure 3). Fractured splaying occurs when the laminate suffers repetitive localized damages in plies due to bending during the crush. An obstacle, like a debris wedge, in the front separating two splaying arms is necessary to maintain this mode. Fragmentation occurs when plies reach the metallic base at right angle and can not slip on the base. Failure appears continuously at microscopic scale. Energy absorbed is much higher: up to 50 kJ/kg, depending on the number of damaged plies (see figure 3). But fragmentation is not a stable crush mode, plies crush evolving from fragmentation to bending, leading to undamaged or fractured splaying crush mode.

NUMERICAL SIMULATION

The numerical simulation is based on a dynamic test. The specimen is a laminated plate made of twenty 0,25 mm thick UD T700/M21 plies, with the stacking sequence [0°,45°,-45°,0°,90°] 2/s . The trigger is a metallic blade. The crush mode obtained for this test configuration is a splaying mode with multi-delamination of the plies. The initial speed of the falling mass is 5 m/s. Load/time curve, load/displacement curve and high speed videos of the tests are available to make precise correlation between test and simulation. The calculation is made with Abaqus explicit code.

Model description

In this model, x is the crushing direction and z the direction perpendicular to the plate (see figure 5). Directions 1, 2 and 3 represent fibre direction, transverse direction and through the thickness direction of the ply. 0° corresponds to the crushing direction (x).

The laminated plate is modelled in three-dimension (figure 5). Only 100 mm length of the plate is modelled, to reduce calculation time, which is long enough to have realistic boundary conditions on the metallic plate and the guides. Each layer is modelled with one square thick shell elements in the ply thickness (continuum shell in Abaqus). Dimensions of elements are 0.23 mm in the thickness, and 0,25 mm in the two other directions. The shell element mid-plan corresponds to the ply mid-plan (plan xy). Layers are separated by 3D cohesive elements. The thickness of the interface is 0.02 mm, close to the tenth of the ply thickness. Details of plate mesh and plate stacking sequence are given in fig 6 .   To reduce calculation time, only one element is used in the width of the plate (y direction). The model is then made of 8000 thick shell elements and 7600 cohesive elements, which give around 100000 degrees of freedom.

The metallic blade is defined as a rigid boundary condition with contact (penalty method) and friction (friction coefficient of 0.1). The guides are represented by two rigid half-planes, with similar contact and friction rules.

Additional masses are added at the top of the plate, to represent the mass of the falling weight used in the test. An initial speed of 5 m/s is applied to the plate and these additional masses.

Damage in ply

The pseudo 2D characteristic of the model leads to define different materials for 0° and 90° plies in one hand, and 45° plies in the other hand.

For 0° and 90° plies, the behaviour is linear elastic damageable. Damage initiation is given by Hashin criteria, reduced to plane stresses. The mechanical characteristics are given in table 1, where X T , X C , Y T , Y C and S C are the tensile and compressive strength in fibre direction, the tensile and compressive strength in the transverse direction, and the in-plane shear strength. Then, for each of the four damage mode, a damage variable is defined: dft, dfc, dmt, dmc for fibre rupture mode in tension and compression, and matrix rupture mode in tension and compression. The element rigidity decrease linearly with the strain, which give the following evolution of damage variables:

( ) ( ) 0 0 i f i i i i f i i d ε ε ε ε ε ε - ⋅ - ⋅ = (1)
Where ε i is the strain concerned by the damage mode i, ε i 0 is the strain when Hashin criterion is reached for the damage mode i, and ε i f is the strain when total rupture is reached in mode i.

The rupture strain ε i f is calculated so that the surface energy dissipated in the element is equal to the energy dissipated in an elementary stable compression/tensile crack test propagation (fibre in compression: G FC , fibre in tension: G FT ). The values used in fibre direction are extrapolated from tests made by Pinho [START_REF] Pinho | Fracture toughness of the tensile and compressive fibre failure modes[END_REF], on T300/913. The values for matrix rupture are extrapolated from interlaminar fracture mechanics test. See table 1 for the values.

The material degradation law, taking into account a few couplings is given below:

[ ] ( ) ( )( ) ( )( ) ( ) ( ) [ ] ε α ν ν α σ ⋅ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - - - - - - - = 12 2 2 12 1 21 1 1 0 0 0 1 1 1 0 1 1 1 1 G d E d E d d E d d E d s m m f m f f (2)
Where d f and d m represent tension or compression mode, depending on σ 11 sign,

( )( )( )( ) c m t m c f t f s d d d d d - - - - - = 1 1 1 1 1 ( )( ) 12 21 1 1 1 ν ν α m f d d - - - =
To represent the pseudo plastic behaviour of 45° plies, an elastic-plastic bilinear law is used for shear stress. Damage law is the same as for 0° and 90° plies, applied when Hashin criterion is reached.

Mechanical characteristics of the ply given in table 1 (except G FC and G FT ) come from characterisation tests.

Delamination law

Delamination is represented thanks to cohesive elements, with the following behaviour law:

⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⋅ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ t s n tt ss nn t s n K K K T D t t t δ δ δ 0 0 0 0 0 0 1 0 (3)
Where (t n , t s , t t ) are stresses in the direction 3, 1 and 2, [K] is the rigidity matrix, (δ n , δ s , δ t ) are separations in the direction 3, 1 and 2, T 0 is the initial thickness of the element, and D is the damage variable.

Damage initiates when one of the three stresses reaches its maximal value. Maximal values are the rupture values calculated from characterisation tests of composite ply: t n 0 = σ 33 = 75 MPa in tension, and t s 0 = t t 0 = σ 13 = σ 23 = 150 MPa in shear.

Propagation of damage is ruled by a bilinear softening law. The linear decrease is calculated so that energy absorbed is equal to the energy release rate at interface in each mode. Coupling between the three modes is chosen linear:

1 = ∑ ic i G G (4)
After total rupture of the interfaces, contact with friction is defined between adjacent plies.

Experimental / numerical comparison

Simulation gives results close to experimental data. Figure 7 shows that curves from experiment and simulation are very similar, which means that the main phenomena are well represented: a first little peak load when the plate reaches the metallic blade, followed by a light increase of load when the two separated parts of the plate slide on the blade, then the high increase of load when they reaches the metallic base plate, and finally, the decrease of load and plateau when splaying is stabilised.

The maximum load reached during the peak (50 kN) is inferior in the simulation (35 kN), but due to the experimental dynamics effects of both the structure and the sensor (which are not completely modelled) and to the simplified rules of delamination initiation, it is to compare test and simulation. The weak value of the plateau is coherent with a crush in splaying mode, where there is low damage in plies.

Vibrations observed on the test curve, after the two peaks, are not simulated because the whole structure (steel cylinder, sensor, falling mass) is not modelled, and therefore not excited by the peak load. These vibrations do not correspond to a damage mode of the plate.

Pictures from high speed videos and simulation, in figure 8 also show a good correlation in terms of global displacement and global behaviour.

Analysis of energy dissipation

The analyse of energy absorption in the crushing plateau, from simulation results, shows that in the case of splaying crush mode, in the plateau, most of the energy is dissipated in interfaces delamination. In this test configuration, it represents 60%, and only 15% for friction, and 15% for damage in plies.

Study of influence

Different parameters were changed to make a study of influence.

Change in friction coefficient between composite and metallic parts of the setup, from 0 to 0.3, shows that increase of friction leads to a change in front geometry (decrease of the curvature radius) and an increase of crush load and crush energy absorption. There is very small effect of composite/composite friction coefficient, varying from 0 to 0.4, as the contact surface and relative movement between plies are small during the plateau. The influence of interface parameters is more complex. The increase of interface toughness by the mean of initiation stresses and critical energy release rate leads to different changes of rupture mode: decrease of the number of delaminated interface (other splaying mode, with smaller energy absorption at this step), and then increase of stresses in the plies, which lead to ruptures of plies (towards fractured splaying or fragmentation mode, with bigger energy absorption). This study of influence shows the important sensitivity of the structure behaviour under crush load to the material characteristics.

CONCLUSION

The aim of this study was to propose a model for splaying mode in plates subjected to crush load, only dependant on the elementary material characteristics of the laminate. Thanks to precise experimental measures and observations, it is shown that simulation of crush splaying mode is possible with reasonably correct results, and that this model also enables to predict the influence of material mechanical characteristics on the behaviour of laminated plates under crush load.

However, it is clear that splaying mode is not the crush mode that absorbs the maximum of energy, and even if experiments show that in most of crush mode, outer plies are in splaying mode and help to stabilise the crush front, it is necessary to focus on simulation of fragmentation, which is the more efficient mode in terms of energy absorption.

Figure 1 :

 1 Figure 1: Test fixture: photo and simplified drawing (without uprights)

  , -Crush speed: 20 mm/min for static test, and 5 m/s impact speed for drop tower tests, -Laminate configuration: three materials were tested: 1) Vicotex 914/42%/G803 carbon-epoxy prepreg fabrics with stacking sequences [(0°/90°),(+45°/-45°)] 4/s and [(0°/90°)] 8/s , 2) Hexply M21/35/268/T700GC thick unidirectional carbon-epoxy prepreg with stacking sequences [0°,45°,-45°,0°,90°] 2/s and [45°,0°,-45°,90°] 4/s , 3) Hexply M21/35/134/T700GC thin unidirectional carbon-epoxy prepreg with stacking sequences [45°,-45°] 8/s and [0°,90°] 8/s.

Figure 2 :Figure 3 :

 23 Figure 2: Trigger: a) chamferb) steeple -c)metallic blade

Figure 4 :

 4 Figure 4: Examples of crush modes observed (from left to right) : undamaged splaying mode: laminate of G803/914 -Chamfer trigger Mixed mode: laminate of T700/M21 -chamfer trigger Fragmentation: laminate of G803/914-chamfer trigger

Figure 5 :Figure 6 :

 56 Figure 5: Finite element model

Figure 7 :

 7 Figure 7: Comparison between experiment and simulation: Load-Displacement curve

Figure 8 :

 8 Figure 8: Comparison between experiment and simulation: visualisation of splaying

Table 1 :

 1 Mechanical characteristics of T700/M21 for simulation

	E 1 (GPa)	E 2 (GPa)	G 12 (GPa)	υ 12
	125	9	5	0,4
	X T (MPa)	X C (MPa)	Y T (MPa)	Y C (MPa)	S c (MPa)
	1950	1950	75	220	150
	G FT (mJ/mm²) G FC (mJ/mm²) G I (mJ/mm²)	G II (mJ/mm²)
	133	80	0,35	1,2