L Charroyer

O Chiello

J-J Sinou

Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction

published or not. The documents may come

Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction

Lucien Charroyer, Olivier Chiello,

Introduction

Even if the problem of friction-induced vibrations has been a topic of great interest for many researchers in the past years [START_REF] Crolla | Brake noise and vibration -state of art[END_REF][START_REF] Kindkaid | Automative disc brake squeal[END_REF]3,[START_REF] Ibrahim | Friction-induced vibration, chatter, squeal and chaos part ii : Dynamics and modelling[END_REF][START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF]6], squeal is still a current issue and a difficult subject since the effects of the system parameters are not completely understood. The comprehension of the phenomena involved in these vibrations and the choice of the contact models at the frictional interface have been much investigated [START_REF] Kindkaid | Automative disc brake squeal[END_REF]3,[START_REF] Ibrahim | Friction-induced vibration, chatter, squeal and chaos part ii : Dynamics and modelling[END_REF][START_REF] Ibrahim | Friction-induced vibration, chatter, squeal and chaos part 1: mechanics of contact and friction[END_REF] and are still an active field of research. For example, Massi et al. [8] highlight the significant impact of both the contact surface topography and the frictional contact model. Many mechanisms have been proposed to explain the emergence of friction-induced instabilities. It is generally assumed that four categories of friction-induced vibration may be distinguished [START_REF] Ibrahim | Friction-induced vibration, chatter, squeal and chaos part ii : Dynamics and modelling[END_REF][START_REF] Ibrahim | Friction-induced vibration, chatter, squeal and chaos part 1: mechanics of contact and friction[END_REF]: stick-slip, variable dynamic friction coefficient, sprag-slip and mode coupling. Stick-slip was the first developed theory. Similarly, a decrease of the friction coefficient with relative sliding velocity (which introduces a negative damping due to tribological properties) can be chosen as the source of instability. In this case, a single structural mode can become unstable. Next, academic studies demonstrated that self-excited vibration can occur with a constant friction coefficient. Spurr [START_REF] Spurr | A theory of brake squeal[END_REF] introduced the well-known sprag-slip phenomenon based on a geometrically induced instability. Then, several authors proposed the mode coupling mechanism as the origin of instability. In this latter case, the instability is the result of both the coupling of at least two structural modes of the mechanical system and the coupling between the normal and the tangential dynamics. It is associated to a Hopf bifurcation. This approach is now widely used to reproduce friction-induced vibrations in mechanical systems. Minimal models highlighting this type of instabilities can be easily found in the literature [START_REF] Hoffmann | A minimal models for studying properties of the mode-coupling type instability in friction induced oscillations[END_REF][START_REF] Von Wagner | Minimal models for disk brake squeal[END_REF][START_REF] Sinou | Mode coupling instability in friction-induced vibrations and its dependency on system parameters including damping[END_REF]. In most of these studies, the contact between the mass and the friction plane is modelled by a classical contact spring. It is a relatively common choice not only for this type of minimal models but also for most elaborated finite element models. This choice can be explained by physical considerations (i.e. local elasticity of the contact) or by numerical considerations (i.e. the friction force is directly linked to the normal displacement in order to simplify the problem formulation and the equations of motion). In addition, the authors often restrained the contact formulation at the frictional interface to one tangential degree-of-freedom (i.e. rectilinear friction). However, in this case, some specific phenomena due to the planar nature of the friction Coulomb law cannot be explained, particularly the fact that the direction of the friction force varies with the structure vibration.

A simple model combining these two characteristics (i.e. no contact spring and planar friction) was proposed by Moirot [START_REF] Moirot | Étude de la stabilité d'un équilibre en présence de frottement de Coulomb[END_REF][START_REF] Moirot | Friction and Instabilities, chapter Some examples of friction-induced vibrations and instabilities[END_REF]. He outlined that differences between planar and rectilinear friction formulations can be highlighted and that destabilisation due to the planar friction may occur. However, Moirot's study focused only on a system without damping. On the contrary, similar studies carried out on more complex finite element models [START_REF] Lorang | Instabilité vibratoire des structures en contact frottant : Application au crissement des freins de TGV[END_REF][START_REF] Lorang | Tgv disc brake squeal[END_REF][START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF] showed that planar friction seems to stabilise the system.

Other planar effects of friction have already been investigated in the general context of friction-induced vibrations or brake squeal. Kinkaid et al. [START_REF] Kinkaid | On the transient dynamics of a multi-degree-of-freedom friction oscillator: a new mechanism for disc brake noise[END_REF] proposed a new mechanism for disc brake noise by taking into account both longitudinal and transverse braking directions for a 4-dof model. They demonstrated that instabilities in the radial vibration can be observed during the transient processes. Oberst and Lai [START_REF] Oberst | Nonlinear friction coupling in disc brake squeal. 18th international Congress on Sound and Vibration[END_REF] also studied nonlinear friction coupling in disc brake squeal and more specifically the influence of the lateral vibration of a 2-dof friction oscillator. Considering the variation of the belt angle, they showed that a perturbation in the radial component in a brake system can also cause instabilities. Zhang et al. [START_REF] Zhang | Instability analysis of friction oscillators with uncertainty in the friction law distribution[END_REF] recently studied the stability of coupled friction oscillators on sliding rigid plate with planar friction. They investigated the influence of the uncertainty arising from the tribological aspect and proposed different types of friction modelling with randomised parameters.

However, the precise role of planar friction on mode-coupling instabilities, compared with rectilinear friction, remains an open question. The objective of the present study is to investigate this scientific issue thoroughly. One of the most interesting phenomena concerns the existence of a non-intuitive effect of damping distribution that can generate unstable vibrations. As previously reported in many recent works, structural damping is of primary importance to mode coupling stability and results indicate that the addition of damping in mechanical systems alters the stable-unstable boundaries. In the literature, the relationship between damping and system propensity to develop instability and the so-called "destabilisation paradox" applied to friction-induced vibration have often been illustrated by several analytical and numerical studies [START_REF] Sinou | Mode coupling instability in friction-induced vibrations and its dependency on system parameters including damping[END_REF][START_REF] Hoffmann | Effects of damping on mode-coupling instability in friction induced oscillations[END_REF][START_REF] Sinou | The role of damping and definition of the robust damping factor (rd-factor) for a self-exciting mechanism with constant friction[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode coupling patterns in case of brake squeal[END_REF]. Many contributions on the destabilisation paradox in conservative systems with damping have also been carried out in a general context by Kirillov [24,[START_REF] Kirillov | Destabilization paradox due to breaking the hamiltonian and reversible symmetry[END_REF][START_REF] Kirillov | The effect of small internal and external damping on the stability of distributed non-conservative systems[END_REF][START_REF] Kirillov | Stabilization and destabilization of a circulatory system by small velocitydependent forces[END_REF]. He showed that the effect of proportional and nonproportional damping on the reversible Hopf bifurcation in systems with rectilinear friction is in qualitative agreement with the general theory of the destabilisation paradox in circulatory systems [START_REF] Kirillov | Destabilization paradox due to breaking the hamiltonian and reversible symmetry[END_REF][START_REF] Kirillov | The effect of small internal and external damping on the stability of distributed non-conservative systems[END_REF][START_REF] Kirillov | Stabilization and destabilization of a circulatory system by small velocitydependent forces[END_REF]. Massi and Giannini [28] also proposed an experimental investigation of the relationship between the distribution of modal damping and the propensity to develop squeal in a beam-on-disk setup. They highlighted that a nonuniform repartition of the modal damping causes an increase of the squeal propensity. Extensive studies have also been conducted including the role of damping with other physical factors. For example Hervé et al. [START_REF] Hervé | Extension of the destabilization paradox to limit cycle amplitudes for a nonlinear self-excited system subject to gyroscopic and circulatory actions[END_REF][START_REF] Hervé | Analysis of squeal noise and mode coupling instabilities including damping and gyroscopic effects[END_REF] used a two-degree-of-freedom nonlinear model of clutch squeal in order to examine in detail the influence of structural damping on the effects of the combined circulatory and gyroscopic actions. Recently, Kirillov et al. [START_REF] Kirillov | Destabilization paradox due to breaking the hamiltonian and reversible symmetry[END_REF][START_REF] Kirillov | The effect of small internal and external damping on the stability of distributed non-conservative systems[END_REF][START_REF] Kirillov | Stabilization and destabilization of a circulatory system by small velocitydependent forces[END_REF] worked on the determination in an explicit form of the stabilising damping configurations for a large class of non-conservative systems. They proposed a theory for the qualitative and quantitative description of the destabilisation paradox in such systems. In the case of a vibrational system, planar friction introduces a new contribution of damping which has not been examined in this perspective. The additional inclusion of structural damping also remains an open question. Thus, an original contribution of the present study is to extend well-known results on the non-intuitive effects of damping distribution for the specific case of a damped mechanical system with planar friction.

First, the mechanical system under study, the background on stability analysis and the Routh-Hurwitz criterion are presented. Second, parametric studies and numerical results for damped and undamped systems with planar or rectilinear friction are investigated in order to discuss the elementary effect of a planar or rectilinear friction and to undertake the extension of the destabilisation paradox in the presence of planar friction.

Description of the mechanical model

The system considered in this study is described in Figure 1. It is composed of a mass m in frictional contact with a rigid plane moving with a constant rectilinear velocity. The value of the imposed velocity is denoted V while its direction is denoted t as indicated in Figure 1. The mass is held against the moving plane by three springs and pressed by an external force F . Damping is also included as shown in Figure 1. The angle between the direction t of the imposed velocity of the rigid plane and the coordinate x of the mechanical system is defined by θ. The equations of motion can be written in matrix form:

M Ü + C U + KU = R + F (1)
where

U = U x U y U z T
, U and Ü denote respectively the displacement, velocity and acceleration vectors along the three directions x, y and z. The mass matrix M is defined by M = mI (with I the identity matrix). The vector F represents the external force. The vector R includes the contact normal reaction and the friction force applied on the centre of mass due to frictional contact with the plane. C and K are respectively the damping and stiffness matrices. They are given by :

C =      c x 0 0 0 c y 0 0 0 c z      (2)
K =      k 1 0 k 1 tan α 0 k 2 k 2 tan β k 1 tan α k 2 tan β k n      (3)
F = 0 0 -F n T with F n > 0 (4)
with

k 1 = k 1 cos 2 α (5) k 2 = k 2 cos 2 β (6)
k n = k n + k 1 tan 2 α + k 2 tan 2 β (7)
The coupling between the normal and the tangential degrees-of-freedom is induced by angles α and β. It is a necessary condition for the occurrence of mode coupling instability. The choice of orthogonal springs k 1 and k 2 is only introduced to simplify the formulation of the dynamic equation of the system and to make the comprehension of various phenomena on the instability occurrence easier. It provides uncoupling equations for the tangential degrees-of-freedom (i.e. the two modes concerned by the mode coupling phenomena) without friction. This choice does not limit the results generality, since other configurations can be found with an appropriate modification of the coordinate system. For the readers' comprehension, there is no geometrical symmetry in the (x, y) plane if stiffness k 1 and k 2 or angles α and β are different.

In this study, we are interested in the occurrence of instability around the sliding equilibrium (i.e. the appearances of self-excited system oscillations due to mode coupling) . The contact is assumed to be bilateral and sliding: the mass remains in contact with the plane and the mass tangential velocity remains lower than the velocity of the rigid plane so that no adhesion can occur. The friction force is given by the sliding part of the Coulomb's law with a constant coefficient µ. The contact equations can be expressed as :

                           R = P t T R t + P b T R b + P n T R n U z = 0    R t R b    = µR n        V -Ut -Ub               V -Ut -Ub        (8)
where

P b = -sin θ cos θ 0 (9)
P t = cos θ sin θ 0 (10)

P n = 0 0 1 (11)
R t , R b and R n are respectively the two components of the friction force and the contact reaction which is applied to the mass, according to the directions t, b and z. Ut = P t U and Ub = P b U are the mass velocities according to the directions t and b. Matrices P t , P b and P n are used to transpose forces from the contact to the global coordinate system. For rectilinear friction, the friction force component in direction b is zero and Coulomb's law in the sliding phase is simply written R t = µR n , which is directly linear. For planar friction, equation ( 8) remains non-linear and has to be linearised in order to achieve the stability analysis.

[Figure 1 about here.]

Stability analysis

In this section, the stability of the system is studied. It consists in analysing the behaviour of small perturbations around the quasi-static sliding equilibrium. In this purpose, the calculation of this equilibrium is firstly explained. Then, the friction equations are linearised around the equilibrium and the problem is projected on the tangential plane. Finally the stability of the equilibrium is determined by the classical Complex Eigenvalues Analysis (CEA) or by the Routh-Hurwitz criterion.

Calculation of the equilibrium

First, the equilibrium U e , such that Üe = 0 and Ue = 0, is calculated. Since the tangential velocities of the mass are zero, we have :

   R e t = µR e n R e b = 0 (12)
where R e t and R e b are the quasi-static frictions forces according to directions t and b respectively and R e n is the quasi-static normal reaction. We now have to solve the following system:

   KU e = F + (µP t + P n ) R e n U e z = 0 (13)
which leads to the following equations :

       U e = µR e n cos θ k 1 sin θ k 2 0 T R e n = F n 1 -µ (cos θ tan α + sin θ tan β) (14)
The last equation of system [START_REF] Moirot | Friction and Instabilities, chapter Some examples of friction-induced vibrations and instabilities[END_REF] explains the relation between R e n and the system factors and requires that R e n > 0 to ensure the contact. Therefore, recalling that F n > 0, a critical friction coefficient µ e below which the quasi-static sliding equilibrium exits is given by :

µ e = 1 cos θ tan α + sin θ tan β (15)
We notice that for the proposed model, only the angles α, β and the friction direction θ have an influence on the equilibrium existence.

Stability around the equilibrium

The stability of the system is studied for small perturbations δU around the equilibrium U e of the system (i.e. δU + U e ).

By linearising the third equation of system (8), the following equation is obtained:

   δR t δR b    = ∂    R t R b    ∂R n U=0, Rn=R e n δR n + ∂    R t R b    ∂ U U=0, Rn=R e n δ U = µ    1 0    δR n - µR e n V     1 0 0 1   -    1 0    ⊗    1 0         δ Ut δ Ub    =    µδR n -µ R e n V δ Ub    (16)
In V and depending on the sliding equilibrium. This contribution, which increases when the external force increases or the sliding velocity decreases, can be taken into account through a damping matrix C b . The system dynamic is now given by :

   mδ Ü + (C + C b ) δ U + KδU = (µP t + P n ) δR n δU z = 0 (17) with C b = µP b T R e n V P b (18)
Finally, a unique formulation using displacements, velocities and accelerations in the tangential plane can be obtained by substituting the normal reaction in the problem . The first step uses the fact that the normal displacement is zero for expressing δR n according to U p = δU x δU y T :

δR n = k 1 tan α k 2 tan β U p (19)
Then, by using this equation, the motion equation projected into the tangential plane can be written as:

m Üp + C + Cb Up + K + K µ U p = 0 (20)
where C, Cb and K are the restrictions of the matrices C, C b and K to the tangential plane, respectively. Matrices m Ĩ, C and K are diagonal. They characterise the two uncoupled modes which determine the sliding system behaviour without friction. Moreover, we have :

K µ = -µ PT t k 1 tan α k 2 tan β (21)
where Pt = cos θ sin θ is the restriction of P t to the tangential plane. The friction effect is linked to matrices Cb and K µ which lead to mode coupling. It can be noted that the matrix Cb is symmetric while matrix K µ is non-symmetric.

Furthermore, mode coupling depends on the system geometry with an inherent normal/tangential coupling due to angles α and β and on the friction direction in relation to the system which is characterised by the angle θ.

The system stability is analysed by calculating the complex eigenvalues λ of the characteristic equation :

det λ 2 m Ĩ + λ C + Cb + K + K µ = 0 (22)
If all the real parts of the complex eigenvalues are negative, then the system equilibrium is stable. However if, one real part becomes positive, the equilibrium becomes unstable.

To determine the stability areas of the system equilibrium according to different factors, the Routh-Hurwitz criterion can also be used [START_REF] Meinsma | Elemantary proof of the routh-hurwitz test[END_REF]. One of the main advantages of the Routh-Hurwitz criterion is to provide analytical expressions versus the stability analysis of the mechanical system under study. The development of equation ( 22) yields to a polynomial equation of degree 4:

λ 4 + a 3 λ 3 + a 2 λ 2 + a 1 λ + a 0 = 0 ( 23
)
where λ are the eigenvalues of the linearised mechanical system. The application of Routh-Hurwitz criterion to this characteristic equation gives the four following coefficients a, b, c and d

a = a 0 (24) b = a 3 (25) c = a 2 a 3 -a 1 (26) d = a 1 (a 2 a 3 -a 1 ) -a 0 a 2 3 ( 27
)
Generally speaking, if all these coefficients are positive, the equilibrium is stable. When at least one of the four coefficients is negative, the equilibrium is unstable. In our physical practical case, it may be observed that a > 0, b > 0 and c > 0. Thus, equilibrium stability is governed only by coefficient d: if d > 0, then the equilibrium is stable (i.e.

∀ i (λ i ) < 0). If d < 0 the equilibrium is unstable. By using the Routh-Hurwitz criterion, the critical friction coefficients µ c such that d (µ c ) = 0 can be determined. The stability of the equilibrium changes when µ = µ c . Depending on different cases, the critical friction coefficient are solutions of a polynomial equation of degree 1 (planar friction without damping), degree 2 (rectilinear friction with or without damping) or degree 4 (planar friction with damping). In the last two cases, if there are two positive solutions, the equilibrium can become stable again after being unstable. All the mathematical developments and results are given in Appendix A.

In the following section, the stability of the mechanical system will be analysed by calculating the complex eigenvalues λ of the characteristic equation defined in equation 22.

Results

This paper presents the results for two specific cases. The first one concerns the undamped model. It not only provides an overview of the influence of some physical parameters but also offers the opportunity to discuss the elementary effect of a planar or rectilinear friction. The second case concerns the damped model and undertake the extension of the destabilisation paradox in the presence of planar friction.

The study is performed for a reference undamped system with the following parameters : m = 0.01 kg,

k i = mω 2 i for i = 1, 2 with ω 1 = 100 × 2π rad s -1 , ω 2 = 75 × 2π rad.s -1 , α = 30 • , β = 60 • , F n = 10 N, θ = -30 • and V = 5 m.s -1 .
Parameters differing from these reference values and addition of damping will be explicitly specified.

4.1 Undamped model with planar or rectilinear friction: parameter influence on stability area and complex eigenvalues

In this section a parametric study on the stability of the undamped system is carried out. The influence of the equilibrium, the system stiffness and the sliding direction on the stability area (i.e. the boundary between the stable and unstable zones) are treated.

Stability area: presentation and understanding

This paragraph introduces the graphical representation of the stability areas used in the paper. It also gives first results showing possible differences in the system stability depending on whether planar or rectilinear friction is considered.

Figure 2 shows the difference between planar and rectilinear friction in the stability area of the undamped system. This graphic has two variable parameters : the direction of the sliding velocity θ and the friction coefficient µ. For each couple (θ, µ), the boundary between the stable and unstable zones is determined (in blue for rectilinear friction and in red for planar friction). In addition, all the stable and unstable zones are labelled with the wording "stable" or "unstable" in blue (rectilinear friction) or red (planar friction). For friction coefficients µ > µ e (θ) (see equation ( 15)), there is no sliding equilibrium: this is indicated by "No equilibrium" in green. To properly understand the reading of this graph, let us illustrate with two examples. At first we consider the couple (-45 • , 0.5) indicated by a circle in Figure 2. It is observed that the mechanical system is unstable in both friction cases. In a second example, we consider the couple (-10 • , 0.2) marked with a star in Figure 2. In this case, the system with planar friction is unstable whereas the system with rectilinear friction is stable.

More generally, regarding the stable area versus all the possible couples (θ, µ), it can be noticed that the undamped system is less stable with planar friction than with rectilinear friction. We can also observed that the boundary of instability is not symmetric with respect to the diagonal line defined by θ = -45 • . This is not surprising since there is no geometrical symmetry in the system (i.e. the stiffness k 1 and k 2 as well as the angles α and β are different).

[Figure 2 about here.]

Influence of stiffness

In this section, the influence of the stiffness is studied. The stiffness sum is fixed using reference value k 1 + k 2 = 156.25 × (2π) 2 N.m -1 whereas three values of stiffness difference are selected:

k 1 -k 2 = 8.75 × (2π) 2 N.m -1 ,
43.75 × (2π) 2 N.m -1 and 78.75 × (2π) 2 N.m -1 , the middle one corresponding to the reference case. Figures 3(a 3(a) and (b) point out that the difference between k 1 and k 2 has a real impact on the system stability. The stable and unstable areas can drastically change according to variations in these parameters. The effect is rather similar with planar or rectilinear friction. We observe that the lower the stiffness difference, the less stable the system equilibrium (i.e. the stable area increases on the graph versus the couple (θ, µ)). This can be easily explained by the fact that the difference between the mode frequencies that can coalesce is directly linked to the difference between the stiffness k 1 and k 2 : if the difference between k 1 and k 2 is low, the difference between the frequencies is also low and coalescence between the associated modes is easier. Thereby, mode coupling instabilities appear at a lower friction coefficient µ as illustrated in Figures 3(c) and(d).

Figures 3(c),(d), (e) and (f) clearly illustrate the difference between planar and rectilinear friction. In the case of rectilinear friction, the classical behaviour of an undamped system is obtained (see Figures 3(d) and (f)). As long as the system is stable, the real parts are equal to zero. Once the system becomes unstable, real part variations versus the friction coefficient are symmetrical with respect to the real part axis (i.e. symmetrical about zero) and a perfect coalescence pattern is obtained (meaning that the frequencies of the unstable mode and the associated stable mode are exactly the same, as previously indicated in [START_REF] Sinou | Mode coupling instability in friction-induced vibrations and its dependency on system parameters including damping[END_REF][START_REF] Hoffmann | Effects of damping on mode-coupling instability in friction induced oscillations[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode coupling patterns in case of brake squeal[END_REF][START_REF] Kirillov | Non conservative stability problems of modern physics[END_REF]). In contrast, in the case of planar friction, the coalescence pattern is not perfect (the unstable mode and the associated stable mode have not the same frequency in Figure 3(c)). Likewise, in Figure 3(c), the real part variations are more complex due to the additional damping term provided by planar friction (see Section 3.2).

Table 1 compares the friction coefficients µ c at the Hopf bifurcation point and the frequencies of the unstable mode for planar or rectilinear friction. It turns out that planar friction also induces a modification of the Hopf bifurcation point.

[Table 1 about here.]

[Figure 3 about here.]

Influence of sliding velocity

The influence of the sliding velocity on the stability of the undamped system with planar friction is now studied. As previously seen in Section 3.2, the consideration of planar friction provides an additional damping term (i.e. matrix Cb ) resulting from the linearisation of the friction force direction. The amplitude of this term depends on sliding velocity V .

In the case of rectilinear friction, this contribution does not exist. To determine the influence of this damping contribution on the system stability, a complex eigenvalue analysis is carried out using three values of the sliding velocity V = 1 m.s -1 , 5 m.s -1 and 10 m.s -1 . For conciseness, the stability areas are not shown and only the eigenvalues graphs for θ = -30 • are presented. Indeed, the variation of the sliding velocity has no effect on the stability areas (i.e. 4(a) shows that the system becomes unstable at a lower friction coefficient with planar friction than with rectilinear friction. Thus, the consideration of planar friction can contribute to the destabilisation paradox: an increase in damping via the additional damping contribution Cb may tend to make the system unstable for a lower value. As already mentioned, the sliding velocity has no influence on the stability boundary: in each planar friction configuration, the system becomes unstable for the same values of the friction coefficient at Hopf bifurcation point µ c (see Figure 4(a) for µ = µ c = 0.168 and Table 2). However, the sliding velocity drastically affects the real part variations. Different behaviours can be observed according to the friction coefficient range.

For µ in [0; µ c ], a decrease of the sliding velocity lowers the real part. The same behaviour is observed for higher friction coefficients (i.e. µ > 0.25). However for friction coefficients in [µ c ; 0.25], a decrease of the sliding velocity increases the value of the real part. In the whole range, when V is large (i.e. c b = R e n V is small), the behaviour of the system tends to the rectilinear friction case.

Finally, Figure 4(b) shows the variation of the associated frequencies according to the friction coefficient µ for the three values of the sliding velocity. A decrease of the sliding velocity V leads to a damping increase provided by planar friction and has an influence on mode coalescence. A decrease of the sliding velocity increases the phenomenon of imperfect coalescence. This is notably true at the Hopf bifurcation point as indicated in Table 2. This is a direct consequence of the influence of damping (and more precisely the "smoothing effect") and supports existing studies on the subject [START_REF] Hoffmann | A minimal models for studying properties of the mode-coupling type instability in friction induced oscillations[END_REF][START_REF] Sinou | Mode coupling instability in friction-induced vibrations and its dependency on system parameters including damping[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode coupling patterns in case of brake squeal[END_REF][START_REF] Hervé | Extension of the destabilization paradox to limit cycle amplitudes for a nonlinear self-excited system subject to gyroscopic and circulatory actions[END_REF].

[Figure 4 about here.]

[Table 2 about here.]

Influence of sliding velocity direction

The influence of the sliding velocity direction is studied more deeply in this section. The variations of the real parts, the frequencies and the complex plane are displayed for two directions (θ = -30 • and -60 • ). The sliding velocity is equal to 5 m.s -1 in both cases. For the reader's comprehension, it is recalled that the stability areas (i.e. the boundaries between stable and unstable zones) for the undamped system with planar friction (and the direction of the sliding velocity between -180 • and 180 • ) are given in Section 4.1.1.

Figure 5(a) shows that the sliding velocity direction affects the system stability: for θ = -30 • , the real parts of the eigenvalues become positive beyond a lower friction coefficient than for θ = -60 • . Figure 5(b) shows the frequency variations and outlines the occurrence of a crossing of both modes involved in the coalescence pattern for θ = -60 • .

Finally, one major effect of the sliding velocity direction can be seen in Figure 5(c) which depicts the variation of the modes in the complex plane (variations of the real and imaginary parts of eigenvalues). For θ = -30 • , the unstable mode at the Hopf bifurcation point has the higher frequency (i.e. the mode with a frequency of 100Hz for µ = 0) whereas the mode with the lower frequency corresponds to the associated stable mode. On the contrary, for θ = -60 • , the unstable mode at the Hopf bifurcation point corresponds to the mode with the lower frequency (the mode with a frequency of 75Hz for µ = 0) and beyond the mode coalescence, the crossing of the stable and unstable modes causes the unstable mode to have higher frequency. This phenomenon is proper to planar friction. In the case of rectilinear friction, the unstable mode corresponds to the mode with the higher frequency for both values of θ.

[Figure 5 about here.]

Influence of damping for a system with planar or rectilinear friction

One main original contribution of this section is to extend well-known results proposed by many of researchers on the effect of structural damping in the field of friction-induced noise and vibration by considering the specific case of a system with planar friction. It will illustrate that the "lowering effect" (that tends to stabilise the mechanical system with rectilinear friction if the two modes involved in the coalescence are equally damped) and the "smoothing effect" (that may act in an unintuitive way by destabilising the system [START_REF] Hoffmann | A minimal models for studying properties of the mode-coupling type instability in friction induced oscillations[END_REF][START_REF] Sinou | Mode coupling instability in friction-induced vibrations and its dependency on system parameters including damping[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode coupling patterns in case of brake squeal[END_REF]) are significantly affected by planar friction. To achieve this goal, the influence of damping on system stability is analysed for planar or rectilinear friction. Comparisons with results of Section 4.1 (i.e. undamped system with planar or rectilinear friction) and differences in the stability of the damped system depending on friction assumption will be undertaken. Damping is defined by default as η = 0.02 and c = mηω 1 . Firstly, the influence of an iso-damping (η 1 = η 2 ) is studied for the mechanical system with planar or rectilinear friction. Effects of a non-iso-damping is discussed then. Finally, the influence of the sliding velocity for the damped system is highlighted.

Influence of iso-damping

A stability analysis of the iso-damped system is carried out for three different damping coefficients values: η = 0.01, 0.02 and 0.05. Figures 6(a) and (b) depict the stability areas for planar friction and a rectilinear friction, respectively. Firstly, we observe that the boundaries between the stable and unstable zones for the damped system with planar or rectilinear friction are quite close even if some differences can appear. For the undamped system, the results on the stability areas obtained with planar or rectilinear friction were very different (see Section 4.1). Thus, the contribution of a structural isodamping in the system completely changes the previous comments about the impact of the choice of friction. Moreover, the extent of the stable zones increases with damping η for all the values of θ between [0; -90] • . This reflects the fact that adding proportional damping stabilises the system for both planar and rectilinear frictions.

To deepen the influence of iso-damping, the variation of the real and imaginary parts of the eigenvalues are plotted in Figures 6(c) and(e), andFigures (d) and (f) for planar and rectilinear frictions respectively and θ = -30 • . First of all, as previously explained by some authors [START_REF] Hoffmann | A minimal models for studying properties of the mode-coupling type instability in friction induced oscillations[END_REF][START_REF] Sinou | Mode coupling instability in friction-induced vibrations and its dependency on system parameters including damping[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode coupling patterns in case of brake squeal[END_REF] for systems with rectilinear friction, increasing damping η makes the eigenvalue real parts decrease, whereas the frequencies remain approximately constant, as illustrated in Figures6(d) and (f). In addition, since the two modes involved in the coalescence are equally damped, the coalescence pattern is perfect. This phenomenon of iso-damping, called the lowering effect, extends to the mechanical system with planar friction : increasing η enlarges the stability range since the real parts are lowered evenly and become positive for a higher value of the friction coefficient, as indicated in Figure 6(c). Likewise, changes in frequencies are invariant with respect to isodamping η. The imperfect coalescence of modes is only due to the presence of the additional damping term Cb that is provided by planar friction. Table 3 assesses the friction coefficient at the Hopf bifurcation point and the frequency of the unstable mode for the different cases.

Thus, the only noticeable difference between rectilinear and planar friction seems to be the sensibility of the boundaries between the stable and unstable areas: an increase of iso-damping η has a larger effect on the stability area with planar friction (see Figures 6(a) and(b)). This can easily be explained by considering the slope of the real parts versus the friction coefficient. In the case of planar friction, the slope of the real part is lower due to the presence of the additional damping term Cb . This leads to a stronger modification of Hopf bifurcation point µ c in comparison with rectilinear friction, as illustrated in Figures 6(c) and(d). However it must be recalled that the additional damping term due to planar friction can also contribute to the destabilisation paradox and may tend to make the system unstable for a lower value of the friction coefficient than with rectilinear friction as previously shown in Section 4.1.3.

In conclusion, due to the addition of iso-damping, the stability areas obtained with planar friction are close to those obtained with rectilinear friction. This can be globally explained by the fact that the added damping matrix C (i.e. contribution of iso-damping) limits the effects of the contribution of matrix Cb due to planar friction. We recall, however, that these results obviously depend on the set of physical parameters. If the contributions of matrices Cb and C are changed and the contribution of Cb becomes much more important than that of C (e.g. by changing sliding velocity V or iso-damping η), larger differences between planar and rectilinear friction may occur (as described in Section 4.1.1).

[Table 3 about here.]

[Figure 6 about here.]

Influence of non-proportional damping

The influence of a non-proportional damping is now studied for planar and rectilinear friction. One main contribution of this section is to investigate the extension of the well-known destabilisation paradox due to the presence of nonproportional damping and the associated phenomena such as the lowering and smoothing effects to the specific case of a system with planar friction. The stability areas are plotted for several values of η x with η y = 0.02 constant such as ηy ηx = 0.5, 1.0 and 2.0 and for several values of η y with η x = 0.02 constant such as ηx ηy = 0.5, 1.0 and 2.0. Figures 7(a) and (c) show the stability areas for planar friction whereas Figures 7(b) and (d) correspond to the stability areas for rectilinear friction. The stability areas for planar or rectilinear friction and several distributions of non-proportional damping are quite close even if some differences are noticeable. For example, in the case of planar friction, it is noticed that the variation of η x or η y has some effect on stability for sliding velocity directions close to x-axis (0 • or 180 • ) or to y-axis (90 • or -90 • ). In the case with rectilinear friction, no changes are observed.

Compared with the results obtained with an iso-damping, the system stability in regard to the pair of variables (θ, ηx ηy ) is not straightforward. Whereas in the case of an iso-damping stability areas changes were similar in regard to adding damping whatever the values of θ, in the case of a non-proportional damping, variation of the system stability is more complex, especially for planar friction (Figures 7(a) and (c)). For example, an increase of the damping ratio ηy ηx with η y constant decreases the stability for θ = [0;-54.4] • but slightly increases the stability for θ = [-54.4; -90] • (see Figure 7(a)). On the contrary, an increase of the ratio damping ηx ηy with η x constant has no noticeable effect on the stability for θ = [0; -45] • whereas it decreases the stability for θ = [-45; -90] • (see Figure 7(c)). Thus, even if it is possible to increase the stability of the non proportional damped system in a specific sliding direction by modifying the damping ratio, this may decrease the stability of the system in another direction.

[Figure 7 about here.]

To explore and develop these initial observations on the role of non-proportional damping more deeply and to investigate the changes brought about by the additional damping contribution through planar friction, we propose to plot the variations of the real parts and frequencies for two specific values of θ (θ = -30 • and θ = -60 • ). The results are displayed in Figures 8. The variation of the eigenvalues in the complex plane are also plotted in Figure 9.

First of all, the classical results on the effect of non-proportional damping for mechanical system with rectilinear friction are found again [START_REF] Hoffmann | A minimal models for studying properties of the mode-coupling type instability in friction induced oscillations[END_REF][START_REF] Sinou | Mode coupling instability in friction-induced vibrations and its dependency on system parameters including damping[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode coupling patterns in case of brake squeal[END_REF]: if damping is spread non-equally over the two modes, shifting and smoothing effects can be seen on the variations of the real part and the coalescence curves: whereas the variation of the real part only occurs at the Hopf bifurcation point in the case of iso-damping η y = η x (see Figure 6(d)), the real parts associated with stable and unstable modes increase and decrease, respectively, with the variation of the friction coefficient in the case of non-proportional damping η y = η x (see Figure 8(e)). Then, if the smoothing effect prevails, added damping acts in an unintuitive way by destabilising the system, as illustrated in Figure 8(e): the system with non-proportional damping (i.e. ηy ηx = 0.5 and η y = 0.02 in Figure 8(e)) becomes unstable from a lower value of the friction coefficient than in the case of the system with iso-damping (i.e. ηy ηx = 1 and η y = 0.02 in Figure 8(e)) while damping η x increases (for the same value of η y ). As a consequence, increasing the gap in damping between the two modes tends to decrease the stability range as previously explained by Hoffmann and Gaul [START_REF] Hoffmann | Effects of damping on mode-coupling instability in friction induced oscillations[END_REF] and Fritz et al. [START_REF] Fritz | Investigation of the relationship between damping and mode coupling patterns in case of brake squeal[END_REF]. Moreover, the merging between the two modes is imperfect in the case of non-proportional damping as illustrated in Figure 8(f): the difference between the frequencies of the stable and unstable modes depends on the structural damping ratio ηy ηx . Table 4 gives the friction coefficients at the Hopf bifurcation point and the corresponding frequency of the unstable mode in all configurations. Figure 9 allows us to track the unstable mode on the coalescence patterns via the complex plane. Results reveal that the unstable mode is the least damped mode. Thus, the unstable mode may be the mode that initially has the highest or lowest frequency value depending on the ratio between the coefficients η x and η y (see [START_REF] Sinou | The role of damping and definition of the robust damping factor (rd-factor) for a self-exciting mechanism with constant friction[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode coupling patterns in case of brake squeal[END_REF] for more details). In our case, when ηy ηx = 0.5 the unstable mode corresponds to the one with a frequency of 75Hz for a friction coefficient equal to zero. When the ratio between η y and η x is reversed, the unstable mode is the mode with a frequency of 100Hz for a coefficient of friction equal to zero.

We now discuss in details the case of the non-proportional damped system with planar friction and the notion of destabilisation paradox with the concept of lowering and smoothing effects in this case. Figures 8(a-d) and 9(a,b) display the real parts and frequencies for two specific values of θ (θ = -30 • and θ = -60 • ) and the corresponding variation of the eigenvalues in the complex plane, respectively. Table 4 gives the friction coefficient at the Hopf bifurcation point and the corresponding frequency of the unstable mode. Firstly, for a sliding velocity direction close to the x-axis (i.e. θ = -30 • ), an increase of the ratio ηy ηx (for a fixed value of η y = 0.02) generates a more unstable system (i.e. the Hopf bifurcation occurs at a lower value of the friction coefficient). This result is consistent with the destabilisation paradox for a system with rectilinear friction. However, when θ = -60 • (which means that the sliding velocity direction is close to the y-axis), the Hopf bifurcation point is substantially independent of the ratio ηy ηx . This result does not follow the trend of the classical notion of the destabilisation paradox found with rectilinear friction. The introduction of the additional damping Cb due to planar friction alters the classical conclusions concerning the effect of the non-proportional damping. The changes on the system stability remain small but are not negligible if one is interested in making a mechanical system more stable through adding modal damping. The sole consideration of the modal damping ratio (which is a classical result for a mechanical system with rectilinear friction) is not enough to ensure a robust design versus the damping effect and the contribution of the damping induced by planar friction has to be taken into account.

We discuss now the mode coupling phenomena more specifically. For ηy ηx = 1, the imperfect merging (i.e. difference between frequencies of the stable and unstable modes), that can be seen in Figure 8(b) and (d), depends only on the additional damping Cb due to planar friction. For all the other ratios (i.e. ηy ηx = 1) this frequency gap is a combination of both the additional damping Cb due to planar friction and the non-proportional damping ηy ηx . Figures 9(a) and (b) reveals an unexpected effect due to the presence of the additional damping Cb . As illustrated for two sliding velocity directions (θ = -30 • and θ = -60 • in Figures 9(a) and (b) respectively), the unstable mode is not directly related to the ratio ηy ηx contrary to a classical result with rectilinear friction. In fact, the sliding direction seems to be the main factor determining which mode becomes unstable. For θ = -30 • , the unstable mode corresponds to the mode with a frequency of 100Hz for µ = 0 (the mode whose frequency is higher for µ = 0). On the contrary, for θ = -60 • the unstable mode is the mode with a frequency of 75Hz for µ = 0 (the mode whose frequency is lower for µ = 0), even if the unstable mode frequency becomes larger than the stable mode frequency after the imperfect coalescence. Thus, it can be concluded that the additional damping contribution due to planar friction can drastically modify the coalescence patterns. The extension of the well-known destabilisation paradox due to the presence of non-proportional damping is not so intuitive and phenomena such as the lowering and smoothing effects are more complex. Considering planar friction, one of the key parameters that alters the conclusions obtained with rectilinear friction is the sliding velocity direction.

All these results show that the classical trends pointed out for a system with rectilinear friction cannot be extrapolated to a system with planar friction. Even if the smoothing effect is still observed (i.e. the fact that added damping may destabilise the system, as previously discussed and illustrated in Figure 8(a)), planar friction can either reinforce or compensate the destabilisation paradox (see for example the comparison between the two configurations shown in 8(a) and (c)). Thus, the effect of non-iso-damping for a system with planar friction is complex and the extension of the destabilisation paradox needs to be done in conjunction with the role of key parameters related to planar friction.

[Figure 8

Conclusion

In this paper the stability of a 3-D model with planar or rectilinear friction has been studied for different configurations by considering the influence of the equilibrium, the system's stiffness and the sliding direction.

In a first phase, the stability results for the undamped system with planar or rectilinear friction have shown significant differences in the stable and unstable areas. The critical friction coefficients at the Hopf bifurcation point are lower in the case of planar friction, which may be interpreted as a destabilisation due to the damping term introduced by the linearisation of the friction force direction. This damping leads to a smoothing of the bifurcations curves, i.e. a nonperfect mode coalescence. By increasing the sliding speed, the damping term is smaller and the complex eigenvalues corresponding to planar friction are close to those obtained with the rectilinear assumption, except in the neighbourhood of the Hopf bifurcation point. The effect of stiffness is similar for planar and rectilinear friction : the mode coupling effect is mainly affected by the proximity of the frequencies of the modes without friction. Compared with the rectilinear friction case, the sliding direction has a greater influence on the coalescence pattern: for some directions, a crossing of the two modes involved in the instability is observed.

In a second phase, the combined effect of the damping in the structure and the planar friction has been discussed. In terms of stability areas, the introduction of damping brings the results obtained with planar or rectilinear friction closer.

However, this obviously depends on the ratio between the different damping terms. An extension of the destabilisation paradox due to the presence of damping and phenomena such as the lowering and smoothing effects has also been undertaken by considering the specific case of a system with planar friction. It was demonstrated that the additional damping contribution due to the consideration of planar friction can drastically modify the coalescence patterns and the associated instability phenomena are not so intuitive. In addition to the damping ratio and the friction coefficient, the sliding velocity direction appeared as one of the predominant factors.

A Calculation of Routh-Hurwitz coefficients

In this section, the Routh-Hurwitz coefficients are calculated in the general case (for the system with planar friction and damping). The simplified expressions for the undamped system with planar friction and the system with rectilinear friction with damping are given.

First of all, after calculation the coefficients a 3 , a 2 , a 1 and a 0 of Equation 23 are given by

a 3 = mc x + mc y + µmc b m 2 (28)
a 2 = µ c x c b cos 2 θ + c y c b sin 2 θ -mk 1 tan α cos θ -mk 2 tan β sin θ + mk 1 + mk 2 + c x c y m 2 (29)
a 1 = -µ 2 (c b k 1 cos θ tan α + c b k 2 sin θ tan β) + µ c b k 2 sin 2 θ + c b k 1 cos 2 θ -c x k 2 tan β sin θ -c y k 1 tan α cos θ) + c x k 2 + c y k 1 ) /m 2 (30)
a 0 = -µ (k 1 k 2 tan β sin θ + k 2 k 1 tan α cos θ) + k 1 k 2 m 2 (31)
where

c b = R e n V . b 1 = ((c x -c y ) k 2 + (c y -c x ) k 1 ) c b k 2 m + (-2c y -2c x ) c x c y c b k 2 + (-c y -c x ) c 2 y c b k 1 sin 2 θ + (2k 2 -2k 1 ) m + c x c y + c 2 x c x c y k 2 tan β sin θ + (((c x -c y ) k 2 + (c y -c x ) k 1 ) c b k 1 m + (-c y -c x ) c 2 x c b k 2 + (-2c y -2c x ) c x c y c b k 1 cos 2 θ + (2k 1 -2k 2 ) m + c 2 y + c x c y c x c y k 1 tan α cos θ + -c x k 2 2 + (c y + c x ) k 1 k 2 -c y k 2 1 c b m -c 2 x c y c b k 2 -c x c 2 y c b k 1 (39) b 0 = -(k 2 -k 1 ) 2 c x c y m + -c x c y -c 2 x c x c y k 2 + -c 2 y -c x c y c x c y k 1 (40)
In the general case, if the coefficients a, b and c are strictly positive, the value of the friction coefficient at the Hopf bifurcation µ c such that d(µ c ) = 0 is a solution of a polynomial equation of degree 4. This yields to more than one solution with a positive value of the friction coefficient which means that the mechanical system can be stable again after being destabilised by a mode coupling instability.

A.2 Planar friction without damping

The case of a system with planar friction and without damping (c x = c y = 0) is now studied. The Routh-Hurwitz coefficients are :

a = -µ (k 1 k 2 tan β sin θ + k 2 k 1 tan α cos θ) + k 1 k 2 m 2 (41) b = µmc b m 2 (42) c = µc b k 2 cos 2 θ + k 1 sin 2 θ m 2 (43) d = µ 2 c 2 b (k 1 -k 2 ) cos θ sin θ m 4 ((k 1 -k 2 ) cos θ sin θ -µ (k 1 tan α sin θ -k 2 tan β cos θ)) (44)
In this case the value of the friction coefficient at the Hopf bifurcation µ c such that d(µ c ) = 0 is the solution of a polynomial equation of degree 1 and does not depend on c b . There is only one value of µ c for each sliding direction θ, which means that contrary to the system with damping, the system without damping remains unstable when it is destabilised by a mode coupling instability.

A.3 Rectilinear friction with damping

The case of a system with rectilinear friction and damping (c b = 0) is now studied. The Routh-Hurwitz coefficients are : In this case the value of the friction coefficient at the Hopf bifurcation µ c is the solution of a polynomial equation of degree 2. If the two solutions are real positive, the system can be stabilised again after being destabilised by a mode coupling instability.

a = -µ (k 1 k 2 tan β sin θ + k 2 k 1 tan α cos θ) + k 1 k 2 m 2 (

A.4 Rectilinear friction without damping

In the case of a system with a rectilinear friction and without damping, the system stability has to be studied by calculating the eigenvalues instead of using the Routh-Hurwitz criterion. Details of this method can be found by instance in [START_REF] Moirot | Étude de la stabilité d'un équilibre en présence de frottement de Coulomb[END_REF]. The friction coefficient at the Hopf bifurcation µ c is solution of a polynomial equation of degree 2. Thus, the system can be stabilised again after being previously destabilised by a mode coupling instability.

  the case of rectilinear friction, component δR b is neglected and Coulomb's law is directly written again. With planar friction, the term -µ R e n V δ Ub , resulting from the linearisation of the friction force direction, provides a velocity contribution. This additional term, occurring only in the case of planar friction, corresponds to a new damping contribution characterised by the viscous damping coefficient c b = R e n

  ) and (b) show the stability areas for the system with planar friction and rectilinear friction respectively. At θ = -30 • , the Hopf bifurcation (i.e. the variation of the real part of the eigenvalues with the friction coefficient) and the variation of the frequencies are represented in Figures3(c), (d), (e) and (f) for both planar and rectilinear friction.

Figures

  Figures 3(a) and (b) point out that the difference between k 1 and k 2 has a real impact on the system stability. The

  µ c does not depend on c b as shown by the expression of coefficient d(µ) in section A.2 of appendix A). Results on the boundaries of stable and unstable zones are identical to those presented in the previous section 4.1.1 (Figures 2).

Figures 4 (

 4 Figures 4(a) and (b) display the variations of the real and imaginary parts for θ = -30 • respectively (versus the friction coefficient for the three values of the sliding velocity). First of all, Figure 4(a) shows that the system becomes unstable

45) b = mc x + mc y m 2 2 = k 2 2 tan 2 β sin 2 θ + k 2 1 tan 2 α cos 2 θ c x c y m + c 2 y + c 2 x k 1 k 2 m 2 x k 2

 222222 (46) c = µ (-k 2 c y tan β sin θ -k 1 c x tan α cos θ) m + (c y k 2 + c x k 1 ) m + c x c y (c x + c y ) m 3 (47) d = µ 2 b 2 + µb 1 + b 0 m 5 (48) where b tan α tan β sin θ cos θ (49) b 1 = -c x c y (k 2 -k 1 ) 2m + c x c y + c tan β sin θ + (k 1 -k 2 ) 2m + c 2 y + c x c y k 1 tan α cos θ (50) b 0 = c x c y (k 1 -k 2 ) 2 m + c x c y + c 2 x k 2 + c 2 y = c x c y k 1 (51)

Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :

 23456789 Figure 2: Stability area using the reference parameters

Table 1 :

 1 Friction coefficient and unstable mode frequency at the Hopf bifurcation point for θ = -30 • (undamped model) Sliding velocity (m.s -1 )

			Planar Friction		Rectilinear Friction
	(N.m -1 )	µ c	unstable frequency (Hz)	µ c	unstable frequency (Hz)
	8.75 × (2π) 2	0.029	89.12	0.043	87.81
	43.75 × (2π) 2 (reference) 0.168	94.22	0.222	88.31
	78.75 × (2π) 2	0.347	96.44	0.412	87.58

• k 1 -k 2

Table 2 :

 2 Friction coefficient and stable and unstable mode frequencies at the Hopf bifurcation point for planar friction, θ = -30 • and different sliding velocities (undamped model)

	Friction	damping	µ c	unstable frequency (Hz)
	Planar	η = 0.01 0.195	92.70
	Planar	η = 0.02 0.209	91.88
	Planar	η = 0.05 0.236	90.74
	Rectilinear η = 0.01 0.222	88.27
	Rectilinear η = 0.02 0.224	88.28
	Rectilinear η = 0.05 0.231	88.31

Table 3 :

 3 Friction coefficient and unstable mode frequency at the Hopf bifurcation point for θ = -30 • (iso-damped system)

	Velocity direction	Friction	Damping ratio η y /η x	µ c	unstable frequency (Hz)
	θ = -30 •	Planar	0.5	0.227	90.50
	θ = -30 •	Planar	1.0	0.209	91.88
	θ = -30 •	Planar	2.0	0.192	93.00
	θ = -30 •	Rectilinear	0.5	0.219	86.30
	θ = -30 •	Rectilinear	1.0	0.224	88.28
	θ = -30 •	Rectilinear	2.0	0.216	90.39
	θ = -60 •	Planar	0.5	0.206	89.30
	θ = -60 •	Planar	1.0	0.208	89.98
	θ = -60 •	Planar	2.0	0.209	90.49

Table 4 :

 4 Friction coefficient and unstable mode frequency at the Hopf bifurcation point for η y = 0.02 (non-iso damped system)

Acknowledgments

The authors would like to thank the financial support provided by Celya (Centre Lyonnais d'Acoustique, Université de Lyon). J-J. Sinou acknowledges the support of the Institut Universitaire de France.

A.1 General case : planar friction with damping

The coefficients of the Routh-Hurwitz criterion are now given in the general case :

List of Figures


Mohand Bentobache
email: mbentobache@yahoo.com

Ahcène Bounceur
email: ahcene.bounceur@univ-brest.fr

Reinhardt Euler
email: reinhardt.euler@univ-brest.fr

Salvador Mir
email: salvador.mir@imag.fr

Yann Kieffer
email: yann.kieffer@esisar.grenoble-inp.fr

Minimizing test frequencies for linear analog circuits: new models and efficient solution methods

Keywords: set covering problem, consecutive-ones property, analog circuit testing, linear programming, interval graphs

This work presents new approaches to minimize the number of test frequencies for linear analog circuits. The cases of single and multiple fault detection regions for multiple test measures are considered. We first address the case when the injected faults have a single detection region in the frequency band. We show that the problem can be formulated as a set covering problem with a matrix having the consecutive-ones property for which the network simplex algorithm turns out to be very efficient. A second approach consists in modeling the problem by means of an interval graph, leading to its solution with a specific polynomial-time algorithm. A case-study of a biquadratic filter is presented for illustration purposes. Numerical simulations demonstrate that the two different approaches solve the optimization problem very fast. Finally, the optimization problems arising from multiple detection regions are modeled and solution approaches are discussed.

Introduction

Testing of analog circuits to verify their functionality is a time consuming and expensive task. In order to reduce costs, testing for the presence of faults that affect the device structure leads to more structured and cheaper test sets. Although the number of faults to be tested may indeed be very large, a reduced set of tests that are simple to apply may be sufficient to achieve high fault coverage values. There is today a pressing need from the Semiconductors Industry to provide these fault coverage measures and optimized tests, especially in the context of SoC (System-on-Chip) devices that embed digital and mixed-signal blocks in a single chip.

Multi-frequency tests (i.e., multi-tone sinusoidal signals) have been classically considered for the test and diagnosis of linear analog devices such as analog filters. Since the effect of parametric and catastrophic faults varies as a function of frequency, it is possible to derive a minimal set of test frequencies either for the detection or for the diagnosis of all potential faults. To optimize the set of test frequencies, the approaches based on sensitivity analysis have typically addressed parametric faults [START_REF] Alippi | Automated selection of test frequencies for fault diagnosis in analog electronic circuits[END_REF][START_REF] Grasso | A method for the automatic selection of test frequencies in analog fault diagnosis[END_REF]. However, these approaches are not accurate for very large deviations, such as those that result from catastrophic faults. On the other hand, the approaches based on fault simulation can handle catastrophic faults, but at the expense of very time consuming simulations when realistic faults at transistor-level are considered. Today, it is evident that fault simulation of analog circuits is becoming essential in order to optimize test sets relying on new techniques to accelerate fault simulation.

In this context, this work proposes a new technique for the optimization of multifrequency tests for linear analog circuits. Fault simulation is used to obtain the frequency intervals for the detection of each fault. New efficient algorithms are then presented for the selection of the optimal set of test frequencies within these intervals for the detection of all faults. A simple case-study is used to illustrate the algorithms. Numerical simulations with randomly generated problem instances demonstrate the good time complexity of the proposed algorithms, with a large improvement over previous approaches [START_REF] Mir | Fault-Based ATPG for linear analog circuits with minimal size multifrequency test sets[END_REF]. We notice that the test optimization algorithms are general, applicable to any case-study requiring an optimization of multi-frequency tests based on fault simulation data. This technique is today feasible for analog filters, but it is also applicable to other analog devices such as analog-to-digital converters or radio-frequency front-ends requiring multi-tone tests, provided that fault simulation data are made available.

This work is organized as follows: in Section 2, the mathematical formulation of the Set Covering Problem (SCP) is reviewed and some definitions are given. In Section 3, we present the mathematical formulation of the problem of minimizing the number of frequency intervals necessary to detect the faults of an analog circuit, and we study the specific structure of the related coefficient matrix. In Section 4, we present two approaches for solving the problem: the Linear Programming (LP) approach and the interval graph approach, and a numerical example is given for illustration purposes. In Section 5, we present a mathematical formulation of the problem to minimize the number of test measures necessary for detecting the faults of an analog circuit. In Section 6, we present a case-study of testing a biquadratic filter. In Section 7, we carry out a large-scale numerical study in order to compare both approaches and to evaluate their time complexity. Finally, Section 8 concludes the work and provides some perspectives.

The Set Covering Problem

In this section, we present the set covering problem, we cite the existing methods for solving this problem and some of its applications. Finally, we present the corresponding mathematical model and some definitions.

Overview

The Set Covering Problem is one of the most important models in combinatorial optimization. Indeed, a wide range of real-world problems can be modeled as SCP, namely: railway crew scheduling, airline crew scheduling, facility location, etc. The practical importance of the SCP has motivated many researchers to develop efficient algorithms and heuristics for finding good solutions in reasonable time. We can cite: exact algorithms based on branch-and-bound or branch-and-cut [START_REF] Balas | A Dynamic Subgradient-Based Branch-and-Bound Procedure for Set Covering[END_REF], greedy heuristics [START_REF] Chvatal | A greedy heuristic for the set-covering problem[END_REF], Lagrangian-based heuristics [START_REF] Beasley | A Lagrangian Heuristic for Set Covering Problems[END_REF], genetic algorithms [START_REF] Beasley | A Genetic Algorithm for the Set Covering Problem[END_REF], etc. In order to test the efficiency of the proposed algorithms, there exists a library of SCP Benchmarks (Beasley's OR Library [START_REF] Beasley | OR-Library: Distributing Test Problems by Electronic Mail[END_REF]).

The set covering problem is known to be NP-hard [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. However, there exist some particular forms of the SCP which are polynomial-time solvable such as the SCP with a constraint matrix having the consecutive-ones property (i.e., the ones in each row appear consecutively). This particular case, denoted in the following by SCP-C1P, can be solved efficiently with LP algorithms.

Some definitions

Let M = {1, . . . , m} and N = {1, . . . , n} be two sets of indices. Let A = (a ij , i ∈ M, j ∈ N ) be a binary (m × n)-matrix and c = (c j , j ∈ N ) be an integer n-vector. The value c j , j ∈ N, may represent the cost of column j. We assume without loss of generality that c j > 0 for all j ∈ N . We say that a column j ∈ N covers

a row i ∈ M if a ij = 1. A set S ⊆ N is called a cover if each row i ∈ M is covered by at least one column j ∈ S.
The problem of finding a cover S of minimum cost can be formulated as an ILP (Integer Linear Programming) problem as follows:

min z = c T x, ( 1a
) subject to Ax ≥ 1 m , x ∈ {0, 1} n , ( 1b
)
where 1 m is the m-vector of ones. When c j = 1, for all j ∈ N , problem (1) is called a unicost set covering problem.

Definition 1. ([10]).

• A block of 1's (block of 0's) in a row of a binary matrix A is a maximal set of consecutive 1-entries (0-entries) in this row.

• A binary matrix has the strong consecutive-ones property (strong C1P) if in every row the 1's appear consecutively, that is, if every row contains at most one block of 1's.

• A binary matrix has the consecutive-ones property (C1P) if its columns can be permuted in such a way that the resulting matrix has the strong C1P. If an ordering for the columns of a binary matrix yields the strong C1P, it is called a C1-ordering.

• A graph G is an interval graph if its vertices can be mapped to intervals on the real line such that two vertices are adjacent if and only if their corresponding intervals have non empty intersection.

• A binary matrix is totally unimodular if every square submatrix has determinant 0, 1, or -1.

Minimization of test frequencies to detect all the faults of an analog circuit

In [START_REF] Mir | Fault-Based ATPG for linear analog circuits with minimal size multifrequency test sets[END_REF], the problem of selecting a minimum number of test frequencies is formally presented and solved using boolean expressions. However, when the number of faults is large, the technique used is not efficient. In this section, we present a mathematical formulation of the general problem as a set covering problem. We show that, under some special conditions on the considered faults, the coefficient matrix of the SCP will have the C1P. Moreover, we suggest an interval graph formulation of the problem.

Mathematical formulation of the problem as an SCP-C1P

Consider a given linear analog circuit C. Let F = {F 1 , F 2 , . . . , F m } be the set of all faults which can occur in C. Let T be a test measure which will be used in the fault detection process. In order to detect a given fault F i using the test measure T , test signals with maximum amplitude and at different frequencies in the interval [f min , f max ] are used as inputs of the analog circuit. The fault is detected if the test measure T exceeds a fixed threshold τ . Fault simulation allows then to compute, for each fault F i , the frequency intervals for which the threshold is exceeded. These intervals are called the detection regions of fault F i . Note that each fault can have one, two or more detection regions. Let n i be the number of detection regions of the fault F i and R i = {R i1 , R i2 , . . . , R in i } be the set of all the detection regions of the fault F i . Note that these detection regions are disjoint: for each two detection regions R ip and R iq , we have R ip ∩ R iq = ∅. After that, we sort in increasing order the bounds of the different detection regions, we find the vector of frequencies:

f = (f 0 , f 1 , . . . , f n ), with f 0 < f 1 < f 2 < • • • < f n-1 < f n . ( 2
)
Then, we compute a set of n frequency intervals I = {I 1 , I 2 , . . . , I n } as follows:

I 1 = [f 0 , f 1 [, I 2 = [f 1 , f 2 [, . . . , I n = [f n-1 , f n [. (3)
We denote by A = (a ij , i = 1, . . . , m, j = 1, . . . , n) the matrix which is defined as follows:

a ij = 1, if ∃k ∈ {1, . . . , n i }: I j ⊆ R ik ; 0, otherwise. (4)
The problem consists in finding a minimal-cardinality set of frequency intervals S ⊆ I which can detect all the possible considered faults. One frequency for each selected frequency interval (typically in the middle of the interval) can be used in the optimized test set. This problem can be formulated as a unicost set covering problem: for each S ⊆ I and each frequency interval I j , j ∈ {1, . . . , n}, we define a binary variable as follows:

x S j = 1, if I j belongs to the set S; 0, otherwise. (5)
Now, if x = (x j , j = 1, . . . , n) is any vector of binary variables, then the mathematical model corresponding to this problem will be given by:

min r = 1 T n x, (6a) s.t. Ax ≥ 1 m , x ∈ {0, 1} n . ( 6b
)
Remark 1. The matrix A cannot have more than

2 m i=1 n i columns. Proposition 1. If n i = 1, for i = 1, . . . m, i.e.
, each fault has a unique detection region, then the matrix A will have the strong consecutive-ones property.

Proof. Let F i , i ∈ M, be a given detected fault and R i1 = [α i , β i ] be the unique detection region of fault i. By construction, R i1 can be written as a union of k consecutive intervals I j1 , I j2 , . . . , I j k , where α i is the lower bound of I j1 and β i is the upper bound of I j k . (Two intervals are consecutive if they have a common bound). Since the intervals I j r correspond to the consecutive columns a j r of A for r = 1, . . . , k, the detection regions correspond to the rows of A, and

I j r ⊆ R i1 , i.e., a ij r = 1 for r = 1, . . . , k.
Note, that the previous proposition is very important because it gives a simple sufficient condition for the matrix arising in our application to have the strong C1P.

In the following we consider only problems with faults having one detection region. Following Proposition 1, ILP problem ( 6) is a unicost set covering problem with a coefficient matrix having the strong C1P (SCP-C1P). More precisely, it is given by:

min z = n j=1 x j , (7a) s.t. rx(i) j=lx(i) x j ≥ 1, i ∈ M, x j ∈ {0, 1}, j ∈ N, ( 7b
)
where for each row i, lx(i) denotes the leftmost index l for which a il = 1 and rx(i) the rightmost index r for which a ir = 1. To get the LP-relaxation of the above problem, we simply exchange the integrality constraints against the non-negativity constraints [START_REF] Matousek | Understanding and Using Linear Programming[END_REF]. Thus, we get the following LP problem:

min z = n j=1 x j , (8a) s.t. rx(i) j=lx(i) x j ≥ 1, i ∈ M, x j ≥ 0, j ∈ N. ( 8b
)

Graph formulation of the problem

Previously, we have formulated the problem of minimizing the number of frequency intervals to detect all the faults of an analog circuit using an ILP model. In this subsection, we suggest a new formulation based on the concept of interval graph: we denote the detection region of fault i by the interval [a i , b i ] for i ∈ M . Let the interval graph G = (F, E F ) be defined as follows:

F = {[a i , b i ], i ∈ M }, E F = {F i F j :]a i , b i [∩]a j , b j [ = ∅}. ( 9
)
In the next section we will suggest a polynomial algorithm using this interval graph for solving the considered problem.

Approaches for solving the minimization problem of frequency intervals

In this section, we suggest two approaches for solving the problem of minimizing the number of frequency intervals necessary to detect all the faults of an analog circuit.

LP Approach

Let us recall the following results:

Theorem 1. ( [START_REF] Ghouila-Houri | Caractérisation des matrices totalement unimodulaires[END_REF]).

An (m × n)-matrix A with entries 0, 1 and -1 is totally unimodular if and only if each collection of columns from

A can be partitioned into two column sets such that in each row the sum of the entries of the first set and the sum of the entries of the second set differ by at most 1.

Theorem 2. ([13]).

Let A be an m × n integral matrix. Then the polyhedron defined by Ax ≤ b and x ≥ 0 is integral for every integral vector b if and only if A is totally unimodular.

Remark 2. ([10]).

Any matrix A having the C1P fulfills the conditions of Theorem 1 and, hence, is totally unimodular.

Following Theorem 2 and Remark 2, any basic feasible solution of the LP problem (8) represents a cover for the SCP-C1P [START_REF] Beasley | A Genetic Algorithm for the Set Covering Problem[END_REF]. Therefore, we can find an optimal solution using LP algorithms such as the primal or dual simplex method [START_REF] Dantzig | Linear Programming and Extensions[END_REF], the support method [START_REF] Bentobache | A two-phase support method for solving linear programs: Numerical experiments[END_REF], the hybrid direction algorithm [START_REF] Bibi | A hybrid direction algorithm for solving linear programs[END_REF], etc. However, it is more efficient to transform the SCP-C1P into a min-cost network flow problem [START_REF] Dom | Recognition, Generation, and Application of Binary Matrices with the Consecutive-Ones Property[END_REF].

Let us make the Veinott-Wagner transformations [START_REF] Veinott | Optimal capacity scheduling[END_REF] for the variables of problem (8): we introduce the variables y j , j = 1, . . . , n + 1 such that x j = -y j + y j+1 , j = 1, . . . , n. Hence, we obtain the following equivalent LP problem:

min z = -y 1 + y n+1 , (10a) s.t. -y lx(i) + y rx(i)+1 ≥ 1, i ∈ M, ( 10b
) -y j + y j+1 ≥ 0, j ∈ N. ( 10c
)
The dual of the above problem has m + n variables v 1 , v 2 , . . . , v m+n and n + 1 constraints and it is given by:

min w = - m j=1 v j , (11a) s.t. m+n j=1 a 1j v j = -1, ( 11b
) m+n j=1 a ij v j = 0, i = 2, . . . , n, ( 11c
) m+n j=1 a (n+1)j v j = 1, v j ≥ 0, j = 1, . . . , m + n. ( 11d
)
where A = (a ij , i = 1, . . . , n + 1, j = 1, . . . , m + n), and A has exactly one 1 and one -1 in each column. Remark that the LP problem ( 11) is a min-cost network flow problem. Thus, it can be solved by the network simplex algorithm. The scheme of the LP approach to solve the problem is described in the following steps:

Step 1. Compute the frequency intervals necessary to detect the different faults using relation ( 2 Step 5. Solve the min-cost network flow problem with the network simplex algorithm.

Interval graph approach

In order to solve the problem using the interval graph G defined by ( 9), we suggest the following algorithm: Algorithm 1 • Sort the vertices F i , i = 1, m by increasing order of their upper bounds b i ;

• Compute the frequencies f k , k = 0, n with (2); • Set S = ∅, label(i) = 0, i = 1, m; • For i = 1 to m • If label(i) = 0, then • Find the index k such that b i = f k and set S = S ∪ {[f k-1 , f k [}; • For j = i + 1 to m If F j is adjacent to F i in G, then set label(j) = 1; endif endfor endif endfor
The input of this polynomial algorithm is the interval graph G and the output is a minimal-cardinality set S of frequency intervals, which detects all the faults.

Note, that it is not necessary to construct the whole interval graph: we can use only its vertices, i.e., detection regions of the different faults, and replace the test "If F j is adjacent to F i in G" by the test "If a j < b i ≤ b j ". Remark 3. Contrarily to the LP approach, the interval graph approach uses the detection regions of the considered faults directly, i.e., without transforming the problem.

Numerical example

Example 1. We consider an analog circuit with five faults. The detection regions of these faults are:

R 1 = [1, 80], R 2 = [160, 1700], R 3 = [1, 1400], R 4 = [1400, 2000] and R 5 = [1000, 1500].
Let us compute the minimal-cardinality set of frequency intervals necessary to detect all the five faults with the two approaches. LP approach: First, we compute the different frequency intervals I j , j = 1, . . . , n: we sort the different bounds of the detection regions in increasing order, we get the vector f = (f i , i = 0, . . . , 7) = (1, 80, 160, 1000, 1400, 1500, 1700, 2000). Thus, the different frequency intervals are computed as follows: In order to answer this question, we solve the following ILP problem:

I 1 = [1,
min z = 1 n T x, (12a) s.t. Ax ≥ b, x ∈ {0, 1} n , ( 12b
)
where

x T = (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 ), b T = (1, 1, 1, 1, 1) and A =       1 0 0 0 0 0 0 0 0 1 1 1 1 0 1 1 1 1 0 0 0 0 0 0 0 1 1 1 0 0 0 1 1 0 0       . ( 13
)
It is easy to remark that the constraint matrix has the consecutive-ones property. Using the Veinott-Wagner transformations (x j = -y j + y j+1 , j = 1, 2, . . . , 7), we get the following equivalent LP problem:

min z = (c ) T y s.t. V y ≥ b , ( 14
)
where y ∈ R 8 , c = (-1, 0, 0, 0, 0, 0, 0, 1) T , b = (1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0) T and

V =                    
-1 1 0 0 0 0 0 0 0 0 -1 0 0 0 1 0 -1 0 0 0 1 0 0 0 0 0 0 0 -1 0 0 1 0 0 0 -1 0 1 0 0 -1 1 0 0 0 1 0 0 0 -1 1 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 -1 1

                    . ( 15
)
The dual problem of the LP problem ( 14) is

min w = -(b ) T v s.t. V T v = c , v ≥ 0 , ( 16
)
where v ∈ R 12 . Remark that the LP problem ( 16) is a min-cost network flow problem. The primal and dual optimal solutions obtained by the network simplex algorithm of CPLEX are:

v * = (0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0) T , y * = (0, 1, Interval graph approach: Let G = (F, E F ) be the graph defined by [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. First, we sort the vertices of the graph G in increasing order of their upper bounds, we get: The different iterations of the interval graph algorithm are as follows:

F 1 = [a 1 , b 1 [= [1, 80|, F 2 = [a 2 , b 2 [= [1, 1400[, F 3 = [a 3 , b 3 [= [
We set S = ∅ and label(i) = 0, i = 1, 5.

For i = 1, b 1 = 80, f k = f 1 = 80, S = {[f 0 , f 1 [} = {[1, 80[}. The unique adjacent vertex of F 1 is F 2 (the edge F 1 F 2 ∈ E F because a 2 < b 1 < b 2 ), so we set label(2) = 1.
For i = 2, we have label(2) = 1, so we pass to i = 3.

For i = 3, b 3 = 1500, f k = f 5 , S = {[f 0 , f 1 [, [f 4 , f 5 [} = {[1,

Case of multiple detection regions

Let us illustrate this case by the following numerical example:

Example 2.

Faults Detection regions

F 1 R 11 = [1, 80] F 2 R 21 = [160, 200], R 22 = [300, 1700] F 3 R 31 = [1, 1400] F 4 R 41 = [1400, 2000] F 5 R 51 = [1000, 1500]
The vector of frequencies computed by ( 2) is then f = (f 0 , f 1 , . . . , f 9 ) = (1, 80, 160, 200, 300, 1000, 1400, 1500, 1700, 2000).

The frequency intervals are:

I 1 = [1,

The constraint matrix of the SCP is then given by

A =      
1 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 1 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 1 1 0 0

      . ( 18
)
Let us remark that this matrix loses the consecutive-ones property because the second row does not have consecutive ones. Therefore, we cannot apply the techniques suggested above to solve it. If this case occurs in practice, the branch-and-cut algorithm of CPLEX can be applied to solve the ILP problem. For this example, CPLEX gives us the optimal solution: x * = (1, 0, 0, 0, 0, 0, 1, 0, 0) T . So a minimal-cardinality set of frequency intervals is

S = {I 1 , I 7 } = {[1, 80[, [1400, 1500[}. ( 19
)
Remark 4. Although the constraint matrix of the previous SCP does not possess the consecutive-ones property, solving the LP relaxation gives us the optimal solution x * = (1, 0, 0, 0, 0, 0, 1, 0, 0) T . This means that for this type of problems the LP relaxation can give a good approximate solution for the optimal one (good upper bound for the optimal value). That is why we expect the branch-and-cut algorithm of CPLEX (or any other solver) to produce an optimal solution in a few iterations in case that the optimal solution of the LP-relaxation is not optimal for the SCP. Further numerical experiments have to be carried out in order to confirm this assertion.

Minimization of the number of test measures necessary for detecting all the faults of an analog circuit

Previously, we have assumed that only one test measure is used in the testing process. However, in practice, often we use several test measures. So if k is the number of the used test measures, we need to solve k SCPs for solving the problem of frequency interval minimization. Therefore, solving the problem of minimization of test measures beforehand can dramatically reduce the CPU time of solving the whole optimization process.

Consider a given analog circuit C. Let F = {F 1 , F 2 , . . . , F m } be the set of all faults which can occur in C. Let T = {T 1 , T 2 , . . . , T k } be the set of k test measures to be used in the fault detection process. We denote by H = (h ij , i = 1, . . . , m, j = 1, . . . , k) the Fault-Test-Measure (FTM) matrix which is defined as follows:

h ij = 1, if the test measure T j detects the fault F i ; 0, otherwise. ( 20
)
The problem consists in finding the minimal-cardinality set of test measures which can detect all the considered faults. This problem can also be formulated as a unicost set covering problem: for each test measure T j , j = 1, 2 . . . , k and any set S ⊆ T , we define a binary variable as follows:

x S j = 1, if the test measure T j ∈ S; 0, otherwise. (21)
Hence, the mathematical model corresponding to this problem is given by:

min p = 1 T k x, s.t. Hx ≥ 1 m , x ∈ {0, 1} k . ( 22
)
The optimal solution for this integer linear programming problem constitutes a desired minimal-cardinality set of test measures which detect all the faults. Note that, in general, the coefficient matrix of this problem does not have a specific structure, so the branchand-cut algorithm of CPLEX can be used to solve this problem. Once a minimal-cardinality set of test measures determined, we can solve the problem of frequency interval minimization for each test measure in this set.

The global optimization scheme for testing a linear analog circuit with multiple test measures is shown in Figure 1.

Case-study

To illustrate our approach, similar to [START_REF] Mir | Fault-Based ATPG for linear analog circuits with minimal size multifrequency test sets[END_REF], we will now present a case-study biquadratic filter as shown in Figure 2. There are 6 test measures for this circuit that correspond to the common-mode signal at the input and at the output of each operational amplifier. For simplicity, only catastrophic (10 MOhm open and 1 Ohm short) faults in the passive We denote the test measures by T 1 , T 2 , . . . , T 6 and the faults by F 1 , F 2 , . . . , F 16 . Figure 3 shows the frequency behavior of test measure T 1 for some of these faults. The detection regions of each fault F i , i = 1, 2, . . . , 16 using test measures T j , j = 1, 2, . . . , 6 are computed using the fault simulator developed in [START_REF] Bounceur | CAO Platform for mixed circuit testing[END_REF][START_REF] Bounceur | CAT platform for analogue and mixed-signal test evaluation and optimization[END_REF]. Note, that the simulation results obtained in [START_REF] Bounceur | CAO Platform for mixed circuit testing[END_REF] indicate that test measures T 2 , T 4 and T 6 do not detect any fault. Hence, we only consider test measures T 1 , T 3 and T 5 . For simplicity also, we have only considered nominal simulations of the catastrophic faults. In practice, Monte Carlo simulations of each catastrophic fault should be considered, and worstcase detection regions be computed (that is, the intersection of the detection regions for each Monte Carlo instance).

First, we start by minimizing the number of test measures necessary for the detection of the sixteen faults. Following simulation results, test measure T 1 detects the faults: F 1 , F 2 , . . . , F 8 ; T 3 detects the faults F 9 , F 10 , F 11 , F 12 and T 5 detects the faults F 13 , F 14 , F 15 , F 16 . So the FTM matrix corresponding to the minimization problem of test measures is

H T =   1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1   . (23)
Observe, that since test measures T 2 , T 4 and T 6 do not detect any faults, the matrix H contains only three columns which correspond to the test measures T 1 , T 3 and T 5 . Without solving the SCP problem [START_REF] Bounceur | CAT platform for analogue and mixed-signal test evaluation and optimization[END_REF], it is clear that a minimal-cardinality set of test measures which detects all the faults is S = {T 1 , T 3 , T 5 }. Now, that a minimal set of test measures is at hand, we can proceed to solve the problem of minimization of test frequencies for each test measure belonging to this set. When we sort in increasing order the bounds of the different detection regions, we find the following vector of frequencies: (1, 85, 336, 647, 1014, 1566, 2685, 3442, 3732, 10 5 ) .

f = (f 0 , f 1 , . . . , f 9 ) =
Therefore, the frequency intervals are The constraint matrix of the SCP ( 6) is then: 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0

I 1 = [1,
A =             1 1
0 0 0 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 1             , ( 24
)
where each row i, i = 1, . . . , 8, corresponds to the detection region R i1 and each column j, j = 1, . . . , 9, to the frequency interval I j . Since the matrix A has the consecutive-ones property (each fault has a unique detection region), the optimal solution of the SCP (6) can be obtained using the network simplex algorithm. In our example, this leads to the following optimal solution: x * = (0, 0, 0, 1, 0, 0, 0, 0) T and r * = 1, i.e., a minimal-cardinality set of intervals detecting all the faults is S = {I 4 }. Therefore, it is sufficient to use a single test frequency belonging to the interval [647Hz, 1014Hz[ to detect all the faults F 1 , F 2 , . . . , F 8 .

Similarly, when we use the test measure T 3 , we detect the faults

The constraint matrix of the SCP is

A =     0 1 1 1 1 1 1 1 0 0 1 1 0 0 0 0 1 1 1 0     . ( 25
)
Hence, for the test measure T 3 , the solution of the SCP [START_REF] Beasley | A Lagrangian Heuristic for Set Covering Problems[END_REF]

The constraint matrix of the SCP is

A =     1 1 1 1 1 1 0 0 1 0 0 0 1 1 1 0     . ( 26
)
Thus, for the test measure T 5 , the solution is obvious again, that is, the frequency interval which detects all the faults F 13 , F 14 , F 15 and F 16 is I 1 = [1Hz, 1412Hz[, and a single test frequency belonging to this interval is required.

Numerical experiments

In order to compare the efficiency of the two approaches presented in Section 4 (LP approach and the interval graph one), we have implemented them in C++ programming language and carried out large-scale numerical experiments on a set of randomly generated test instances using an Intel(R) Core(TM)2 Duo CPU P8600 @ 2.40GHz machine with 4 GB of RAM.

We have generated 60 problems with number of faults m = 500, 000; 600, 000; 700, 000; 800, 000; 900, 000; 1, 000, 000 .

The detection region bounds are generated in the interval [1Hz, 10 5 Hz]. For each class of test problems with m faults, we generate ten problems. We have solved the different instances with the LP approach (LPA) using the network simplex method of the LP and ILP solver CPLEX [START_REF] Cplex | [END_REF] and the suggested interval graph algorithm (IGA). The CPU time of the two approaches IGA and LPA are reported in Table 1.

Let [a i , b i ] be the detection regions of the faults F i , i = 1, . . . , m. In the following, we give the implementation details for the LP approach:

Step 1. Compute the set of frequency intervals I as follows:

-sort the bounds of the detection regions in increasing order: let f be the vector of these sorted bounds; -delete duplicate elements from f , and let f = (f 0 , f 1 , . . . , f n );

-set I = {I j = [f j-1 , f j [, j = 1, . . .

, n};

Step 2. Compute the constraint (m × n)-matrix of the SCP-C1P as follows:

-set f max = f (n) and compute the vector t of dimension f max as follows:

-set t(k) = 0, for k = 0, . . . , f max; -set t(f (j)) = j, for j = 0, . . . , n; -set lx(i) = t(a i ) + 1 and rx(i) = t(b i ) for i = 1, . . . , m; Step 3. Make the Veinott-Wagner transformations: compute the constraint matrix V of the LP problem (10) as follows:

-for i = 1, . . . , m; for j = 1, . . . , m + n:

V ij =    -1, if j = lx(i); 1, if j = rx(i) + 1; 0, otherwise; (28)
-for i = m + 1, . . . , n + 1; for j = 1, . . . , m + n:

V ij =    -1, if j = i -m; 1, if j = i -m + 1; 0, otherwise; (29)
Step 4. Compute the constraint matrix A , the (n + 1)-vector of right-hand-sides b and the cost (m + n)-vector c of the min-cost network flow problem [START_REF] Matousek | Understanding and Using Linear Programming[END_REF]:

A = V T , b = (-1, 0, . . . , 0, 1) T , c = -1 m 0 n ;
Step 5. Solve the min-cost network flow problem with the network simplex algorithm.

Note, that the efficient implementation presented above computes the constraint matrices of the SCP-C1P and the min-cost network flow problem of the LP approach in small CPU times (less than 1 second on average for all the test problems). That is why we have not reported the CPU times of computing the constraint matrices of the SCP-C1P and the min-cost network flow problem in Table 1.

The average optimal values found by the two approaches IGA and LPA for the different test problems are shown in column OptVal. We compute the ratios of the CPU time of IGA over LPA. These ratios are shown in column CPU Ratio. Finally, the CPU time of the two approaches are plotted in Figure 4. The LPA shows a time complexity O(m), while IGA shows a time complexity O(m 2 ), where m is the number of faults.

The graphs of Figure 4 indicate that transforming the problem into a min-cost network flow problem and solving it with the network simplex algorithm is more efficient than the interval graph algorithm. Indeed, LPA is, on average, two times faster than IGA. However, the computation time is very small even for problems with an extremely large number of faults. Therefore, the interval graph algorithm can also be used in practice by test engineers because of its simple implementation; the good CPU times (an average of 16 seconds for optimizing the test of analog circuits with 1,000,000 faults); and the fact that it solves the original problem directly.

Conclusion

In this work, we have formulated as an SCP the problem of minimizing the number of test frequencies necessary to detect a set of faults injected into an analog circuit. We have shown that when the considered faults have a unique detection region, the constraint matrix of the SCP will have the strong consecutive-ones property. After that, we have reformulated this special case using interval graphs and an algorithm working directly with this graph is suggested. In order to solve the problem efficiently, two approaches are compared: an LP approach and an interval graph approach. The obtained numerical results show that the approach which transforms the problem into a mincost network flow problem to solve it by the network simplex algorithm is the most efficient. However, the interval graph approach can also be used by test engineers because it solves the original problem directly; its implementation is very simple and it is extremely fast with CPU times of a few tenths of seconds even for large-scale problems. Moreover, the optimization problem associated with multiple detection regions and multiple test measures is studied. Future work will focus on developing a branchand-cut algorithm for solving efficiently the set covering problem corresponding to the case of multiple detection regions.

); Step 2 . 4 .

 24 Compute the constraint matrix of the SCP-C1P using relation (4); Step 3. Make the Veinott-Wagner transformations. Let V be the constraint matrix of the LP problem (10); Step Compute the constraint matrix A , the (n + 1)-vector of right-hand-sides b and the cost (m + n)-vector c of the min-cost network flow problem (11): A = V T , b = (-1, 0, . . . , 0, 1) T , c = -1 m 0 n ;

  80[, I 2 = [80, 160[, I 3 = [160, 1000[, I 4 = [1000, 1400[, I 5 = [1400, 1500[, I 6 = [1500, 1700[ and I 7 = [1700, 2000[. The question is then: among all these seven intervals, which one must be included in a minimal-cardinality set?

T . ( 17 )

 17 Hence x * = (1, 0, 0, 0, 1, 0, 0) T , which means that a minimal-cardinality set of frequency intervals is S = {I 1 , I 5 } = {[1, 80[, [1400, 1500[}.

  1000, 1500[, F 4 = [a 4 , b 4 [= [160, 1700[ and F 5 = [a 5 , b 5 [= [1400, 2000[. The vector of frequencies computed by (2) is thenf = (f 0 , f 1 , . . . , f 7 ) =[START_REF] Alippi | Automated selection of test frequencies for fault diagnosis in analog electronic circuits[END_REF] 80, 160, 1000, 1400, 1500, 1700, 2000).

  80[, [1400, 1500[}. The adjacent vertices of F 3 are F 4 and F 5 , so label(4) = 1 and label(5) = 1. For i = 4, we have label(4) = 1, so we pass to i = 5. For i = 5, we have label(5) = 1. Since i = m, a minimal-cardinality set of frequency intervals is S = {[1, 80[, [1400, 1500[}.

  80[, I 2 = [80, 160[, I 3 = [160, 200[, I 4 = [200, 300[, I 5 = [300, 1000[, I 6 = [1000, 1400[, I 7 = [1400, 1500[, I 8 = [1500, 1700[, I 9 = [1700, 2000[ .

Fig. 1 .Fig. 3 .

 13 Fig. 1. Global optimization scheme for testing analog circuits with multiple test measures

  First, we start by minimizing the frequency intervals under the test measure T 1 : the faults detected using test measure T 1 are F 1 , F 2 , . . . , F 8 . The detection regions of the different faults are R 11 = [1, 10 5 ], R 21 = [85, 3732], R 31 = [85, 3732], R 41 = [1, 2685], R 51 = [1, 3442], R 61 = [336, 1566], R 71 = [1, 1014] and R 81 = [647, 10 5 ] .

  85[, I 2 = [85, 336[, I 3 = [336, 647[, I 4 = [647, 1014[, I 5 = [1014, 1566[, I 6 = [1566, 2685[, I 7 = [2685, 3442[, I 8 = [3442, 3732[ and I 9 = [3732, 10 5 [ .

F 9 ,

 9 F 10 , F 11 and F 12 . The corresponding detection regions are R 91 = [159, 7957[, R 10,1 = [1, 1740[, R 11,1 = [1, 1739[ and R 12,1 = [159, 7940[ . The different frequency intervals are I 1 = [1, 159[, I 2 = [159, 1739[, I 3 = [1739, 1740[, I 4 = [1740, 7940[ and I 5 = [7940, 7959[ .

R 13 , 1 =

 131 [1, 2798[, R 14,1 = [1, 1413[, R 15,1 = [1, 1412[ and R 16,1 = [1, 2794[ . The different frequency intervals are I 1 = [1, 1412[, I 2 = [1412, 1413[, I 3 = [1413, 2794[ and I 4 = [2794, 2798[ .

  is obvious, that is, a single test frequency belonging to the interval I 2 = [159Hz, 1739Hz[ detects all the faults F 9 , F 10 , F 11 and F 12 .Finally, with test measure T 5 , we detect the faults F 13 , F 14 , F 15 and F 16 . The corresponding detection regions are

Table 1 .

 1 CPU time of the different approaches

	m	OptVal IGA LPA CPU Ratio
	500,000 560.80 4.10 3.18	1.29
	600,000 611.90 5.30 3.52	1.51
	700,000 672.50 6.70 3.85	1.74
	800,000 704.90 9.40 4.30	2.19
	900,000 762.90 12.00 4.53	2.65
	1,000,000 809.10 16.30 4.76	3.42
	Mean	8.97 4.02	2.13

Fig. 4. CPU time of the two approaches: IGA and LPA