Sofiane Ghenna

Eric Vezzoli

Member, IEEE Christophe Giraud-Audine
email: christophe.giraud-audine@ensam.eu

Member, IEEE Frederic Giraud

Michel Amberg

Member, IEEE Betty Lemaire-Semail

B L Semail

Enhancing Variable Friction Tactile Display using an ultrasonic travelling wave

Keywords: Travelling wave, Friction Modulation, tactile stimulator and display. !

In Variable Friction Tactile Displays, an ultrasonic standing wave can be used to reduce the friction coefficient between a user's finger sliding and a vibrating surface. However, by principle, the effect is limited by a saturation due to the contact mechanics, and very low friction levels require very high vibration amplitudes. Besides, to be effective, the user's finger has to move. We present a device which uses a travelling wave rather than a standing wave. We present a control that allows to realize such a travelling wave in a robust way, and thus can be implemented on various plane surfaces. We show experimentally that the force produced by the travelling wave has two superimposed contributions. The first one is equal to the friction reduction produced by a standing of the same vibration amplitude. The second produces a driving force in the opposite direction of the travelling wave. As a result, the modulation range of the tangential force on the finger can be extended to zero and even negative values. Moreover, the effect is dependant on the relative direction of exploration with regards to the travelling wave, which is perceivable and confirmed by a psycho-physical study.

INTRODUCTION

T HE control of the friction between a finger and a flat surface can be performed through different techniques, mainly electrovibration and ultrasonic vibration. The first relies on the application of a high voltage between the finger and an insulated electrode. The electrostatic attraction which is produced increases the normal force and, consequently, the friction force [START_REF] Kaczmarek | Polarity effect in electrovibration for tactile display[END_REF], [START_REF] Bau | TeslaTouch: electrovibration for touch surfaces[END_REF]. The second uses a standing wave obtained at a resonance of a plate within the ultrasonic range. The finger sliding on it experiences a reduction of the friction coefficient if the amplitude of vibration is in the micrometer range. The effect has been firstly explained as a squeeze film effect [START_REF] Watanabe | A method for controlling tactile sensation of surface roughness using ultrasonic vibration[END_REF], [START_REF] Biet | Discrimination of Virtual Square Gratings by Dynamic Touch on Friction Based Tactile Displays[END_REF]; late evidence based on the fluidodynamic interaction explains the features of the friction reduction [START_REF] Sednaoui | Experimental evaluation of friction reduction in ultrasonic devices[END_REF] and points toward a mechanical explanation [START_REF] Vezzoli | Role of fingerprint mechanics and non-Coulombic friction in ultrasonic devices[END_REF] involving intermittent contact [START_REF] Wiertlewski | Partial squeeze film levitation modulates fingertip friction[END_REF]. The induced sensations of the two techniques are similar [START_REF] Vezzoli | Physical and Perceptual Independence of Ultrasonic Vibration and Electrovibration for Friction Modulation[END_REF] and the implementation of position tracking techniques allows to simulate simple texture rendering through the spatiotemporal modulation of the vibration [START_REF] Biet | Discrimination of Virtual Square Gratings by Dynamic Touch on Friction Based Tactile Displays[END_REF].

So far, these techniques do not produce an active force on the finger. It means that they need a relative motion of the fingertip on the plate to increase or decrease the friction. To address this limitation [START_REF] Dai | LateralPaD: A surfacehaptic device that produces lateral forces on a bare finger[END_REF] uses a surface haptic device which vibrates the touched surface simultaneously in both normal and lateral directions. The resulting elliptical surface motion produces the active force. The normal-bending and lateral-longitudinal resonances of a glass plate, and its driving piezoelectric system, are arranged to occur at the same frequency. The relative phase of these vibration modes is used to control the magnitude and direction of the lateral force. Alternatively, a travelling wave can be used instead of a standing wave. In [START_REF] Biet | A Piezoelectric Tactile Display Using Travelling Lamb Wave[END_REF], a travelling Lamb wave tactile display has been proposed, to indicate smooth or braking sensation, using the stator of a travelling wave motor. The authors modelled the lateral stretching force as proportional to the relative velocity between finger and the surface without measuring the generating force.

In this paper, we propose a method to create a travelling lamb wave in a finite beam to enhance the tactile interaction with a variable friction display. In the first part, we describe how to obtain the travelling wave, and how it can produce an active force. In the second part, experimental study and comparison between the use of standing and travelling wave in tactile display is provided, exploring the influence of vibration amplitude, applied normal force and speed. Finally, a psychophysical validation is presented.

THEORETICAL CONSIDERATIONS

Travelling wave generation

In [START_REF] Sashida | An introduction to ultrasonic motors[END_REF] two transducers were used to generate a travelling wave in a finite beam. One transducer is used to produce the travelling wave, while the other absorbs it to prevent the formation of standing wave. This method requires impedance matching between the beam and a transducer, as described in [START_REF] Kurosawa | High speed ultrasonic linear motor with high transmission efficiency[END_REF].

The second method [START_REF] Loh | An object transport system using flexural ultrasonic progressive waves generated by two-mode excitation[END_REF] is based on the excitation of two successive flexure modes of the beam. Indeed, the frequency response of a beam features several modes such as represented on Fig. 1 which highlights four modes of vibration.

By exciting the beam at an intermediate frequency between two modes (denoted f middle on the same figure) the Fig. 1. Simulation of the frequency response of an aluminium beam of 600 × 6 × 6 mm 3 for an excitation force amplitude of 0.5N at 19mm of the edge modes A and B will be promoted. By imposing the same amplitude and a phase shift of ±90 • on the two modes vibrations, a travelling wave is obtained [START_REF] Minikes | On the sensing and tuning of progressive structural vibration waves[END_REF]. In this work, we control the amplitude of each transducer in closed loop as explained in [START_REF] Ghenna | Modelling, identification and control of a Langevin transducer[END_REF]; the advantages are to control the direction amplitude and direction of the travelling wave without any preliminary fitting procedures as in [START_REF] Minikes | On the sensing and tuning of progressive structural vibration waves[END_REF], and to mitigate external disturbances such as the pressure of a user's finger [START_REF] Giraud | Vector control method applied to a traveling wave in a finite beam[END_REF], [START_REF] Ghenna | Modelling and control of a travelling wave in a finite beam, using multi-modal approach and vector control method[END_REF]. Figure 2 depicts a typical implementation of a system generating an ultrasonic travelling wave. In this implementation, the two piezoelectric transducers excite each mode in the normal direction. Therefore we proceed by giving the method to determine the vibration amplitude of the actuators. In steady state, the vibration amplitude on the beam w, depends on both the normalized mode shapes φ n (x), and their modal amplitude or modal participation factor W n where n is an integer, according to:

w(x, t) = W e jωt = ∞ n=1 W n φ n (x)e jωt (1)
W includes real and imaginary part of the vibration, at location x, ω is the angular frequency of vibration, and j 2 = -1. Considering two dominant modes of vibration A and B and the positions of the actuators on the beam x 1 and x 2 , the vibrations of the two actuators W 1 and W 2 are equal to the vibration of the beam at their position, which are given by:

W 1 = W (x 1 ) = W A φ A (x 1 ) + W B φ B (x 1 ) W 2 = W (x 2 ) = W A φ A (x 2 ) + W B φ B (x 2 ) (2)
Introducing matrix notation, (2) becomes:

W 1 W 2 = φ A (x 1 ) φ B (x 1 ) φ A (x 2 ) φ B (x 2 ) W A W B = φ x1,x2 W A W B (3)
To obtain a travelling wave, the two modes A and B have to be superimposed, with the same modal amplitude, and phase shifted by 90 • , as explained in figure 3. This leads to the following condition to obtain a travelling wave :

W A = ±jW B (4)
Hence, the reference signal for the transducers for an amplitude W writes:

W 1 (t) = real W (±jφ a (x 1 ) + φ b (x 1 ))e jωt W 2 (t) = real W (±jφ a (x 2 ) + φ b (x 2 ))e jωt (5)
The equations 5 shows that the references of each actuator do not necessarily have the same amplitude nor is their phase necessarily shifted by 90 • to obtain a travelling wave. The next section presents the modelling of the contact force produced by the travelling wave.

Contact force modelling produced by the travelling wave

The contact mechanism between the finger and the surface is modeled based on Fig. 5(a). When a pure travelling wave is generated along the beam, each point of the surface follows an elliptical trajectory as illustrated in Fig. 5(b). The horizontal component of the elliptical motion of particle leads to the generation of relative tangential speed between an object in contact and the beam as in [START_REF] Sashida | An introduction to ultrasonic motors[END_REF], while the perpendicular component provides conditions for generating a squeeze film effect. The maximum particle speed at the top of the rail accordingly to [START_REF] Maas | Simulation model for ultrasonic motors powered by resonant converters[END_REF] is given by:

S h = πW hω λ (6)
where W , h and λ are respectively the vibration amplitude, thickness of the beam and wavelength. For the • a driving force denoted F p which propels the object in the direction opposite to the travelling wave,

• a friction force which depends on the velocity of the object, and which is opposite to the direction of its motion, and is denoted F f . These considerations allow us to write the resulting tangential force F t of the beam on a finger according to the relative direction of motion with the regards to the direction of the travelling wave. In this paper, the direction of the travelling wave remains the same, and we define F t positive when it acts against the movement of the finger. Therefore, F p and F f are always positive, and F t should write differently according to the direction of motion:

F t+ = F f -F p (7) F t-= F f + F p (8)
where F t+ and F t-are tangential forces when the finger moves respectively in opposite direction and in direction of the travelling wave.

The contributions of each force is experimentally studied in the sequel and we present the experimental set-up in the next section.

EXPERIMENTAL PROTOCOL

Experimental set-up

To implement the generation of the travelling wave, a 600 × 6 × 6 mm 3 aluminum beam was chosen. With those dimensions and material, the working frequency was 28.3 kHz, corresponding to the midpoint between the frequency of flexural mode A and mode B, and the wavelength of the travelling wave was λ = 41 mm. The experimental test bench is presented in Fig. 4(a).

The beam is actuated by two Langevin transducers (FBL28252SSF-FC, Fujicera, Japan). The vibration amplitude of the two actuators is controlled through a DSP (TI 2812, Texas Instrument, USA) and two amplifiers (4051, HSA, Japan), the references for the control being previously calculated using Eq [START_REF] Sednaoui | Experimental evaluation of friction reduction in ultrasonic devices[END_REF]. The vibration amplitude of the beam was measured using an interferometer laser (OFV -5000, Polytech, Germany). The whole system (beam with transducers) is attached to the support as shown in Fig. 4(b). This support is fixed to the moving part of a Tribometer (TRB, CSI, Switzerland) allowing lateral friction force measurement. The parameters that can be selected are: the normal load and the linear sliding velocity. The apparatus is supported by a computer (data acquisition, monitoring, recording and analysis) software. The exploring length of the tribometer has been set to one wavelength due to its range limit of 60 mm, with a sinusoidal exploration speed. In the paper, the friction is not measured with a fingertip, but with a probe which is validated in the next section.

Probe Validation

The probe interacting with the travelling wave has to behave as the human finger in the ultrasonic domain. To reflect the fact that the fingertip is a soft material, the probe which size is 16 × 13 × 10 mm 3 , is made of a silicone core with a grooved surface and rounded edges, covered with surgery tape to reproduce the mechanical behaviour of the fingertip [START_REF] Bueno | Pile Surface Tactile Simulation: Role of the Slider Shape, Texture Close to Fingerprints, and the Joint Stiffness[END_REF].

In order to evaluate the used probe, we compared the friction modulation produced by a standing wave on it to the ones on a human finger [START_REF] Sednaoui | Experimental evaluation of friction reduction in ultrasonic devices[END_REF] and a TangoPlus probe [START_REF] Friesen | Bioinspired artificial fingertips that exhibit friction reduction when subjected to transverse ultrasonic vibrations[END_REF]. The results are presented in Fig. 7 showing a similar friction modulation in the ultrasonic domain.

F f [N] Friction force F f 0.5 N 1 N 1.5 N (

Ft

Probe Human finger TangoPlus Once the probe has been assessed, it has been used as a substitute to the human finger for the tribology measurements presented in the next section.

Identification of F p and F f

The experimental set-up allows to measure the friction force for several speed (v = 12.9 mm/s, v = 25.8 mm/s and v = 38.6 mm/s) and normal force (F = 0.5 N, F = 1 N and F = 1.5 N) conditions. We can also set different vibration amplitude from W = 0 µm to W = 1.6 µm. For each operating condition, we averaged the measurements over 10 test runs.

Using ( 7) and ( 8) it is possible to extract the forces F p and F f from F t+ and F t-by:

F p = F t--F t+ 2 (9) F f = F t-+ F t+ 2 (10)
A tribological investigation has been done in order to evaluate the lateral friction forces, by changing experimental conditions: vibration amplitude, applied normal force and the exploration speed. The obtained results show that a direct influence on the lateral friction forces can be affected by these parameters.

EXPERIMENTAL RESULTS

In figure 6

Influence of the vibration amplitude

The figures 6(b), 6(c), 6(e), 6(f) show that the driving force F p increases with W , while F f decreases. The figures 6(a) 6(d) clearly show an influence of the progressive wave on the tangential force. The effect can be so high that the probe is propelled by the travelling wave. This happens for vibration amplitude W above 1 µm, and for a displacement in opposite direction to the travelling wave. By contrast, when the probe moves in the same direction as the travelling wave, the lateral force remains almost unchanged, because F p compensate the drop of the friction force F f produced by the vibrations. This highlights the fact that two antagonist effects exist in the contact mechanism: for one mechanism, vibrations reduce the friction as demonstrated by figure 6(c) and 6(f), while the other increases the driving force which can be added or subtracted by choosing direction of propagation. Although we could not obtain a higher apparent friction than the original one of the beam, we could cancel it. This is an improvement on a device based on standing waves only, which has a saturation on the friction reduction, as demonstrated in [START_REF] Vezzoli | Role of fingerprint mechanics and non-Coulombic friction in ultrasonic devices[END_REF].

Influence of the normal force

As shown in [START_REF] Sednaoui | Experimental evaluation of friction reduction in ultrasonic devices[END_REF], the normal force applied by the finger has a direct influence on the friction reduction for a device based on a standing wave. This also the case here, and the driving force F p and the friction force F f are both concerned, however, the main effect is on the friction force which is consistent with a coulombic force.

Influence of the velocity

The velocity of the probe has an influence on the tangential force F T in both directions. However, F p is almost unchanged, and thus, it can be concluded that the effect is mainly on the friction F f which again shows a coulombic force, but also presents a viscous component.

Psychophysical study

The perception of the consequences of the force applied by the travelling wave has been evaluated with a psychophysical experiment, performed using two-alternatives forcedchoice method described in [START_REF] Jones | Application of Psychophysical Techniques to Haptic Research[END_REF]. The assumption under test is that this force induces different friction reduction according to the relative speed of the travelling wave and the finger. Six participants (5 males, 1 female; age range: 19 -32 years) were asked to wash their hands and properly dry them before the experiment. Also, they were asked to touch the beam where a standing wave was produced in order to train them to friction reduction. Afterwards, during the test, only travelling waves were used, and their directions were randomly set. The exploration was allowed for 10 seconds for each case with a vibration amplitude of 1.2 µm. The participants were asked to indicate "which of the two directions had less friction". An answer was positively scored if the direction with the lowest tangential force was detected. The experience was repeated 10 times with a randomization of the trials. The results are summed up on figure 8. The participants were able to perceive the force applied by the travelling wave, with an average success rate of 81.66% and a standard deviation of 16%. The hypothesis test against the random answer was refused with a p < 0.05 leading to accept the conclusion that the direction of the applied force by a travelling wave was perceived. Indeed, Travelling wave

F t- F t+ F t
Standing wave Normalized tangential force Fig. 9. Low friction level with a travelling wave and with a standing wave at 1.6 µm the exercise is harder than telling the difference between "on" and "off", because the participants had to detect the change of friction which occurs with direction and also to detect the direction of propagation of the travelling wave. Hence, 100% is difficult to achieve. Considering this, we can conclude that when the finger explores the vibrating surface, it appears that one of the directions is more slippery in comparison with the opposite direction contrarily to the standing wave.

DISCUSSION

The first feature of the travelling wave used in a variable friction tactile display is the higher contrast between low and high friction levels. In Fig. 9, the average of F t over the different experimental runs at maximum vibration amplitude are presented to highlight this point. For comparison purposes, the case of a standing wave is also presented on the same figure. The standing wave is obtained on the same experimental setup by setting W B = 0, and W A is controlled so as to always have the same amplitudes for standing and travelling wave. The control used in the experimental test bench exactly allows to compare the two cases together. The difference between the possible range between the extreme values of the net force on the probe approaches 70% of the force without vibration. It is also worth noticing that a negative net force can be produced by the travelling wave.

Moreover, the travelling wave offers the possibility to differentiate the feeling according to the direction of the exploration, this is shown on the same figure and confirmed by the psychophysical study.

So far the existence of the two contributions of F p and F f has been postulated. has been clearly linked to a coulombic friction with viscous components, due to its dependence on the normal force and the velocity of the exploration. To further investigate this point, the figures 10 superimpose the force F f obtained previously and compare it to tangential force obtained with a pure standing wave.

The figures clearly show that the estimated force F f on a travelling wave is equal to the tangential force produced by a standing for similar operating conditions. This demonstrates that the tangential force in the case of a travelling wave is the result of the superimposition of two contributions :

•

The friction force which is equal to the effect of the standing wave which has been extensively studied,

•

The driving force which is specific to the travelling wave and can be attributed to the elliptical motion at the surface of the beam. This decomposition into two phenomena give some insight that could guide the design of tactile devices using travelling waves. Properties of the design such as material, geometry, wavelength, frequency could be exploited to promote the driving force F p over F f . For example Eq [START_REF] Vezzoli | Role of fingerprint mechanics and non-Coulombic friction in ultrasonic devices[END_REF] shows that it is necessary to increase beams thickness to promote the effect in the tangential axis. These points, however, require further investigation.

CONCLUSION

In this paper, the effect of an ultrasonic travelling wave on the friction force has been measured and the estimation of the applied force has been proposed for a finger pressure and an exploration speed compatible with touch exploration. The results confirmed the possibility to obtain an asymmetry of forces for a probe exploring the surface according to the direction of the travelling wave. The contrast of the friction is more important in the case of a travelling wave, compared to the standing wave. For extreme vibration amplitude, the tangential force can even reach a negative value and thus drag the finger. The influence of the velocity and the applied force has been reported.

Future work should focus on the design of new devices using a travelling exploiting two modes of vibration, with a larger exploring area, and the tracking of a user's finger to create spatio-temporal stimuli.

Fig. 2 .

 2 Fig.2. A typical implementation of a system, and the two modes produced by the langevin transducers.

Fig. 3 .

 3 Fig.3. The superposition of mode A and mode B at three different time stamps t 0 , t 1 , t 2 , and the resulting travelling wave.

Fig. 4 .

 4 Fig. 4. Experiment set-up schematic, a beam is actuated by two Langevin transducers, fixed to the moving part of a Tribometer (a) and picture (b)

  Fig. 5. (a) Principle of the force application by travelling wave. (b) Elliptical trajectory of the particles

Fig. 7 .

 7 Fig. 7. (a) Comparison between the measured friction force for the used probe and the human finger in consistent condition of applied force and exploration speed[START_REF] Sednaoui | Experimental evaluation of friction reduction in ultrasonic devices[END_REF]. (b) Comparison of the asymptotic friction reduction for the developed probe, human finger and the Tango plus probe[START_REF] Friesen | Bioinspired artificial fingertips that exhibit friction reduction when subjected to transverse ultrasonic vibrations[END_REF]

  we present the results for the different operating conditions; figures 6(a), 6(b) and 6(c) were obtained at a probe velocity of 25.8 mm/s and normal forces of F = 0.5 N, F = 1 N and F = 1.5 N while figures 6(d), 6(e) and 6(f) were obtained at a normal force of 0.5 N and probe velocity of v = 12.9 mm/s, v = 25.8 mm/s and v = 38.6 mm/s. Figures 6(a) and 6(d) represent the evolution of the Lateral force F t measured by the tribometer (plain line: opposite direction to the travelling wave, dashed line: same direction as the travelling wave), 6(b) and 6(e) represent the actual driving force F p obtained using equation 9 and 6(c), 6(f) represent the calculated force F f using equation 10.

  Fig. 8. Results of psychophysical testing

  The behavior of the friction force F f Fig.10. Tangential forces for a standing wave (plain line) compared to the friction force F f of the travelling wave (dashed line) for several normal force (top) and velocities (bottom) conditions.

			Friction force F f and F t
		0.5
						0.5 N	1 N	1.5 N
		0.4
	f and F t [N]	0.2 0.3
	F
		0.1
		0 0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6
			Vibration amplitude [µm]
				Friction force F f and F t
		0.25
				12.9 mm/s	25.8 mm/s	38.6 mm/s
		0.2
	f and F t [N]	0.1 0.15
	F
		0.05
		0 0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6
			Vibration amplitude [µm]

ACKNOWLEDGMENTS

This work has been carried out within the framework of the project StimTac of IRCICA (institut de recherche sur les composants logiciels et matériel pour la communication avancé), and the Project Mint of Inria. Authors thank FP7 Marie Curie Initial Training Network PROTOTOUCH, grant agreement No. 317100 for founding


Farhang Radjai

H Ochstleistungsrechenzentrum

Stéphane Roux

Friction-induced self-organization of a one-dimensional array of particles

We study the frictional motion of a linear array of parallel identical cylinders on a plane where all interparticle and particle-plane contacts can be described by a Coulomb friction law. The entire system can be characterized in its steady state by an effective coefficient of friction and an effective inertia, both varying with the external driving force. In the steady-state regime, the organization of the array appears as regions involving different rotation modes (sliding with no rotation, rolling on the support, alternate rotations, etc.). The mechanisms leading to this self-organization in modes are studied analytically.

I. INTRODUCTION

In the context of standard continuum mechanics, the only field needed to describe the elastic deformation of a solid body is the displacement field alone. Rotations do not appear. They can be computed indirectly from the rotational of the displacement field. Other elasticity theories have been developed in order to account for the rotational degrees of freedom of a solid material, such as "Cosserat" or "micropolar" theory [START_REF] Nowacki | Theory of Micro-polar Elasticity {Springer[END_REF]. Dimensional analysis shows that an intrinsic length scale À has to be introduced. This length may be interpreted as the scale below which rotations can significantly differ from those deduced from the rotation of the displacement field. At a larger scale, the difference vanishes, and thus upon coarse graining the standard theory is recovered.

In the case of granular materials, the discontinuous geometry of the solid phase, and thus of the displacement field, makes it difficult to establish a direct connection between a discrete description at the particle level and a continuum homogenization at a large scale. In particular, the effect of individual rotations appears to play a determining role [START_REF] Bardet | [END_REF]. Although the framework of Cosserat theory seems more suited to such media, it does not appear to provide a simple key to the analysis. Sorne attempts in this direction have been proposed [3]; however, basic questions remain, such as the identification of this intrinsic length scale, which is often postulated to be a typical particle size. Other effects may motivate the study of rotations in granular media, such as the propagation of sound [START_REF] Schwartz | a one-dimensional geometry, we expect similar qualitative effects in higher dimensionalities. A two-dimensional case is presently under study[END_REF]6].

• Also at Laboratoire de Physique et Mécanique des Milieux Héterogènes, ESPCI, 10 rue Vauquelin, F-75231 Paris Cedex 05, France.

Considering only the rotations of particles, one major difficulty appears immediately: the generally disordered geometry of the packing of particles induces a "frustration" of rotations. To mention only the two-dimensional case, as soon as a "loop" of contacting particles contains an odd number of elements, rotations of the particles in the loop are not possible without friction being mobilized on at least one contact. We employ the word frustration in analogy to the physics of spin glasses, where the presence of antiferromagnetic coupling between neighbors makes it impossible to satisfy ali couplings simultaneously. We will not try here to develop the comparison between these two fields any further , although sorne suggestions along these lines have been proposed [4].

In this paper we consider a very simple system consisting of an array of disks in mutual contact supported by a plane. Two neighboring particles are in contact with the plane and thus a pair of particles and the support form a loop of three solids such that rotations are frustrated. Since frustration involves sliding on at least one contact, it is important to incorporate the friction law accurately in the description of the system. We will consider here the most basic law of friction, namely, Coulomb's law of friction.

When the array is pushed in the direction of its mean orientation, particles will move and rotate and a collective organization of rotations will emerge. This organization is rather complex. In particular, it involves typical length scales w hi ch may be mu ch larger thau the parti cie size, depending on the coefficients of friction (interparticle and particle-plane) and the confining pressure. This observation might motivate revisiting discrete-continuum transition through the analysis of intrinsic length scales different from the particle size.

In order to present our approach, we will first describe in Sec. II the simple case of a single particle being pushed on a plane. In Sec. III, we will consider the more general case of a long array and report on sorne numerical simulations. In Sec. IV, we will study analytically the collective organization of the array and describe the different phases observed. We will discuss the transition between the different phases and the expected effects on the global motion of the entire system.

II. THE CASE OF A SINGLE PARTICL E

Consider first the geometry illustrated in Fig. 1: A single disk supported by a basal plane is pushed by a constant horizontal force N. The mass and radius of the particle are assumed to be unity. The acceleration of gravity is set to unity as well. All quantities are normalized with respect to these natural units of the problem and, in particular, the moment of inertia I is a mere geometrie constant. The friction between the particle and the support is assumed to obey Coulomb's law with a friction coefficient J.L 1 • The friction between the particle and the pushing block also follows Coulomb's friction law with a coefficient J.L• The translational and angular velocities of the particle are called v and w, respectively. We are interested in the steady-state regime, where the acceleration is constant. All solids are assumed to be rigid, and thus they cannot interpenetrate. Figure 1 gives the convention used for the orientation of the forces. The tangent force at the parti de-block contact is called T, whereas the friction force on the support is S. The normal reaction of the support is R.

The equations of the dynamics are the following:

v = N-S, 0 = R -1-T , lw = S -T, ( 1
)
where the fust two equations refer to Newton's law of motion in the horizontal and vertical directions, and the last one accounts for the rotation.

To find a complete solution to the equations of dynamics, we need two more equations prescribed by the friction law at the two contact points. A friction law is a relation between the friction force and the relative tangential velocity at the contact point between two solids. The graph of this relation for Coulomb's law is displayed in Fig. 2. In two dimensions, it is written

T J.1.N v -J.1.N 1 ------ FIG.
Vr = 0 ::} T E [-J.LN 1 J.LN ] 1 Vr > 0::} T = -J,LN, Vr < 0 ::} T = J.LN 1 (2)
where Vr is the relative tangential velocity and T and N are the tangential and normal components of the contact force. ln the case where Vr = 0 , there is no equation relating the tangential force to the normal contact force. In the present case of a single particle, we will consider different configurations with different relative velocities at the contact points and we will show that the friction force is always uniquely determined. In the next section, we will discuss a more general formulation of Coulomb's law adapted to dynamical systems [START_REF] Jean | Mechanics of Geometrical Interfaces[END_REF].

There are two particular values of rotation velocities:

w 1 = 0 and w2 = v. Depending on the relative position of w with respect to these particular values, the tangent force at the contacts will assume different values. However, up to a change in the origin of time, the position of w with respect to the key values Wi can be equivalently analyzed through the relative positions of the accelerations w with respect to Wi, since we are interested in the steady state, where the accelerations remain constant in time. Five different cases can be distinguished.

(1) w > v. In this case, we have T = J.LN and S = -J.L' R. Since both tangent forces, T and S, contribute to a torque opposite to w, the angular acceleration has to turn negative. Therefore, no steady-state solution can be found.

(2) w =v. In this case, T = J.LN, but S cannot be deduced directly from Coulomb's law. However, the kinematic constraint w = v suffices to determine the system entirely. The complete solution is written

T = J.LN, R = 1 + J.LN, s = Ut'î )N, (3)
.

-

(1-~-') N W - l+I .
For this solution to be acceptable, the condition ISI < J.L 1 R has to be satisfied. This gives a constraint on the value of N:

J.L'(1 + I) N < N 1 = ---''----'-----,-=------,- I + J.L -J.LJ.L'(1 + I ) (4)
(3) 0 < w < v. ln this case, T = J.LN and S = J.L 1 Ras given by Eq. ( 2). The complete solution of the system is s FIG.

2. Coulomb's law graph. (5) Compatibility of this solution with the conditions 0 < w < v sets the following constraints on the values of N:

T = J.LN, R = 1 + J.LN, S = J.L' (1 + J.LN) , v = -J.L' + (1 -J.L 1 J.L)N, w = t (J.L' -(1 -J.L 1 )J.LN) .
(6) (4) w =O. ln this case, S = J.L 1 R while T is not determined by the friction law directly. The condition w = 0, however, gives an additional constraint which leads to

T = _L_ 1 1-J-1.'' v = N-_E_ R 1 = 1-f'' 1-1-''' w = o. s = _E_ 1-1-''' (7)
The condition JTJ < J.LN implies [START_REF] Jean | Proceedings of the Contact Mechanics International Symposium[END_REF] (5) w < O. In this case, S = J.L 1 R and T = -J.LN. As in the first case, the two tangent forces contribute to a torque opposite to w and thus no steady-state solution is possible.

ln this way, the dynamics of the block-particle system show three possible steady-state phases: rolling (phase 1), rotating and sliding on the plane (phase 2), and sliding on the plane without rotating (phase 3). Figure 3 shows the friction forceS at the particle-plane contact as a function of the driving force N. The former increases linearly with the applied force in the two fust phases to saturate to a value greater than J.L 1 in the last one.

As long as the translational motion of the particle is of interest, these results can be compared to the dynamics of a rigid block of mass m, pushed by a force N. The respective equation of motion is written 3. Evolution of the particle-plane friction force as a function of the driving force N. Ali forces are normalized with respect to the weight of one particle (see text) . The first regime N < N1 corresponds to rolling on the support, the intermediate regime N 1 < N < N2 corresponds to partial rotation where both contacts are sliding contacts, and finally in the last regime N > N2, the particle slides on the support with no rotation. Here, J.L = 0.1 and J.L 1 = 0.2.

N-Mm = mv, (9) where M is the coefficient of friction at the black-plane contact. Now, a look at the steady-state solutions of the single-particle problem [Eqs. {3), [START_REF] Schwartz | a one-dimensional geometry, we expect similar qualitative effects in higher dimensionalities. A two-dimensional case is presently under study[END_REF] , and ( 7)) in different regimes shows that they can ail be written in the canon-ica} form of Eq. ( 9) if we introduce an effective inertia m* and an effective coefficient of friction M*:

{ i!±fl_ (T=Il) , m*- 1 - (1-1-''1-'), 1, N<N1 N1 < N < N2 N2 <N.
(10)

The "effective friction force" on the plane is given by M*m*. While it is obvious that the inertia m• should depend on whether the particle slides or rolls, it is surprising that the coefficients of friction J.L and J.L 1 enter the expression for the inertia. Note that the effective inertia is piecewise constant, and decreases by steps with N.

Similarly, the effective coefficient of friction M* also depends on the driving force in a piecewise constant manner, increasing with N. The fust case, which corresponds to a low applied force, is somewhat counterintuitive. The effective coefficient of friction is exactly zero, even though there is a nonzero friction force in the system. This friction force being proportional to the acceleration, it only contributes to the effective inertia rather than to the coefficient of friction M*. We will see below that this behavior remains qualitatively valid for an array of partides, with essentially the same low force behavior (zero coefficient of friction).

The steady state can be characterized equivalently by the rotation of the particle. The relevant dimensionless variable is then a = wjv. The steady-state diagram in the N -a space is shawn in Fig. 4. For a given value of the driving force , whatever the initial velocities, a is attracted to a point on the steady-state diagram. In phase 1, the steady-state value of a is equal to unity, while it is equal to zero in phase 3. ln phase 2, a decreases with the driving force as (11) It can also be shown that the particular motion of the particle, for a given driving force, is the one that minimizes the dissipation.

In the equations of the problem, the particular value of unity for the coefficients of friction appears to play a key role. So far, we have tacitly assumed that the two coefficients of friction were smaller than unity. However, it is now straightforward to understand what happens when this is not the case. If J.L < 1 and J.L 1 > 1, then the only steady state is phase 1: the particle can only roll on the plane. If J.L > 1 and J.L 1 < 1, then it is phase 3 that dominates the steady state: the particle can only slide on the plane. In the case where both J.L and J.L' are greater than unity, no motion at ali is possible! The particle simply sticks to the plane. It is easy to see that if the driving force N had a component downward, then the same situation would occur even for smaller values of the coefficients of friction. In this respect, unity has nothing special and reflects only the particular geometry we consider here. In a long array of particles this same effect takes place when the coefficients of friction are greater than unity, and the applied force is screened by the particle in contact with the pushing block. This observation suggests that in granular media, things will change radically when the interparticle coefficient of friction is high ènough. In the forthcoming sections we will assume that the coefficients of friction are smaller than unity.

III. A ONE-DIMENSIONAL ARRAY OF PARTICLES

The study of the motion of a single particle provides us already with sorne insight into the question we address in the general case of a long array. How are the rotations of the particles organized in the steady state? How does the motion of the system depend on the applied force? Our aim is to propose a global description of the frictional motion of this system as a result of the collective organization of rotations.

We have set up a computer program to simulate the motions of particles on a plane. The case of a single particle (Sec. II) shows that the friction law is actually the most important ingredient of the problem. Our program prescribes the exact Coulomb law (with no regularization) at the particle-particle and particle-plane contacts. We are interested in the situation where each particle remains in contact with its neighbors, so that the whole system can be considered as a single abject. This requirement is motivated by our purpose here to focus on the role played by the rotations in the translational motion of the array. Only if ali interparticle contacts are closed can the global translation be defined. Hence, the system we simulate is collision free and the relative velocities at contact points are only tangential. Moreover, our system involves no cohesion force and no elasticity.

As we shall see later in this section, the close-contact condition is satisfied only for a limited region in the parameter space. We add a second block to confine the array, where a constant force opposite to the direction of motion is applied to the first particle of the array. Although this new parameter is neither necessary nor suffident to prevent particle-particle contacts from opening, it is, in general, a very useful control parame ter, especially when the driving force is weak. So, the parameters are the following: J.L, J.L 1 , and J.L 11 are the coefficients of interparticle, plane-particle, and block-particle friction, respectively. The driving force NL, the confining force No, and the number of particles L are the other parameters, as shown in Fig. 5.

A. Governing equations Figure 5 shows a scheme of the array with the sign conventions for forces used throughout this paper. With these conventions, the equations of dynamics for each particle i are written

T(i)-T(i-1) + R(i) = 1, N(i)-N(i-1) + S(i) = T(i) + T(i-1)-S(i) v, Iw(i).
The boundary conditions are

N(L) = NL, N(O) =No. (12) (13)
There are 5L + 1 variables to be determined, whereas we have only 3L equations given by dynamics. The remaining 2L + 1 equations are prescribed by the contact law at the 2L + 1 contact points.

In order to avoid any artifact in the solution of our equations, we would like to use the exact Coulomb law of friction. However, it is not straightforward to implement it in an algorithm, since, as shown in Fig. 2, the (vr, T) should belong to a set of admissible values.

1.0 0.5

•8 0.0 -0.5 -l.O 0 '"'"' 1 <) --~2 -1-1-3 *----*4 0 5 .o • .o o ••• ~ i ,'\ 10 ' 1 1
' ' ~ ' .

15

. A mathematical framework for handling such problems, known as nonsmooth analysis, has been developed [START_REF] Jean | Proceedings of the Contact Mechanics International Symposium[END_REF]. In static equilibrium, where all the relative velocities are zero, Coulomb's law results indeed in an indeterminate state of forces (the so-called hyperstaticity of a sandpile at equilibrium). Nevertheless, in a dynamical system where the particles are moving and the states of contacts are permanently subject to change, the friction force on the vertical branch of Coulomb's law can be deter mined in a unique way. We have shown this explicitly in the case of the single-particle problem. The main point is that one has to distinguish the contacts where the relative tangential velocity remains equal to zero from those where this is not so. ln terms of the relative accelerations, we would say that there are two classes of contacts where the relative tangential velocity is zero: one where only the relative tangential velocity is equal to zero and one where both the relative tangential velocity and acceleration are equal to zero. We shall caU them active contacts and nonsliding contacts, respectively. Whenever a contact is active, i.e., when v,. = 0, there are three different alternatives distinguished by the relative tangential acceleration ù,. at the contact.

(1) v,. =O. In this case, the contact is nonsliding. The friction force T has to be in Coulomb's limits [-J..LN,J.LN],

where N is the normal force at the contact. Here we have an equaiton and an inequality. When the equation v,. = 0 supplements the equations of dynamics, we get the values of T and N and we can check for the inequality. If the inequality is satisfied as well, then we have the solution. If this is not so, we have to switch to one of the two other alternatives.

(2) v,. > O. ln this case, the contact is sliding although it is active. So, we have to put T = -J,LN following the sign convention of Fig. 2. We supplement the equations of dynamics with the latter equation, from which we calculate v,. among others, and we check for the inequality.

If the inequality is not satisfied, then we turn to another alternative.

(3) v,. < O. This is the same as the second case except that the corresponding equation is T = J,LN.

In this way, Coulomb's law in its most general formulation takes the following form: r = 0 and TE [-J.LN, J..LN], v,. = 0 => ~r > 0 and T = -J,LN, v,. < 0 and T = J,LN, (14

)
Vr > 0=> T= -J..LN, v,. < 0=> T=J..LN. ln order to solve the system of equations when a contact becomes active, the three alternative configurations are to be tried successively until the solution is found. Actually, there may be a great number of nonsliding contacts in a granular medium at the same time. Whenever a single contact becomes active, all the other nonsliding contacts are to be considered as weil. For p active contacts, this implies 3P possible configurations. Depending on the situation, different phenomena may take place. When a contact becomes active, it may become nonslid-ing, but cause other nonsliding contacts to turn sliding at the same time; or it may simply turn out to be a sliding contact and force other nonsliding contacts to become sliding ones as well, and so on. This basic image is too time consuming to be implemented in an algorithm, since the number of possible configurations grows exponentially with the number of active contacts. An important ingredient of our computer program is a relaxation scheme that converges rapidly to the solution in the space of configurations. The method used is to start with an arbitrary configuration when, in the process of evolution of the system, sorne contact becomes active. Then, if the selected configuration is not the right solution, it is mapped onto the space of configurations in ar der to build a new configuration to be tested. This iterative process converges very rapidly to the solution with the following mapping: v,. = 0 and T < -J.LN ----+ T = -J.LN, v,.= 0 and T > +J.LN ----+ T = +J.LN, (15) T =-J..LN and Vr < 0 ----+ v,. = 0, T = +J.LN and Vr > 0 ----+ v,.= o.

This procedure has to be applied simultaneously to all active contacts.

The evolution of the system can be studied starting from a random initial state. The easiest way to generate initial conditions is to choose initial velocities in such a way that ali contacts are sliding. For a given value of the driving force NL, all accelerations and forces can then be determined. These values do not change unless sorne contact becomes active. Meanwhile, the motion of each particle is uniformly accelerated and the time fJh needed for sorne contact to become active is simply given by . { . v,.(i) fJh =mm fJh(t) =-v,.(i)' oh(i) > o} . (16) At this moment, a relative tangential velocity at sorne contact vanishes and new values of forces and accelerations are to be calculated. Since the states of sorne contacts change, the values of forces and accelerations undergo discontinuous changes. The motion is again uniformly accelerated. This event-driven process continues until the system achieves the steady state where at every contact the relative velocity and acceleration are of the same sign or both are zero. In the steady state, all forces and accelerations rem ain the same forever.

B. Num erical res ults

Our simulations show that, depending on the driving force, three different regimes may occur. For small enough values of the driving force and for most initial conditions, sorne interparticle contact opens and the array separates in two or more identical arrays. On the other band, for large enough values of the driving force, sorne particle-plane contact may be lost and the geometry of the array is modified because of the rising of one particle. ln between these two limits, all contacts are preserved. We consider only this intermediate regime, ters are such that the steady-state value of Mis smaller than the particle-plane coefficient of friction.

The global steady-state coefficient of friction M of the array, although independent of the initial conditions, differs from the ordinary coefficient of friction in that it is a function of the driving force. Figure 9 shows this dependence for a system of 20 particles. M increases from zero linearly with N L and satura tes to a value equal to the particle-plane coefficient of friction. In other words, for weak values of N L, the driving force is accelerating the system less than it would if the particles were replaced by a single block. Another effect of the driving force is to regulate the internai organization of forces so as to increase its frictional resistance. For weak enough values of the driving force, the system is in phase 1. Dissipation takes place only at interparticle contacts. In this regime, M increases linearly with the driving force because, as we shall see in the next section, in phase 1 all forces scale with the driving force except for the normal force on the plane. As the driving force increases, more and more particles of phase 2 appear in the tail of phase 1 and hence the number of particles in phase 1 decreases. Then, Mis no more linear in NL. Finally, when the driving force is high enough, phase 1 disappears completely from the array so that ail particles are either in phase 2 or in phase 3. In this regime, M is simply equal to the particle-plane coefficient of friction. Thus, the number of rolling particles indicates the capacity of the system to increase its friction. Figure 10 shows the configuration of rotations in the steady state of a system of 40 particles for different values of the driving force. The connection between the rotation modes and the global behavior of the array in translation will be studied in more detail in the next section.

IV. THEORETICAL ANALYSIS OF THE STEADY STATE

In this section we will focus on the steady state and the mechanisms leading to the spatial patterns of rotations. These patterns provide the key to the global description of the array as a single object in displacement.

A. Different phases

The most interesting feature of the rotation patterns is that pure phases of rotation come one after another without mixing. This scheme suggests that as one goes along the array, sorne quantitative changes in the state of forces are followed by qualitative changes resulting in quite different behaviors. Let us consider first the case of rolling (phase 1).

Since ali of the particles in phase 1 are rolling in the same direction, friction forces at interparticle contacts are fully mobilized so that the equations of dynamics can be supplemented in this case by

T(i) = -J.LN(i). ( 18
)
Moreover, rolling implies

w(i) =v. ( 19
)
From Eqs. ( 12), (18), and (19) we get the following expression for the friction force at particle-plane contacts:

) ( ) i 2J.L 1 + I . 1 + J.L S(i) =v---(No + --v --, 1 + J.L 2J.L 1 -J.L (20)
as well as the following limit on the value of the particleplane force:

-S(i) < J. L' (1 + J.LV).

-

1 + J.LJ.L' (21)
Equation (20) shows that the particle-plane friction is mobilized progressively along the array and its direction is opposite to that of the motion. The absolute value of the friction force has an upper bound given by (21). From (20) and ( 21) one gets the maximum number of parti des in phase 1:

. {1 ( 1 1 -J.L ( . ')) t :::;

L 1 = 1 + n - , v + J.L 2J.L 1 -J.LJ.L ( 1+1 .)}/ (1+J.L) -ln N 0 + ~v ln 1 _ J.L , ( 22
)
where L 1 is the length of phase 1. So, the simple mechanism behind this phase is the progressive mobilization of the friction force at particle-plane contacts up to the maximum possible value of the friction force. Simulations in fact confirm this simple picture; see the functions S(i) and w(i) in Fig. 7 from the simulation of an array of 40 particles. The contacts of the particles coming after the particle i = L 1 with the plane can no longer be nonsliding and in phase 2, even if the particles continue rotating in the same direction, all contacts are sliding contacts. Now, let us consider the situation where all contacts are sliding contacts and all of the particles rotate in the positive direction (phase 2). The following equations are then to be used to supplement the equations of dynamics:

T(j) -11-N (j), S(j) = -11-' R(j) , w(j) > o. ( 23
)
From these and the equations of dynamics we get the following expression for the angular accelerations:

l . ( ') 211-(11-' +v) . !1-(1 + 11-') . W ) = - J+ V 1 -11-J-L' 1 -I-Ll-L' +11-'( 1 + 11-) -2 'N(. = 0) 1 1 Il- J , -11-/1- ( 24
)
where j refers to the jth particle of phase 2 and N(j = 0)

is the normal interparticle force on the fust partricle in the phase. At the same time, we have

S(j) = _11-'(1 + 11-v). 1 -11-11-' (25)
These equations show that the angular acceleration is decreasing linearly with the particle number while the friction force at the particle-plane contact remains constant. However, the positivity of the angular acceleration, which is a consistency condition, gives the maximum number of particles in the phase:

. <

L = _1_ [1 + 11-' v + 11-'(1 +11-) J - 2 11-' + v 2 2!1- -(1-J-L~-L')N(j = o)], (26)
where L2 is the length of phase 2. This simple mechanism behind phase 2 is illustrated in Fig. 7, where T(i) and w ( i ) are shown for a simulation of an arr a y of 40 partiel es.

From the end of phase 2, the behavior changes again radically as the particles can no longer rotate consistently in the same direction.

Let us now consider phase 3 , where the interparticle contacts can no longer be sliding. The equations to be used to supplement the equations of the dynamics are the following:

S(k) = -11-' R(k), w(k) = -w(k + 1) =w. (27)
Assuming that the fust particle in this phase in contact with the last particle of phase 2 bas a negative rotation velocity, we immediately get the following expression for the tangential interparticle force:

T(k) = T(k = 0) + -+ - [ 11-' lw] 2 2j.L 1 ( ~)k 1 -1-L' k lw 11-' -(-1) ---. 211-' 2 (28)
The interparticle friction force oscillates with a period of two particles and grows exponentially. We remark that this behavior of the interparticle friction forc e is only related to the fact that the interparticle contacts here are nonsliding. Even in the case where w = 0, which is a particular case of phase 3, this oscillating behavior is obtained. ln the same way, we get the following expression for the normal interparticle force:

N(k) = N(k = 0) + k(v + 11-') [ 11-' lw] + 11-' T(k = 0) + 2 + 2 11-, [ ( 1 + ') kl x 1--1-:1 . ( 29
)
From Eqs. ( 28) and ( 29) and the boundary condition

T(O) = -11-N(O) , it is easy to see that JT(k) l $ 11-N(k)
for all k. The length of phase 3 is thus not limited by the mobilization of the interparticle friction forces. However, sorne limit arises from the normal particle-plane force which, for consistency of the equations, bas to be positive. The expression of the normal particle-plane force is the following:

R(k) = 1 + (-1)k 1 ~ 11- +--T(k = 0) + !:!:__ + __:::._ 2 [ 1 l'] 1 -11-' 2 211-' x [-1 + 11-']k-1 1 -11-' (30)
The positivity of R(k) sets an upper limit on the length of phase 3:

k $ L3 = 1 + {ln [ 1 ~ 11-' ( 1 + :, ) ] [ 1-L' lw ] } 1 ( 1 + 11-' ) -ln T(O) + 2 + 2 11-, ln 1 _ 11-' , (31)
where L 3 is the maximum length of the phase and T(O)

is to be replaced with the last particle in phase 2. The particle coming just after the particle L 1 + L 2 + L3 will rise as the reaction force on the plane vanishes.

What is the relation between w and v in phase 3?

There are not enough equations to relate these two accelerations. The answer is to be looked for in the boundary conditions. For the array to be stable (no contact opening), we should have

L $ L 1 + L 2 + L 3 .
Suppose that there is an even number of particles in phase 3 so that the last particle is rotating in the positive direction. As the latter is in contact with a black, with a coefficient of friction generally different from J.L, its angular acceleration is different from that of phase 3. Simulations show that if J, . L 11 < J.L, then the angular acceleration of the last particle is greater than that of phase 3. In this case, the contact between the last particle of the array and the last particle of phase 3 is sliding and we have the consistency condition T(k = L~-1) = -J,.LN(k = L~-1), where L'

is the length of phase 3. This gives the relation between w and v. As seen in Fig. 10, w decreases with the driving force, while v increases. Figure 11 shows the angular accelerations of the particles in two simulations: one with J.L 11 = 0 , the other with J.L 11 > J.L. ln this latter case, the rotations in phase 3 simply vanish.

B. Global dynamics

The global behavior of the array results from the juxtaposition of the three rotation phases. By global behavior, we mean the relation between the applied force NL -No , as input to the system, and the acceleration v, which is a global output. ln this section we will show that, as in the case of a single particle (Sec. II) , this relation can be cast in the canonical form of the equation of frictional motion of a single black:

NL-No-F* = m*v , ( 32
)
where m* , the effective inertia and F*, the effective friction force on the plane, are functions of the parameters and the applied force. Figure 12 shows v as a function of N L -No for an array of 20 particles. The inverse of the slope at each point corresponds to the effective inertia for the corresponding value of the applied force, while the intersection with the axis of accelerations gives the effective friction force on the array. Three different regimes can be distinguished.

(1) From zero up to a critical value N' of the applied force, the effective mass is constant. In this regime ali particles are in phase 1. From Eqs. (12) , {18), and (19) it can be shawn that the normal interparticle force in phase 1 is given by

N( .) _ (N 1 + I . ) ( 1 + J.L) i 1 + l . ' 1. - o + --v -----v 2J,.L 1 -J.L 2J,.L ' ( 33
)
so that the effective mass is equal to

m • = ~ [ ( 1 + J.L) L -1] 1 2J,.L 1 -J.L (34)
Thus, the effective inertia in this regime can be rouch greater than the real mass L of the system. lt increases rapidly with J.L and for small values of J.L , is at least 1 + I times the real mass of the array. From Eq. (33) we get also the following expression for the effective friction force:

(35) ln Sec. III, we introduced a global coefficient of friction M, which is the "total'' friction force F on the plane divided by the real mass L of the system. Since a global description implies analytically an effective inertia to be introduced, a description in terms of "effective" quantities is more consistent. Hence, we introduce the effective renormalized coefficient of friction by

M* = F*. m* (36)
We recall that the acceleration of gravity bas been set to 1. From Eqs. (34) and {35) , the expression of M* in the regime of collective rolling is given by (37) This is a simple expression in that it depends neither on the applied force nor on the number of particles. In this respect, M{ shows exactly ali the features of a Coulombian coefficient of friction. The confining force No in this regime is a control parameter of M{. When No = 0, the effective coefficient of friction vanishes! The system is bence reduced to a black with no friction on the plane and an effective mass greater than its real mass. This has interesting physical implications. For instance, since the effective inertia is greater than the real mass, the speed of sound through the system, due to a small elasticity at contact points, should be less by at least a factor of v'1 + I than when rotations are absent. On the other hand, since the friction force for small values of 1-L is proportional to the number of particles L, a small change in the confining force is amplified by a factor L to give a much larger friction force.

(2) When the applied force is greater than N', the effect ive mass begins to increase with the applied force. At the same time, the number of particles in phase 1 decreases and there are more and more particles in phases 2 and 3. As long as there are particles in phase 1, the effective mass continues to increase. From Eqs. ( 12) and ( 23) , it can be shown that the effective mass and coefficient of friction in phase 2 are given by m• -

1 L 2 - 1-J.<p.' Mi = 1-L' ( 38
)
where Lis the number of particles in phase 2. In the same way, for phase 3, from Eqs. ( 12) and ( 27), we obtain

m3 = L , M* = u' -~ [(~)L -1] 3 ,_ L 1 -p.' x ( -~.LN(k = 0) + ~ + J;,). (39)
where Lis the length of phase 3. Here we have taken w to be an independent parameter which is related essentially on the particle-block coefficient of friction. The effective coefficient of friction in phase 3 is approximately equal to f.L 1 • The deviation from this value is of the order of I.L' 3 when the angular acceleration is zero, and of the order of As the number of particles in phase 1 decreases, the effective coefficient of friction increases and, at the same time, the effective inertia of the array tends to its real mass.

(3) For a certain value N" of the applied force, there will be no particle in phase 1 and the relation between acceleration and applied force becomes linear again. In this regime, the effective mass is practically equal to the real mass of the system, except for small oscillations coming from phase 3. The effective coefficient of friction in this regime is equal to the particle-plane coefficient of friction. Figures 13 and 14 display the variation of effective coefficient of friction and effective inertia as a function of the applied force for an array of 20 particles. N 0 and f.L are the control parameters of the first regime.

It can be seen from Fig. 13 that the necessary condition for observing the fust regime for any number of particles is given by a simple inequality: When the confining force is small, the system is always in the fust regime for small enough values of t he applied force. Collective rolling of particles should occur for example at the free surface of sandpiles. Therefore, physical effects due to this mode should be observable in real sandpiles.

V . CONCLU SION

We have introduced a simple one-dimensional model of particles where the interplay of the dynamics and the contact law (Coulomb's law of fric tion) leads to a well defined organization of the rotations of particles in the steady state. This self-organization has been studied numerically and analyzed in terms of a juxtaposition of "pure modes." The latter have been fully characterized analytically and their respective lengths have been computed. The agreement between the numerical simulations and the theoretical analysis is perfect. The translational motion of the array on the plane in the steady state can be described in terms of a global friction force, increasing with the external driving force. However, direct analysis shows that an effective inertia has to be introduced. This in turn requires an "effective" coefficient of friction to be defined, which is different from the "global" coefficient of friction in many respects.

Finally, let us stress the importance of length scales intermediate between the particle size and the system size appearing in rotation modes and force patterns, which may play a significant role in a continuum description of granular media. While this effect is limited here to

  FIG.3. Evolution of the particle-plane friction force as a function of the driving force N. Ali forces are normalized with respect to the weight of one particle (see text) . The

NFIG . 4 .

 4 FIG .4. Steady state of the single-particle problem inN -a space, where a = % • Ail forces and accelerations are normalized ( see text). The steady state depends only on parameters and not on initial conditions. Here, J.L = 0.1 and J.L 1 = 0.2.

FIG. 5 .

 5 FIG. 5. {a) An array of particles on a plane. (b) Forces exerted on the particle i and sign conventions.

.FIG. 10 .

 10 FIG. 9. Variation of the global coefficient of friction of an array of 20 particles as a function of the driving force. Parameters are J.L = 0.01, J.L 1 = 0.1 , and No = 0.05.

  FIG. 11. Angular accelerations of particles in an array of ten particles for two different values of the particle-block coefficient of friction J.l-11 • The other parameters are JI-= 0.1 ,

FIG. 12 .

 12 FIG.12. Translational acceleration of an array of 20 partides as a function of the applied force N L -N 0 . Parameters are J.L = 0.01, J.L 1 = 0.1, and No= 0.05.

  The effective inertia and coefficient of friction of the array can be calculated from those of the phases:m • mi + m~ + mj , M * = Mimi+Mâmi + M;m;m.i +m;+m;(40)

  FIG.13. Variation of the effective coefficient of friction as a function of the applied force for the system of Fig.12. Mt = 2J.t/ (l + I)No is the effective coefficient of friction in phase 1. M3 = J.t 1 is the mean value of the coefficient of friction in phases 2 and 3.

FIG

  FIG.14. Variation of the effective massas a function of the applied force for the system of Fig.12.

ACKNOWLEDGMENTS

It is a pleasure to acknowledge useful discussions with D. Bideau, H. J. Herrmann, J. P. Troadec, and D. Wolf. This work was supported by the Groupement de Recherche "Physique des Milieux Hétérogènes Complexes" of the CNRS.


Student Member, IEEE Andres Mauricio Lopez

Hubert Piquet

Member, IEEE Diego Patino

Member, IEEE Rafael Diez

Xavier Bonnin

Parameters Identification and Gas Behavior Characterization of DBD Systems

Keywords: Dielectric barrier discharge (DBD), gas conductance, gray box, modeling, parameter identification

This paper proposes an efficient modeling and an identification method for dielectric barrier discharge (DBD) systems, based on input-output (current-voltage) experimental measurements. The DBD is modeled using an equivalent electric circuit associated with a differential equation that describes the dynamics of its conductance. This equation assumes a homogeneous behavior of the gas. This paper introduces a series of polynomial terms of the current of the gas into the conductance equation. These terms, after identification, are a very useful tool to analyze the physical mechanisms that take place in the gas. The identification process also returns the numerical values of other DBD parameters, such as associated capacitances and the breakdown voltage. In addition, an asymmetric model for the gas, which considers the direction of the current, is proposed to consider the possible geometrical dissimilarity between the two electrodes of the DBD setup. Experimental measurements taken on two different DBD applications are used for validating the proposed approach.

I. INTRODUCTION

D IELECTRIC barrier discharges (DBDs) offer specific properties that are very interesting for a large number of processes. The obtained plasma is maintained at low temperatures, the barriers contribute to spread the discharge on the whole surface of the electrodes and they avoid direct contact between the metallic electrodes and the gas. Thus, it increases lifetimes and reduces costs. For these reasons, DBD are widely used in industrial applications such as ozone generators, surface treatments, thin film deposits and coating, lighting, and displays [START_REF] Kogelschatz | Dielectric-barrier discharges: Their history, discharge physics, and industrial applications[END_REF]- [START_REF] Wagner | The barrier discharge: Basic properties and applications to surface treatment[END_REF].

Barriers used in DBD setups are made of dielectric materials (ceramics, glass for instance), which forbid the use of a DC power supply to inject power into the system. Therefore, the electric power source must be especially designed for A. M. Lopez is with Pontificia Universidad Javeriana, Bogota 1100100, Colombia (e-mail: lopez.andres@javeriana.edu.co).

H. Piquet and X. Bonnin are with the LAPLACE laboratory, Toulouse 31071, France (e-mail: hubert.piquet@laplace.univ-tlse.fr; bonnin@ laplace.univ-tlse.fr). R. Diez and D. Patino are with the Department of Electronics, Pontificia Universidad Javeriana, Bogota 12345, Colombia (e-mail: rdiez@ javeriana.edu.co; patino-d@javeriana.edu.co). Fig. 1. Detailed electrical model of a DBD. V s and I s represent the generator voltage and current, C gas and C d the capacitance of the gas and the equivalent capacitance of the dielectric, V gas and I gas the voltage and current of the gas, and G the conductance of the gas.

each DBD process if precise operating points and effective control of the power transfers are desired. To make this possible, an electric model of the DBD, relating the current and voltage at the electrodes, is very useful to correctly dimension the components of the power converters at the design stage.

In this scope, the availability of such a model permits the achievement of numerical simulations (with software such as PSIM or PSPICE) associating the circuit of the electrical generator to the DBD system which has to supply. These simulations, as described in [START_REF] Piquet | DBD lamp converter design using an electrical model of the load[END_REF], allow the prediction of the operating conditions for the whole system. It also gives the possibility to know the power that is transferred to the gas and therefore, to the process in the plasma. In addition, the waveforms and magnitudes of the electrical variables of the electronic components in the power supply can be known as well.

The current-voltage relationship that describes the behavior of the DBD device can be found, using experimental waveforms with various approaches. Among them, the black-box technique is able to define an accurate model of the system. Nevertheless, it is meaningless from the physics and electrical point of view [START_REF] Sjöberg | Nonlinear black-box modeling in system identification: A unified overview[END_REF]- [START_REF] Castano | A fast approximation algorithm for set-membership system identification[END_REF]. On the other hand, the electrical or circuit-based model, displays a capacitance for the dielectrics, a capacitance for the gas, and a conductance for the gas. This model gives more insight into the physics and eases the implementation in simulation software for electric circuits. In this paper, the second approach is used; a nonlinear differential equation to model the dynamics of the gas conductance is proposed. Therefore, nonlinear identification method is used to determine the parameters of this model.

The equivalent circuit of the DBD used in this paper is derived from the well-known model mentioned in [START_REF] Diez | Predictive model of a DBD lamp for power supply design and method for the automatic identification of its parameters[END_REF] (see Fig. 1). This model includes the capacitances of the dielectric barriers, the capacitance of the gas gap, and the conductance of the gas.

One of the most relevant contributions of this paper is obtained because of the introduction of polynomial series terms of the gas current into the equation governing the evolution of the conductance. These terms provide crucial information concerning the behavior of the discharge in the gas gap. They may vary depending on the operating condition, which are not always a priori known. Those terms also help in making the balance between the most significant phenomena, which can take place in the plasma, e.g. thermal influence, ionization, breakdown, and recombination of the excited species [START_REF] Liu | A PSPICE model for the electrical characteristics of fluorescent lamps[END_REF], [START_REF] Thomas | The simulation of circuit breaker switching using a composite cassie-modified Mayr model[END_REF].

The proposed approach is tested on two completely different DBD applications: 1) an experimental reduced scale setup for plastic film surface treatments and 2) a 90 W excimer UV lamp. The first application is a coplanar DBD in nitrogen at atmospheric pressure, oriented to film-coating, as shown in Fig. 2 (left). Both the electrodes are covered with a ceramic layer, acting as dielectric barrier. The second application concerns excimer lamps, which are high performance UV radiation generators [START_REF] Tarasenko | High power VUV and UV excilamps[END_REF]- [START_REF] Lomaev | Capacitive and barrier discharge excilamps and their applications (review)[END_REF]. This UV is produced, in the gas mixture contained in the lamp during the decay of excited molecules returning to their ground state after an electric excitation. In the selected bulb, the DBD is produced in a Kr-Cl gas mixture at low pressure. Hence, the setup presents an asymmetric configuration because the outer electrode is larger than the inner electrode [see Fig. 2 (right)], with different radius of curvature and thus different local electric fields.

Considering this geometry, this paper also deals with the possible dissimilarities between the electrodes on a DBD setup. Expecting that DBD setups with different sizes or shapes of electrodes, can exhibit a specific behavior of the gas depending on the current direction.

This paper is organized as follows: Section II enounces the electrical model of the DBD and its mathematical representation of the gas behavior for the symmetrical and asymmetrical cases. Section III presents the identification technique. Finally Section IV validates the model and the identification method for the two proposed applications, with a detailed analysis of the terms in the conductance equation.

II. DBD ELECTRICAL MODEL

A well-known equivalent circuit of the DBD setup consist of one or two dielectric barriers and a gas volume as shown in Fig. 1 (C d represents the equivalent series capacitance of the dielectric barriers). C gas is the capacitance of the gas when no electric discharge takes place in the latter. G is the conductance of the gas when a plasma is installed. One of the main assumption and limitation of this modeling approach is that a homogeneous behavior of the gas is assumed.

Parameters such as the capacitances of the dielectrics, capacitance of the gas, and constants related to the gas are needed to implement the equivalent model in a circuit simulator. Therefore, they have to be estimated on the basis of the measurable variables. In an experimental setup, the only measurable variables are the current through the DBD I s and its voltage V s .

The structure of this equivalent model of the DBD is similar to the one proposed in [START_REF] Diez | Predictive model of a DBD lamp for power supply design and method for the automatic identification of its parameters[END_REF]. According to the circuit of Fig. 1, the system is described by the following differential state equations:

dV gas dt = I s -V gas G C gas ( 1
)
dV C d dt = I s C d . ( 2
)
The voltage of the DBD setup is easily computed as follows:

V s = V C d + V gas .
(3)

A. Symmetrical Model

For the case of a symmetrical DBD setup such as a plastic films surface treatment reactor, the equation of the conductance G is proposed to be defined with several aspects of the behavior of the electric discharge. Polynomial terms of the current of the gas are added as follows:

dG dt = K 1 u(V th -|V gas |) -K 2 G + ∞ n=1 T n |I gas | n (4)
where the function |.| is the absolute value of the argument and u(.) the Heaviside function and V th is the breakdown voltage of the gas. It can be observed that: 1) the first term represents the ionization process at the breakdown and K 1 is a generation coefficient; 2) the second term represents the recombination process after breakdown and K 2 is the extinction coefficient; 3) In the last term, T n coefficients reflect the relationship between the conductance and different powers of the conduction current into the gas I gas . The relative importance of these terms reflects the properties of the gas and the physical phenomena that take place in the latter.

B. Asymmetric Model

Some DBD systems, such as DBD lamps have an asymmetric structure, with two coaxial cylinders. Therefore, the inner and outer electrodes also have different radii. Depending on the direction of the current of the lamp, the radius of the anode is bigger or smaller than the radius of the cathode (see Fig. 3) and the magnitude of the local electrical field in the neighborhood of the cathode is modified.

To consider this asymmetry in the mathematical model and to improve the identification results, the equation of the conductance of the gas is proposed to be modified as follows:

dG dt = K 1 u(V th -|V gas |)-K 2 G + ∞ n=1 T n I n gas , if I gas ≥ 0 K † 1 u(V † th -|V gas |)-K † 2 G + ∞ n=1 T † n (-I gas ) n , if I gas < 0. ( 5
)
Different parameters are used for each direction of the current. The values of the dielectric capacitance C d in (1) and the gas capacitance C gas in (2) are assumed to be independent on the current direction.

III. IDENTIFICATION OF THE MODEL PARAMETERS

To obtain the electrical and physical constants of the DBD system, a gray-box identification technique is considered. This method maintains the physical and electrical meaning for the

A. Gray-Box Identification Algorithm

As seen in Section II, the system is modeled as a set of differential equations with a set of unknown parameters p to be identified, as follows:

ẋ = f (x, p, u), y = h(x, p, u) ( 6
)
where x is the state space vector, p is the vector of unknown parameters to be identified, u the inputs and y the outputs. For this case u = I s and y = V s . The functions f (x, p, u) and h(x, p, u) are known: f is described in (1)-( 4) and h is computed according to (3). The inputs u and outputs y are measurable and known. The objective of the identification process is to minimize the norm of the quadratic error, varying the parameters p as follows:

min p ( e( p) 2 ) = min p ( y -ŷ 2 ) (7)
where * 2 is the square of the norm of the argument, e( p) is the error between ŷ the measured output and y the estimated one. In a general case, the error is a nonlinear function of p [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF], [START_REF] More | The Levenberg-Marquardt algorithm: Implementation and theory[END_REF].

The algorithm approximates the error function around p 0 (initial condition of the vector p) with a simpler quadratic function, then the minimization problem ( 7) is reduced to the following:

min s 1 2 s T H s + s T g, s < . ( 8
)
A trial step s is computed to try to minimize the quadratic approximation of the error in [START_REF] Diez | Predictive model of a DBD lamp for power supply design and method for the automatic identification of its parameters[END_REF]. Then, it is verified if e( p) is reduced. In this case, the parameter p 0 is updated as p new 0 = Fig. 5. Response of the identified model (FIT = 90.5%) to the identification waveform. p 0 + s if e( p 0 + s) < e( p 0 ). Otherwise, the parameter is reduced [START_REF] Coleman | An interior trust region approach for nonlinear minimization subject to bounds[END_REF].

This algorithm is also capable to deal with constrains in the parameters p (see [START_REF] Coleman | An interior trust region approach for nonlinear minimization subject to bounds[END_REF] and [START_REF] Coleman | On the convergence of interior-reflective Newton methods for nonlinear minimization subject to bounds[END_REF] for further information).

B. Summary of DBD Setup Identification Problem

For the DBD system, the set of nonlinear differential equations ( 6) is defined in (1)-( 4).

The components of the state space vector x are V d , V gas , and G. The input variable u is the DBD current I s and the output y is the DBD voltage V s .

The set of parameters to be identified are p = [C d , C gas ,

V th , K 1 , K 2 , T n with n = 1, 2, 3 . . . , ∞].

IV. IDENTIFICATION RESULTS

The identification method is used for two different DBD systems, considered in this section: a DBD for film treatment with a planar structure and a DBD UV lamp with a coaxial structure (see Fig. 2).

A. DBD for Film Treatment

Fig. 4 shows the measured waveforms used for the identification, current I s and voltage V s in the DBD system. The specific shape of the current has been obtained with an especially designed current-controlled power supply [START_REF] Bonnin | Alimentation électrique des dispositifs à décharge à barrière diélectrique (DBD)[END_REF]. Initial guess of parameters set is deduced from the characteristics of the DBD shown in Table I.

1) Identification Results: The algorithm needs 90 iterations to converge. Table II shows the obtained parameters.

Fig. 5 shows a comparison of the actual response of the DBD (V s voltage) and the response obtained with the identified model. The FIT index is used to measure the quality of the identified model. It is calculated as follows:

FIT = 100 1 - ŷ -y y -ȳ % ( 9 )
Fig. 6. Estimated current I gas and conductance G of a DBD system for film treatment.

Fig. 7. Behavior of the gas (I gas versus V gas ) of a DBD system for film treatment.

where * is the norm of the argument, y the results obtained with the theoretical model, ŷ the experimental measurements, and ȳ the average of the vector ŷ. The model clearly reflects the behavior of the system for the identification data with FIT = 90.5%. One of the advantages of the equivalent circuit approach is the possibility to estimate electrical variables such as the current in the gas I gas or the conductance G (Fig. 6). It is possible to see that the value of the gas conductance G is stable as long as the current flows through the gas.

Fig. 7 shows the behavior of the gas current I gas versus the gas voltage V gas . When the gas is ignited, its voltage remains almost constant for different values of current, suggesting a normal glow operating regime.

2) Comparison of Polynomial Terms:

To simplify the identified model, the relative importance of each polynomial term (i.e, T n |I gas | n ) is analyzed. Fig. 8 shows a comparison of the polynomial terms for the film treatment DBD setup. For the comparison, the RMS which represents the average power of Fig. 8. Comparison of RMS value of the polynomial terms of the current of the gas I gas for the identification waveforms. RMS values of terms with fourth power and greater are in magnitude order ×10 -9 and are not shown. the waveform is used. It is computed as follows:

T n |I gas | n RMS = 1 T T 0 (T n |I gas (τ )| n ) 2 dτ ( 10
)
where T is an integer multiple of the operation period. For this DBD the most important term is T 1 |I gas |. The RMS of this term is more than 10 times greater than for T 2 |I gas | 2 . Fig. 8 confirms that the gas conductances can be sustained by means of the gas current. For instance, let us consider only the first term of the polynomial powers in the gas current T 1 |I g | as important. If a slowly varying current is applied during certain time after breakdown, (4) states that the conductance is proportional to the absolute value of the current amplitude in steady state (as seen in Fig. 6). This confirms the normal glow regime.

3) Predictability of the Model:

The identified model is validated with five other waveforms, with different frequencies and amplitudes, obtaining an average FIT of 85.3%. Fig. 9 shows the waveform (measured and estimated) with the worst FIT.

B. DBD Lamp

The identification algorithm is used for the characterization of a DBD Kr-Cl excimer UV lamp. The waveforms of the measured current I s (input) and voltage of the lamp V s (output) are shown in Fig. 10. These figures are obtained with a specifically designed power supply, designed for the research of the optimal operating conditions of the lamp [START_REF] Florez | DBD excimer lamp power supply with fully controlled operating conditions[END_REF]. This power supply is capable of changing the frequency f , amplitude I , and duty cycle d of the current waveform (see Fig. 10).

1) Identification Results (Symmetrical Model):

To compare the results obtained with symmetrical and asymmetrical model in the DBD lamp, the identification algorithm is run at first for the symmetrical case. The identification algorithm converges in 80 iterations. Table III shows the numerical results of the parameters of the model. Fig. 11 shows a comparison between the measured DBD lamp voltage and the response obtained with the identified model. The FIT index is used to measure the quality of the identified model as in the case of the DBD for film treatment. The model clearly reflects the behavior of the DBD lamp for the identification data with FIT = 91.7%. Response of the identified models (FIT = 91.7%) to the identification waveform. Fig. 12. Comparison of RMS value of the polynomial terms of the current of the gas I gas for the identification waveforms using the polynomial model. RMS values of terms with seventh power and greater are in magnitude order ×10 -9 and are not shown.

2) Comparison of Polynomial Terms:

To simplify the model, the relative importance of each of the polynomial terms (i.e, T n |I gas |) is analyzed. Fig. 12 shows a comparison of the polynomial terms for the identification waveform. For this comparison the RMS value ( 10) is used. It can be seen that in this case the two first powers T 1 |I gas | and T 2 I 2 gas are the most important indicating that thermal processes are involved in the discharge and approaching to the Mayr model used for low current arcs [START_REF] Liu | A PSPICE model for the electrical characteristics of fluorescent lamps[END_REF], [START_REF] Thomas | The simulation of circuit breaker switching using a composite cassie-modified Mayr model[END_REF].

3) Predictability of the Model: To test the accuracy of the identified model for different waveforms a validation process is carried out. Waveforms with different frequencies and duty cycles are used. Fig. 14 shows the obtained FIT. It can be seen that a lowest FIT is obtained for lower values of duty cycle. However, the lower value FIT = 80% still shows a satisfactory response of the model. Fig. 13 shows the response of the model to a different current (I s ) waveform generated by a resonant converter. For this case the FIT is 90.1%.

4) Identification Results (Asymmetrical Model):

Table IV shows the results of the algorithm after 90 iterations (Fig. 10). The results show that for the asymmetric model the parameters vary depending on the current direction. In addition, the parameters that are independent on the conductance equation C gas and C d vary less than 5% with respect to the results of Table III.

Fig. 15 shows a comparison of the measured response of the DBD lamp and the response obtained with the asymmetric and nonasymmetric models. Again, the FIT index is used to measure the quality of the identified model. Concerning the predictability of the asymmetric model, Fig. 16 shows the values of the FIT according to various operating conditions. It can be seen that the lowest value of the FIT is 85%. This is an improvement with respect to the results of Fig. 14 using the nonasymmetric model, where the lowest FIT is 80.5%.

Similar to the coplanar DBD (surface treatment process), with the numerical parameters of the model and the mathematical equations, internal variables such as gas current I gas , gas voltage V gas , voltage of the dielectric V d , and gas conductance G can be estimated (see Fig. 17).

V. CONCLUSION

In this paper, a nonlinear identification technique for DBD systems is proposed and evaluated. This technique allows obtaining physical parameters of the lamp based on experimental measurements of current and voltage. In addition, the identified model allows estimating unmeasurable quantities such as the current of the gas.

The identification technique is tested on a DBD for film treatment and on a DBD UV lamp, obtaining very satisfactory results. On this basis the importance of each of the polynomial terms (powers of the current of the gas used to describe the gas conductance) is compared, and it is possible to carry on an analysis of the prominent phenomena, which occur in the plasma. In addition, a model that considers the asymmetry of some DBD systems such as DBD lamps is proposed and tested. This model shows an improvement in the validation results with respect to the original symmetric model.

Fig. 2 .

 2 Fig. 2. Physical structures of the considered DBD setups.

Fig. 3 .

 3 Fig. 3. Lateral view of the lamp with current directions.

Fig. 4 .

 4 Fig. 4. Measured current (input) and voltage (output) of the DBD for film treatment.TABLEI

Fig. 9 .

 9 Fig. 9. Measured current I s and measured and estimated voltage V s of the DBD for film treatment FIT = 81.1%.

Fig. 10 .

 10 Fig. 10. Measurements of the DBD lamp current I s (input) and voltage V s (output) used for identification.

Fig. 11 .

 11 Fig. 11.

Fig. 13 .

 13 Fig. 13. Measured current I s and measured and estimated voltage V s of the DBD lamp using a resonant converter FIT = 90.1%.

Fig. 14 .Fig. 15 .

 1415 Fig. 14. FIT index for validation waveforms using the polynomial model.

Fig. 16 .

 16 Fig.[START_REF] More | The Levenberg-Marquardt algorithm: Implementation and theory[END_REF]. FIT index for validation waveforms using the asymmetric model.

Fig. 17 .

 17 Fig. 17. Estimated variables of the DBD lamp using the identified model.

TABLE I CHARACTERISTICS

 I OF THE DBD FOR FILM TREATMENT

TABLE II NUMERICAL

 II RESULTS OF THE PARAMETERS OF THE DBD FOR FILMTREATMENT UP TO T 10 identified parameters. It requires previous knowledge of the dynamical system behavior.

TABLE III NUMERICAL

 III RESULTS OF THE PARAMETERS OF THE IDENTIFIED MODELSFOR THE DBD LAMP UP TO T 10

TABLE IV NUMERICAL

 IV RESULTS OF THE PARAMETERS OF THE ASYMMETRIC MODEL FOR THE DBD LAMP UP TO T 10

He is currently a Full Professor with ENSEEIHT Ing. School, Toulouse. He teaches power electronics and systemic approach in electrical engineering. His research activity takes place in the LAPLACE laboratory, Energy Conversion and Plasma Laboratory, Toulouse. His current research interests include quality and stability in embedded networks as well as power supplies for plasma applications.

Diego Patino (M'00) was born in Manizales, Colombia. He received the B.Sc. degree in electronic engineering from Universidad Nacional de Colombia, Bogotá, Colombia, the M.Sc. degree from Universidad de Los Andes, Bogotá, and the Ph.D. degree in automatic control from Nancy University, Nancy, France.

He is currently a Full-Time Professor with Pontificia Universidad Javeriana, Bogotá. He is the author of several articles. His current research interests include hybrid dynamical systems, power converter, and nonlinear control theory. His work concerns both the study of the coupling between electrical power sources and plasma devices and the methodologies in designing power electronics components for plasma applications.

Rafael