Benjamin Dollet

Claire Bocher

Flow of foam through a convergent channel

Keywords: XX, XX, XX No code given

We study experimentally the flow of a foam confined as a bubble monolayer between two plates through a convergent channel. We quantify the velocity, the distribution and orientation of plastic events, and the elastic stress, using image analysis. We use two different soap solutions: a sodium dodecyl sulfate (SDS) solution, with a negligible wall friction between the bubbles and the confining plates, and a mixture containing a fatty acid, giving a large wall friction. We show that for SDS solutions, the velocity profile obeys a self-similar form which results from the superposition of plastic events, and the elastic deformation is uniform. For the other solution, the velocity field differs and the elastic deformation increases towards the exit of the channel. We discuss and quantify the role of wall friction on the velocity profile, the elastic deformation, and the rate of plastic events.

Introduction

Aqueous foams are metastable dispersions of gas bubbles within a continuous liquid phase, which is a solution of surfactants, which stabilise the liquid films separating adjacent bubbles from rupture [START_REF] Cantat | Foams, Structure and Dynamics[END_REF]. Two physical parameters are crucial in the description of foams: the liquid fraction ϕ ℓ , and the mean bubble area A (or, equivalently, the mean bubble radius).

By opposition to bubbly liquids, foams have generally a liquid fraction low enough to be a jammed system: the bubbles are packed and share thin films. Hence, bubbles have to deform elastically by a substantial amount before flowing one past each other, and foams then acquire elastoplastic properties. They thus belong to the vast class of complex fluids, like polymers, pastes or slurries. Contrary to most other complex fluids, the microstructural item of foams, i.e. the bubble, is easy to image and its elementary micromechanics is defined unambiguously. More precisely, the elastic stress can be directly deduced from the geometry of the bubble boundaries [START_REF] Batchelor | [END_REF], and there is a clear elementary plastic event, the so-called T1 event: it is the neighbour-swapping of four adjacent bubbles. The macroscopic fields (elastic stress, plastic deformation, and velocity) describing the mechanics of a foam as a whole then directly result from averaging the information obtained at the bubble scale. Therefore, foams (as well as dense, jammed emulsions) are model systems to understand the flow of complex fluids across the scales.

Because foams are opaque in bulk, obtaining the information at the bubble scale is difficult on three-dimensional (3D) foams and the local description of 3D foam flows has been performed only very recently, using fast X-ray tomography [3]. Most experimental studies have resorted to foams confined in two dimensions (2D), e.g. between two glass plates, which renders the structure of the bubbles transparent (Fig. 1b). For these confined systems, various flow geometries have been studied: parallel flows driven by the motion of side walls (Couette flows) [4,5] or by a pressure gradient (Poiseuille flows) [6], or flows in more complex geometries such as past obstacles [7][8][9] or through constrictions [10,11]. Similarly, most numerical simulations of foams in pure shear [12][13][14], past obstacles [15][16][17][18], or through constrictions [19,20] have been performed in 2D, with some exceptions in 3D, limited to monodisperse foams in quasistatic regimes [21,22].

A specificity of foams confined in 2D by plates is the presence of friction (which we shall henceforth call "wall friction") between the moving bubbles and the fixed walls, which increases at increasing foam velocity [23]. If the velocity is low enough, the effect of wall friction remains negligible; in particular, the foam flow is quasistatic. On the other hand, at large velocity, wall friction competes with the internal stress of the foam, of elastic and viscous origin, and can even lead to the catastrophic failure of the foam, reminiscent of brittle fracture [24][25][26].

So far, the flow of foams or dense emulsions in a convergent channel has not been much studied, except by the groups of Earnshaw and Weeks, with a special emphasis on the local description of plastic events [27][28][29]. However, this configuration is particularly interesting, because the convergent flow is complex enough to include both shear and elongation, contrary to parallel flows, but simple enough to be amenable to some analytical predictions. Indeed, in Newtonian fluid mechanics, it is one of the very few configurations where the full Navier-Stokes equations admit an analytical solution, at least for a steady and purely radial flow [30,31], with in particular an explicit expression of the velocity profile in the boundary layer at high Reynolds number [32,[START_REF] Landau | Mécanique des fluides[END_REF]. This comes from the fact that for incompressible flows at velocity v, the symmetry of the convergent channel (neglecting end effects) giving a radial flow: v = v(r, θ)e r (see Fig. 1b for the definition of the polar coordinates), combined with the continuity equation ∇ • v = 0, yields a particularly simple expression for the velocity field:

v = f (θ) r , (1)
irrespective of the constitutive equation of the fluid. In this paper, our goal is two-fold. First, we document the quasistatic flow of a foam in a convergent channel, when wall friction plays a negligible role. In particular, we show that the velocity profile (1) holds, we show the clear connection between the velocity field and the distribution of plastic events, and we quantify the efficiency of plastic events as increments of deformation. Second, we study non-quasistatic flows at large wall friction, we quantify how the latter modifies both the velocity profile (1) and the elastic stress, and we discuss these effects in the context of recent models.

Materials and methods

Setup

We have adapted the setup described in [11]. The foam flows in a Hele-Shaw cell, made of two horizontal glass plates of length 170 cm and width 32 cm, separated by a gap h = 2 mm thin enough that the foam is confined as a bubble monolayer (Fig. 1b). Two slanted plastic plates of thickness 2 mm are inserted aside the Hele-Shaw cell, so that the remaining space within the channel constitutes the convergent channel, of opening half-angle α 0 = 26 • (Fig. 1b). The exit of the convergent channel has a width of 2.2 cm, much smaller than the width at its entrance (32 cm).

The channel is connected upstream to a vertical chamber (Fig. 1a) in which a soap solution is fed at a prescribed flow rate Q ℓ thanks to a pump. Nitrogen is continuously blown through injectors at the bottom of this chamber. The flow rate in each injector is independently controlled with a flow-rate controller. We identify the liquid fraction ϕ ℓ as the ratio of the liquid flow rate to the total flow rate:

ϕ ℓ = Q ℓ /(Q g + Q ℓ )
, with Q g the gas flow rate. The resulting foam accumulates on top of the chamber, then flows through the channel. The transit time through the whole channel is less than 10 minutes in all experiments; we do not observe significant changes of bubble size during this time, hence coarsening is negligible.

We have used two soap solutions. The first is a solution of sodium dodecyl sulfate (SDS) in ultrapure water at a concentration of 10 g/L, much larger than the critical micellar concentration of SDS (2.3 g/L): the foam was perfectly stable and no film rupture was observed. The surface tension of this solution is γ = 36.8 mN/m. The second solution is a mixture of sodium lauryl-dioxyethylene sulfate (SLES), cocoamidopropyl betaine (CAPB) and lauric acid (LAc), following the protocol described in [START_REF] Golemanov | [END_REF]. Its surface tension is γ = 23.8 mN/m. The solution of SDS has a negligible surface viscoelasticity, and the wall friction is relatively weak. On the other hand, the SLES/CAPB/LAc mixture has a large surface viscoelasticity, and the wall friction is large. We will see in Secs. 3 and 4 that this has a dramatic effect on the flow.

The contraction region is lit by a circular neon tube, giving an isotropic and nearly homogeneous illumination over most of the image. Movies of 1000 images of the foam flow are recorded with a CCD camera.

Image analysis

The movies are analysed using a home-made procedure fully described in [6]. Briefly, the velocity field is obtained after averaging of all the displacements of all individual Fig. 2. Snapshot of a plastic event. Two successive skeletonised images of the foam (grey, then black) are superimposed. The link between the centres of the two bubbles which are about to lose contact is figured by the red plain line, and the link between the centres of the two bubbles that have established contact is figured by the green dashed line. The angle β of the appearing contact with respect to the radial direction is also figured.

bubbles between consecutive frames. The plastic events, characterised by the neighbour-swapping of four adjacent bubbles (Fig. 2), are also automatically detected. This is not an easy task, especially because they are rare events, which are easily overwhelmed by artifacts if image analysis is not performed with extreme care. We estimated the relative uncertainty on the number of plastic events reported in this paper by visual inspection on a selection of images. Moreover, for dry 2D foams, the elastic stress σ equals ( [START_REF] Batchelor | [END_REF]; see [6] for a full discussion):

σ = 2γhρ ℓ ⟨ ℓ ⊗ ℓ ℓ ⟩ , (2)
where ℓ designs an edge vector, i.e. a vector joining the two vertices at the extremities of a bubble edge, and ρ ℓ is the areal density of such edges (its order of magnitude is thus 1/A). The brackets stand for the average over the time interval and space zone over which the elastic stress is computed. The ability to automatically measure plastic events and elastic stress by image analysis is a major advantage of dry 2D foams for studies of local rheology.

To compute averaged quantities like the fields of velocity, elastic stress and frequency of plastic events, we use a system of polar coordinates which reflects the symmetry of the channel. The origin is the convergence point of the channel (Fig. 1b), and the angular domain covering the channel is -α 0 ≤ θ ≤ α 0 , hence θ = 0 represents the axis of symmetry of the channel. The convergent channel occupies the range r ≥ 2.2/(2 sin 26 • ) = 2.5 cm and r ≤ 32/(2 sin 26 • ) = 36 cm (see Sec. 2.1 for the dimensions of the channel). We mesh the flow domain by boxes shaped as angular sectors, each with a radial length ∆r = 7.4 mm and a angular opening ∆θ = 3.4 • .

Results

Reference experiment with a SDS solution

We first report the measurements on a foam made of SDS solution, with the following control parameters: Q g = 100 mL/min, Q ℓ = 9.8 mL/min, hence ϕ ℓ = 8.9%, and A = 33.0 mm 2 . We henceforth call this experiment the "reference" experiment. We plot the velocity field in Fig. 3 and the elastic stress field in Fig. 4. As expected, the velocity field looks radial, and the velocity increases as the foam advances in the convergent channel. The elastic stress also obeys the radial symmetry: the long axes of the ellipses are aligned with the radial direction. The elastic stress does not seem to depend on the distance r along the channel.

We now plot the velocity components along the symmetry axis θ = 0 in Fig. 5. This plot shows that the velocity component v θ remains negligible compared to -v r , and that the absolute value |v r | = -v r of the radial velocity decreases at increasing r; notice that v r < 0 because the flow is convergent. Moreover, a log-log plot of -v r (r, θ = 0) (inset of Fig. 5) shows that its r-dependence is very close to a power law -v r ∝ 1/r, which is compatible with the velocity profile [START_REF] Cantat | Foams, Structure and Dynamics[END_REF]. Concerning end effects, the graph shows a saturation of -v r (r, θ = 0) at the two shortest distances r that we were able to measure (r = 18 and 26 mm): this is the signature of the transition from the convergent channel to the straight exit channel, which takes place at r = 25 mm (see Sec. 2.2). On the other hand, the entrance of the channel is at r = 36 cm (see Sec. 2.2), far larger than the largest distance r presented on Fig. 5, hence entrance effects are probably negligible.

We have checked (data not shown) that the two main features of Fig. 5: v θ negligible compared to -v r , and -v r ∝ 1/r, are recovered at various angular positions θ within the channel. To isolate a possible angular dependence, we therefore plot -rv r (averaged over 21 boxes) as a function of θ, in Fig. 6. It shows that the angular velocity profile is flat, with no significant deceleration towards the side walls.

Next, we plot the elastic stress components along the symmetry axis θ = 0 in Fig. 7. It shows that the shear component σ rθ is negligible compared to the diagonal components σ rr and σ θθ , which was expected due to symmetry. The radial component σ rr is always greater than the orthoradial one σ θθ , which quantifies the fact that the bubbles get elongated towards the exit of the channel (Fig. 4). More interestingly, the two diagonal components remain roughly constant across the channel. This saturation of the elastic stress means that the foam has yielded within our observation field, which suggests again that entrance effects are negligible. Indeed, entrance effects would be characterised by an elastic loading of the foam, hence we would measure an increase of σ rr and a decrease of σ θθ .

We now turn to plastic events. These are relatively rare events: we were able to detect 3 × 10 3 of them on the full movie, whereas there are about 10 6 individual data (about 10 3 bubbles on each of the 1000 images of the run) over which to average to determine velocity or elastic stress fields. Contrary to velocity and elastic stress, we will thus not represent a map of plastic events. We have checked that they are equally distributed along θ (data not shown); we thus directly average over θ, and define the T1 frequency f T 1 , in mm -2 •s -1 , as the number of T1s per unit time and per unit area: f T 1 (r, θ)r ∆r ∆θ is thus the frequency of T1s occurring in the box centred on (r, θ). We estimated a relative uncertainty of 10% on this measure. The T1 frequency is plotted as a function of r in Fig. 8. It increases at decreasing r, more sharply than the velocity, until the exit of the convergent channel is reached, then strongly decreases for the shortest investigated r which is in the straight exit channel. Excluding this end effect, a log-log plot of the T1 frequency shows that its r-dependence is compatible with a power law of exponent -2 (inset of Fig. 8). Not only the spatial distribution of the plastic events matters, but also does their orientation. To quantify it, we measure the angles between the radial direction e r and the direction of disappearing contacts, i.e. the direction linking two centres of neighbouring bubbles which are about to lose contact, and the direction of appearing contacts, i.e. the direction linking two centres of neighbouring bubbles which have just established contact (Fig. 2). For simplicity, we call these angles appearing and disappearing angles. Their range can be restricted to the interval [0, π/2]. Their distribution are plotted as histograms in Fig. 9. The appearing angles have their maximum at π/2, and the disappearing angles at 0. Hence, the bubbles that detach tend to be aligned streamwise, and the bubbles that attach tend to be aligned spanwise. The distribution of the appearing angles is significantly wider than the one of disappearing angles.

To conclude this presentation of the results with the SDS solution, we mention that we have performed extra experiments, changing systematically and independently the following parameters: the liquid flow rate Q ℓ (up to 500 mL/min), to change the average velocity; the gas flow rate (down to 1 mL/min), to change the liquid fraction down to 0.5%; and the bubble size, in the range 18-40 mm 2 . Within this range of parameters, we have not measured any significant modification to the main features of the foam flow, namely a velocity field of the form (1), and a saturation of the elastic stress as in Fig. 7.

Experiments with the SLES/CAPB/LAc mixture

The results obtained with the SDS solutions suggest that wall friction remains negligible. To quantify the effects of wall friction, we now present some measurements done with the SLES/CAPB/LAc mixture. It has a large interfacial dilatational viscoelasticity [START_REF] Golemanov | [END_REF], hence the gas/liquid interfaces are barely compressible. As a result, the bubbles tend to entrain these interfaces at their own velocity. Hence, the wetting films between the bubbles and the top, bottom and side plates are strongly sheared, and these foams develop a strong wall frictional shear stress [23].

We present an experiment with the same liquid fraction: ϕ ℓ = 8.9% as the experiment detailed in Sec. 3.1 with the SDS solution, and a bubble area: A = 36.7 mm 2 very similar (10% larger) to that experiment. Once this foam is formed, we push it at a larger flow rate (equal to 510 mL/min) to enhance the effect of wall friction. This flow rate is larger than the one of the experiment presented in Sec. 3.1 (110 mL/min). However, these two experiments remain comparable: as mentioned in Sec. 3.1, we have checked that the results presented there for SDS solutions do not change significantly when the flow rate is increased to 500 mL/min. We first plot the velocity as a function of r, for two angular positions: the symmetry axis θ = 0, and an angle θ = -25 • close to the left boundary of the channel, in Fig. 10. As expected, the velocity is a decreasing function of r. A first difference compared with the SDS solution (Fig. 6) is that -v r is significantly lower close to the side walls compared to the symmetry axis. Moreover, the decrease of -v r as a function of r is stronger at the symmetry axis. Close to the boundary, the decrease of -v r as a function of r clearly deviates from a power law behaviour, and slows down close to the exit of the convergent channel. Moreover, the angular profile depends on r: the inset of Fig. 10 shows that the larger r, the flatter the angular velocity profile v r (θ). Therefore, the velocity field shows significant deviations from the scaling (1).

As we did for the SDS solution (Fig. 7), we plot the elastic stress components along the symmetry axis θ = 0 in Fig. 11. Once again, the shear component σ rθ is negligible compared to the diagonal components σ rr and σ θθ , and the radial component σ rr is greater than the orthoradial one σ θθ . But contrary to the SDS solution, these components do not remain constant across the channel: σ rr tends to decrease, and σ θθ to increase, at increasing r. As a consequence, the deformation of the bubbles tends to increase as they flow towards the exit of the convergent channel. The bubbles are always more deformed than with the SDS solution (Fig. 12).

Finally, we recorded the T1 frequency for this experiment. It turned out to be significantly noisier than for the SDS solution; the relative uncertainty was estimated to be 30%, and the data for r ≥ 150 mm are too noisy to be exploited. We present and discuss these data in Sec. 4.2.

Discussion

A model plastic flow?

Our experimental results show that for a foam experiencing a negligible wall friction (see Sec. 4.2 for more discussion on the influence of wall friction), the kinematics of the flow in a convergent channel follows the very simple law (1) (Fig. 3 and5).

This result does not come as a surprise, since the foam is quasi-incompressible. Since it is mostly constituted by gas, the order of magnitude of its bulk modulus is given by the atmospheric pressure, equal to 10 5 Pa. The order of magnitude of its shear modulus, given by γ/ √ A [START_REF] Cantat | Foams, Structure and Dynamics[END_REF], is much lower: with γ = 36 mN/m and A = 33 mm 2 , γ/ √ A ≈ 10 Pa. The deviation of the velocity profile from (1) in the presence of wall friction is discussed in Sec. 4.2.

The fact that the velocity field of the reference experiment obeys (1) makes it a model to study the interplay between plasticity and flow. More precisely, at constant elastic stress, the deformation rate must be directly correlated to the rate of plastic events. Our experimental results strongly suggest that it is indeed true: the rate of plastic events, measured independently from the velocity field, depends on r with a power law compatible with 1/r 2 (inset of Fig. 8), which is precisely the scaling of the deformation rate associated with the velocity field (1). Furthermore, the preferential orientations of the plastic events (Fig. 9) also support this idea: bubbles preferentially detach along the direction of elongation e r , whereas they preferentially attach along the direction of compression e θ .

Let us go further and compare the prefactors. A single plastic event is expected to separate the two bubbles that detach by a distance of the order of the bubble size, as Fig. 2 already suggests. Hence, we may expect that the streamwise deformation rate, ∂v r /∂r, is of the order of the rate of plastic events times the bubble size. To test this idea, we fit the data again, but this time we fix the exponents. We fit the velocity by a law: log(-v r ) = log c 1 -log r, and the T1 frequency by a law: log f T 1 = log c 2 -2 log r. We find as best fitting parameters c 1 = (1.03 ± 0.01) × 10 3 mm 2 /s, and c 2 = 25.9 ± 2.2 s -1 . The ratio of these two parameters gives a characteristic area A T 1 = 40 ± 4 mm 2 , which is indeed very close to the bubble area A = 33.0 mm 2 , such that:

∂v r ∂r = A T 1 f T 1 .
This confirms quantitatively the kinematic consequence of T1s, namely to incrementally shift bubbles streamwise, by inserting additional bubbles in between (Fig. 2).

It is interesting to point out that this characteristic area is a somewhat different concept from the range of plastic events which has been abundantly discussed in the context of plastic flows [35], or from the size of a plastic event which can be estimated by indirect methods like diffusive wave spectroscopy [36]. The area A T 1 characterises the "efficiency" of a given plastic event as an increment of deformation, but does not describe its effect in its neighbourhood and the way it redistributes the elastic stress around [37,38]. These two kinds of information are thus complementary.

Influence of wall friction

Velocity field

The reference velocity field to which our results have been compared is given by (1): v = f (θ)e r /r. We have seen that while the flow of a foam constituted by SDS solution complies with such a profile, the foam made by a SLS/CAPB/LAc mixture shows significant deviations. This suggests that wall friction affects the velocity profile.

To start the discussion, we notice that (1) is a selfsimilar profile. Self-similarity is expected to hold when no characteristic length appears in the flow problem. For the flow of a Newtonian fluid in a convergent channel, the flow profile depends only on the Reynolds number Q 2D /ν, with Q 2D the flow rate per unit length, and self-similarity holds. In particular, the no-slip boundary condition does not introduce extra parameters. For the foam with the SDS solution, the wall friction is negligible also at the side walls, hence free shear boundary condition: ∂v r /∂θ = 0 at θ = ±α 0 applies to a good approximation (Fig. 6). Like the no-slip boundary condition, the free shear boundary condition is compatible with the self-similar profile (1). However, as soon as wall friction is introduced, another characteristic coefficient appears, and the self-similarity is lost.

Predicting the velocity profile for such a viscoelastoplastic material is a daunting task. However, it has been shown in pure shear with no-slip boundary condition [5] that elasticity may have a negligible effect on the velocity profile, which can then be predicted analytically from a balance between viscous stress and wall friction. On the other hand, such a simplified approach fails to predict correctly the slip velocity along smooth walls in a Poiseuille flow [6]. In our geometry, we derive the following velocity profile neglecting elasticity in Appendix A:

v r = -Q/h 2rα 0 -2hℓ/(h + 2ℓ) ( 1 - 2h h + 2ℓ
e -rα0/ℓ cosh rθ ℓ

) ,

(3) where Q = Q g + Q ℓ is the foam flow rate and ℓ is a characteristic length comparing the effect of wall friction and viscous stress; this expression is valid only if ℓ/r ≪ 1, which we argue is the case for SLES/CAPB/LAc mixtures in the Appendix A.

The velocity profile (3) shows that wall friction is responsible for a deviation from [START_REF] Cantat | Foams, Structure and Dynamics[END_REF]. It also captures the fact that the velocity at the side walls is lower than at the centre of the channel, because of the extra friction at the side walls. However, it predicts that v r (r, θ = α 0 )/v r (r, θ = 0) is a decreasing function of r (see Appendix A for a proof), whereas the opposite occurs in experiments (Fig. 10) for the SLES/CAPB/LAc mixture, namely the fact that v r (r, θ = α 0 )/v r (r, θ = 0) increases at increasing r. Indeed, it is the same failure as in the Poiseuille flow of foams, namely that the slip velocity is not well predicted [6]. There can be three reasons for this discrepancy: (i) the dependencies of wall friction on velocity on one hand, and of viscous stresses on deformation rate on the other hand, may differ, (ii) the hypothesis ℓ ≪ r breaks down before the channel exit, (iii) nonlocal effects, more precisely the redistribution of elastic stresses by plastic events, may matter.

Elastic stress

Fig. 11 shows that for the SLES/CAPB/LAc mixture, the normal elastic stress difference σ rr -σ θθ increases as the bubbles flow faster towards the exit of the channel. This effect has been previously mentioned in experiments of a foam flowing through a constriction [11]. Similarly, simulations of foam including non-quasistatic rate effects, either a variation of surface tension in the so-called vertex model [14], or an external velocity-dependent friction in the socalled viscous froth model [START_REF] Embley | [END_REF], report an increase of the elastic deformation as the velocity or the deformation rate increases.

To rationalise this fact, it has been proposed that the elastic stress tends to increase as the relaxation time τ r following a plastic rearrangement becomes significant compared to the deformation rate ε [43]. However, this dependence of an elastic stress on the deformation rate has been mostly discussed in simple shear, and only in theory or in simulations. We thus propose to quantify this dependence in our experiments.

First of all, we need an estimate of the deformation rate ε relevant to our experiments. The flow kinematics is more complicated than simple shear or pure elongation, but a representative deformation rate is given by the axial elongation rate, i.e. the spatial variation of the velocity along the symmetry axis. Hence, we will take ε = |∂v r /∂r| θ=0 . To estimate this quantity from our measurements of velocity, we first fit v r (r, θ = 0) by a power law: v r (r, θ = 0) = βr α , and we compute the deformation rate as the derivative of this fit: ε = βαr α-1 . An advantage of our setup is that the deformation rate varies by two orders of magnitude on a given experiment between the entrance (r = 32 cm) and the exit (r = 2.5 cm) of the convergent channel. Since the flow is dominated by elongation, the relevant elastic quantity to study is the normal stress difference σ rr -σ θθ . To make the comparison easier with other experiments and other studies, we divide the normal elastic stress difference by γ/ √ A as in [43], and we call the resulting dimensionless quantity (σ rr -σ θθ /(γ/ √ A) the "elastic deformation". Next, we quantify the relaxation time τ r . We proceed as described in [44], by measuring the evolution of the length of one newly created film after a T1. Even if the relaxation process might be complicated [44,45], Fig. 13 shows that the length of such a new film follows closely an exponential relaxation towards a "final" length (although this length drifts on longer times, as the bubbles concerned by the T1 are advected away), and the characteristic time of this exponential behaviour yields: τ r = 0.019 s for the SDS solution, and 0.13 s for the SLES/CAPB/LAc mixture. We have checked on some other examples of newly created films that these values remain correct within 30%. Hence, we take these two values as good estimates of the relaxation time following a T1.

For the experiment with the SDS solution presented in Sec. 3.1, the deformation rate equals at most 1.7 s -1 , hence ετ r does not exceed 0.03. Hence, the relaxation time after a T1 is everywhere negligible compared to the deformation rate, and the flow remains quasistatic, which is consistent with the absence of variation of the elastic stress components along the channel (Fig. 7). On the other hand, for the experiment with the SLES/CAPB/LAc mixture presented in Sec 3.2, we compute that ε increases from 0.10 to 11 s -1 , hence ετ r from 0.013 to 1.4, along the channel. Hence, as r decreases, the relaxation after plastic rearrangements is no longer short compared to the loading by the flow, which is consistent with the increase of the elastic deformation on Fig. 11.

To go further, we can test whether the elastic deformation can be rescaled by ετ r irrespective of the location r within the channel. To improve the comparison, we com-Fig. 13. Evolution of the length of a film created after a T1, for the SDS solution (squares) and the SLES/CAPB/LAc mixture (circles), and three snapshots of this film, represented by a double arrow. The time origin was chosen as the first image where this length could be measured reliably. The data were fitted by an exponential relaxation of the form ℓ(t) = ℓ1 + ℓ2e -(t-t 0 )/τr , with ℓ1, ℓ2, t0 and τr as fitting parameters; the latter parameter is identified as the relaxation time following the T1. It equals 0.019 ± 0.001 s for the SDS solution, and 0.13 ± 0.01 s for the SLES/CAPB/LAc mixture.

pile the results of five experiments at different gas flow rates (100, 250 and 500 mL/min), but at the same liquid fraction ϕ ℓ = 8.9% and the same area A = 35 mm 2 within 10%, and we plot the elastic deformation as a function of ετ r in Fig. 14. This plot shows that all the data for the SLES/CAPB/LAc mixture collapse on the same master curve which is well fitted by a Herschel-Bulkleylike law: (σ rr -σ θθ )/(γ/ √ A) = σ * [1 + ( ετ r /ε * ) n ] with a characteristic exponent n = 0.40 ± 0.04. This may remind the Herschel-Bulkley law followed by the flow curves of SLES/CAPB/LAc foams [START_REF] Golemanov | [END_REF], which exhibit a slightly lower exponent (about 0.25). However, it is irrelevant to push forward such a comparison, since we have no access to the viscous stress. The limit σ * of this law at vanishing deformation rate is also compatible with the data for the SDS solution for ετ r lower than a few 10 -3 , suggesting that this regime is the quasistatic elastoplastic limit, independent on the solution except on its surface tension. However, the data for the SDS solution deviate from those with the SLES/CAPB/LAc mixture at larger ετ r , and show little dependence on ετ r . This suggests that the departure from the quasistatic elastic deformation is not universal, but depends on the viscous properties of the solution, and not only through the relaxation time after T1s.

Rate of plastic events

In Sec. 4.1, we showed that at constant elastic deformation, the rate of plastic events is proportional to the deformation rate. Since the experiment with the SLES/CAPB/LAc mixture presents a significant increase of the elastic deformation towards the channel exit (Fig. 11), it is interesting to investigate the relationship between the T1 frequency and the deformation rate for this experiment too.

The velocity fields depend on wall friction (Figs. 5 and 10), but overall, the deformation rate is proportional to the foam flow rate Q. Hence, a simple way to test whether the elastic deformation has an influence on the rate of plastic events at given deformation rate is to plot the ratio f T 1 /Q as a function of r (Fig. 15). Since f T 1 is averaged over θ, this approach is not sensitive to the angular dependence of the radial velocity discussed for the SLES/CAPB/LAc mixture (inset of Fig. 10). Fig. 15 shows that f T 1 /Q is lower for this experiment than for the SDS solution. This clearly evidences an influence of the elastic deformation on the rate of plastic events.

Recent models [35,46] claim that the rate of plastic events should be proportional, not to the sole deformation rate, but to the so-called fluidity, which is the ratio of the deformation rate to the stress. Our results are indeed in qualitative agreement with this statement. However, it is not possible to make a quantitative comparison, for the following reasons. (i) It is the total stress (elastic plus viscous) that appears in the fluidity, while we can measure only the elastic contribution. (ii) Fluidity models have been devised for pure shear only, and our flow combines shear and elongation. (iii) The extension of fluidity models in the presence of wall friction is not trivial [47,48].

Conclusions

In this paper, we have quantified in details the 2D flow of a foam through a convergent channel. We have showed that if wall friction is negligible, the velocity field obeys a simple self-similar form. It makes this experiment a model of plastic flow: we could establish a clear connection between the velocity field and the distribution of plastic events, showing that all of them act efficiently as increments of flow.

The use of the SLES/CAPB/LAc mixture enabled us to quantify the departure from this simple picture as wall friction becomes significant. We quantified for the first time in an experiment how the elastic deformation increases when the deformation rate increases to become of the same order as the inverse of the relaxation time following plastic events. We showed how the velocity field is modified, with a significant deceleration of the foam close to the smooth side walls, and a breakdown of the self-similar form. This qualitative picture is predicted by a model balancing wall friction and viscous stress, but this model fails to capture the variation of the ratio of the slip velocity to the centreline velocity. Finally, we showed that the elastic deformation influences the relationship between the rate of plastic events and the deformation rate.

We believe that our work opens several interesting perspectives. First, it would be interesting to find the missing ingredient to offer a full quantitative prediction of the velocity profile, including elastic effects, either in the frame of a viscoelastoplastic model [49] or as an effective nonlocal term [50]. In the same spirit but more on the experimental side, we could investigate the spatial effects of plastic events, and how they modify the elastic stress in their surroundings [38]. Moreover, if we showed that the elastic deformation was an increasing function of the deformation rate for the SLES/CAPB/LAc mixture, it was not as clear for the SDS solution: it would be interesting to understand this difference, and more precisely to iden-tify unambiguously the parameter that controls the departure from the quasistatic limit. Finally, we have observed in some cases that increasing even more wall friction e.g. by increasing the velocity leads to brittle fracture of the foam. Fast confined flows could therefore be an interesting alternative way to study fracture in foams.

We acknowledge discussions with Axelle Amon, Isabelle Cantat and Kirsten Martens. We thank Luca Biferale, Massimo Cencini, Alessandra Sabina Lanotte and Mauro Sbragaglia for the organisation of the workshop "Flowing Matter Across the Scales" which motivated this study.

A Derivation of the velocity profile neglecting elasticity

In the Appendix, we predict the velocity profile of a foam flowing through a convergent channel when neglecting elastic stresses. We revisit the model by Janiaud et al. [39], who first studied the interplay between viscous stress and wall friction in the context of a plane Couette flow. Because of its low density and large effective viscosity, foam has a negligible inertia in the studied flow, which is furthermore steady. With these assumptions, the model by Janiaud et al. writes:

0 = ∇ • τ -∇p + 2 h f , ( 4
)
where τ is the stress tensor within the foam, p the pressure field and f the wall friction force per unit area.

For simplicity, we assume that the viscous stresses scale linearly with the deformation rate: τ = η(∇v + ∇v T ), with η the foam effective viscosity, and that the wall friction force scales linearly with the velocity: f = -hβv/2, with β a friction coefficient. In reality, these scalings are sublinear [23], but nonlinear corrections were shown not to modify much flow profiles in Couette flows [40] and Poiseuille flows [6], and we expect the same here. With these assumptions, (4) becomes the Stokes equation with an additional term accounting for wall friction:

0 = η∆v -∇p -βv. ( 5
)
Because wall friction introduces a new term and a new coefficient β, self-similarity does not hold: the velocity profile (1) is not solution of (5).

As boundary conditions, we require that θ = 0 is an axis of symmetry, hence ∂v r /∂θ = 0 at θ = 0. Furthermore, extra friction at the side walls imposes that τ rθ = ±f • e r at θ = ±α 0 , hence:

1 2 hβv r = ∓ η r ∂v r ∂θ , θ = ±α 0 . (6)
Finally, conservation of flow rate imposes the integral condition:

Q = -h ∫ α0 -α0 rv r dθ, ( 7
)
where Q is the volumetric flow rate, as defined in the experiments. The minus sign reminds that when Q > 0, then v r < 0.

To solve the problem, we take the curl of (5) to remove the pressure. For a 2D flow, the vorticity ω = ∇ × v is a scalar field: ω = ωe z , hence we get:

η∆ω -βω = 0. ( 8
)
This is a Helmholtz equation which can be solved analytically [START_REF] Prosperetti | Advanced Mathematics for Applications[END_REF], but we will content ourselves with a qualitative discussion and an approximate solution.

First, if wall friction is negligible as is the case for SDS solution, then Eq. ( 8) reduces to ∆ω = 0, and the boundary condition (6) reduces to ∂v r /∂θ = 0 at θ = ±α 0 . Together with the symmetry condition ∂v r /∂θ = 0 at θ = 0, this imposes that v r does not depend on θ. Therefore, the conservation of flow rate (7) imposes the velocity profile:

v r = - Q 2α 0 hr ,
which is a particular case of the self-similar law [START_REF] Cantat | Foams, Structure and Dynamics[END_REF]. If wall friction is not negligible, Eq. ( 8) contains a characteristic length: ℓ = √ η/β which is the shear localisation length for the pure shear flow of confined foams [39,5]. In the convergent geometry, it is relevant to compare it with the distance r. In our experiments, we have seen that the exit of the channel is at r = 2.5 cm. In previous experiments with comparable foams [6], the length ℓ was found to be around 0.5 mm. Hence, we now focus on the case ℓ/r ≪ 1.

A naive inspection of ( 8) suggests the ratio of the viscous term η∆ω ≈ ηω/r 2 to the wall friction term is of order ℓ 2 /r 2 ≪ 1, hence wall friction dominates and (8) reduces simply to:

ω = 0, (9)
i.e. an irrotational flow. However, this is incompatible with the boundary condition (6), which shows that sidewall friction is a source of vorticity. Therefore, there is a shear boundary layer of thickness ℓ close to the side walls where friction and viscous stresses balance. In such boundary layers, since ℓ ≪ r, the axial variation of the vorticity is negligible compared to its angular variation, and (8) simplifies into:

1 r 2 ∂ 2 ω ∂θ 2 - ω ℓ 2 = 0,
where r can be treated as a constant. Hence, the vorticity has an exponential profile: ω = ∓ω 0 e -r(α0∓θ)/ℓ close to the side wall θ = ±α 0 , which matches the irrotational flow (9) as r(α 0 ∓ θ)/ℓ ≫ 1. Hence, the vorticity profile equals: ω = ω 0 [-e -r(α0-θ)/ℓ +e -r(α0+θ)/ℓ ] = -2ω 0 e -rα0/ℓ sinh rθ ℓ .

Similar to the Laplacian term, the general relation between the vorticity and the velocity:

ω = 1 r [ ∂ ∂r (rv θ ) - ∂v r ∂θ ]
simplifies to ω = -1 r ∂vr ∂θ , hence v r = v 0 + 2ℓω 0 e -rα0/ℓ × cosh(rθ/ℓ). Inserting this expression in the boundary condition (6) yields a first relation between the two unknowns ω 0 and v 0 : ω 0 = -hv 0 /(2ℓ 2 + ℓh). The second relation comes from the conservation of flow rate (7). Finally, the velocity field is:

v r = -Q/h 2rα 0 -2hℓ/(h + 2ℓ) ( 1 - 2h h + 2ℓ
e -rα0/ℓ cosh rθ ℓ

) , which strictly speaking is valid only if ℓ/r ≪ 1, and up to corrections of order ℓ/r. Using this equation, it is straightforward to show that: ×(1 -e -rα0/ℓ )(2ℓ -he -rα0/ℓ ), which is negative for all values of r > 0 if 2ℓ > h. If 2ℓ < h, this derivative is negative if e -rα0/ℓ < 2ℓ/h, i.e. r > -ℓ α0 ln 2ℓ h . The latter quantity does not exceed 0.368h/2α 0 = 0.8 mm (the factor 0.368 is the opposite of the minimum value taken by the function x ln x for x between 0 and 1). Hence, over our experimental range of r, the velocity profile (3) implies that v r (r, θ = α 0 )/v r (r, θ = 0) is a decreasing function of r.

Fig. 1 .

 1 Fig. 1. (a) Sketch of the side view of the setup. For clarity, the slanted plates limiting the convergent channel have not been drawn. (b) Top view of the convergent channel, with the definition of the polar coordinates used throughout the text.

Fig. 3 .

 3 Fig. 3. Velocity field for the reference experiment.

Fig. 4 .

 4 Fig. 4. Elastic stress field for the reference experiment. The elastic stress (2) is a symmetric and positive tensor, hence it possesses two positive eigenvalues with orthogonal eigenvectors. It can thus represented as en ellipse, with a major (minor) axis along the direction of maximal elongation (compression).

Fig. 5 .Fig. 6 .

 56 Fig. 5. Velocity components -vr (circles) and v θ (squares) as a function of r, along the central axis θ = 0. Inset: log(-vr) as a function of log r, and its best linear fit, which slope is -1.03 ± 0.02.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Elastic stress components σrr (circles), σ θθ (lozenges) and σ rθ (squares) as a function of r, along the central axis θ = 0.

Fig. 9 .

 9 Fig. 9. Histograms of the distributions of the appearing angles (in green) and of the disappearing angles (in red).

Fig. 10 .

 10 Fig. 10. Plot of log vr as a function of log r, for θ = 0 (circles) and θ = -25 • (squares). The dashed line is a guide for the eye of slope -1. Inset: plot of vr/⟨vr⟩ θ as a function of θ, where ⟨vr⟩ θ is the average of the data over the angle θ, for r = 26 mm (lozenges) and 122 mm (triangles).

Fig. 11 .

 11 Fig. 11. Elastic stress components σrr (circles), σ θθ (lozenges) and σ rθ (squares) as a function of r, along the central axis θ = 0.

Fig. 12 .

 12 Fig. 12. Snapshots of the foam at the exit of the convergent channel, (a) for the SDS solution, (b) for the SLES/CAPB/LAc mixture. The bubbles are visibly more elongated streamwise in the latter case.

2 )Fig. 14 .

 214 Fig.14. Plot of the elastic deformation (σrr -σ θθ )/(γ/ √ A) as a function of ετr for foam flows with a SDS solution (full symbols) or a SLES/CAPB/LAc mixture (empty symbols), for Qg = 100 (circles), 250 (lozenges) and 500 mL/min (squares). The curve is a single fit of all data for the SLES/CAPB/LAc mixture by a Herschel-Bulkley-like law: (σrr-σ θθ )/(γ/ √ A) = σ * [1 + ( ετr/ε * ) n ],with best fitting parameters σ * = 0.058 ± 0.004, ε * = 0.17 ± 0.06 and n = 0.40 ± 0.04.

Fig. 15 .

 15 Fig. 15. Plot of the ratio of the T1 frequency to the foam flow rate as a function of r, for the experiment with the SDS solution (circles) and the SLES/CAPB/LAc mixture (squares).

v

  r (r, θ = α 0 ) v r (r, θ = 0) = h + 2ℓ -h(1 + e -2rα0/ℓ ) h + 2ℓ -2he -rα0/ℓ ,and that:d dr [ v r (r, θ = α 0 ) v r (r, θ = 0) ] = -2α 0 h ℓ e rα0/ℓ [(h + 2ℓ)e rα0/ℓ -2h] 2


Kazusa Yamamoto

Damien Koenig

Olivier Sename

Pascal Moulaire

Driver torque estimation in Electric Power Steering system using an H ∞ /H 2 Proportional Integral Observer

This paper deals with the design of a Proportional Integral (PI) observer to estimate the driver torque in an Electric Power Steering (EPS) system. The PI observer is obtained by solving a multi-objective optimization problem: it should both be barely sensitive to road disturbances and sensor noise, and converge swiftly. The performance of the proposed observer is illustrated by simulation results using experimental data.

I. INTRODUCTION

Currently, most of the vehicles uses Electric Power Steering (EPS) systems. Indeed, compared to hydraulic steering systems, EPS systems improve fuel efficiency (with a fuel economy up to 5%), quality of feedback to the driver and ease of integration in the vehicle.

In modern vehicles, steering systems help the driver turn the vehicle, improving both safety and comfort. In order to do so, an assistance torque is provided by an electric assistance motor, to reduce the amount of torque required from the driver to turn the wheels in the desired direction. The amount of supplied power is defined by the Electronic Control Unit (ECU) according to a motor torque control policy, containing an assistance rule depending on the vehicle speed and measurements from a torque sensor [START_REF] Sulakhe | Electric Power Assisted Steering[END_REF].

A. Problem statement

Reducing production costs (e.g by removing a sensor) and improving the performance of EPS control system is a challenge for EPS system suppliers in the automotive market. One of the main requirements in EPS systems is to determine the amount of assist torque to provide. This is usually done using a torque sensor to measure applied steering torque. However, failure of this sensor could lead to a sudden loss of steering assistance. For driving comfort and safety reasons, such an event should be avoided. Moreover, it may be noted that the output of the torque sensor is reliable while the driver torque is constant and the assist motor is not providing an additional torque [START_REF] Marouf | Sensorless control of Electric Power Assisted Steering system, 20th Mediterr. Conf. Control Autom[END_REF], [START_REF] Chabaan | Torque Estimation in Electrical Power Steering Systems[END_REF]. Developing a controller based on sensorless driver torque estimation may offer improvements in both production costs and performance by providing a better estimation of the actual value of the driver torque [START_REF] Chabaan | Torque Estimation in Electrical Power Steering Systems[END_REF], [START_REF] Marouf | A new Control Strategy of an Electric Power Assisted Steering system[END_REF]. 1 Kazusa Yamamoto, Damien Koenig, Olivier Sename are with GIPSA-Lab, Control system department, Université Grenoble-Alpes, 11 rue des Mathematiques, BP46, 38402 Saint Martin d'Heres Cedex, France (kazusa.yamamoto,damien.koenig,olivier.sename) @gipsa-lab.grenoble-inp.fr

B. State of the art

Several studies have already been carried out on driver torque estimation in EPS systems. Marouf et al. proposed in [START_REF] Marouf | Sensorless control of Electric Power Assisted Steering system, 20th Mediterr. Conf. Control Autom[END_REF], [START_REF] Marouf | A new Control Strategy of an Electric Power Assisted Steering system[END_REF], [START_REF] Marouf | Driver torque and road reaction force estimation of an electric power assisted steering using sliding mode observer with unknown inputs[END_REF], a Sliding Mode Observer (SMO) to estimate both the driver torque and the road reaction force in an EPS system. However, it does not satisfy the observer matching condition, since additional inputs for the observer are generated using High Order Sliding Mode Differentiator (HOSMD). The SMO is based on a simplified EPS statespace representation and uses as inputs the steering wheel angle and the absolute angular position of the motor ("motor angle" in the remainder of this paper). [START_REF] Marouf | Sensorless control of Electric Power Assisted Steering system, 20th Mediterr. Conf. Control Autom[END_REF] studies a resolverless EPS system, in which motor angle is computed using HOSMD on the measured motor current, whereas in [START_REF] Marouf | A new Control Strategy of an Electric Power Assisted Steering system[END_REF], [START_REF] Marouf | Driver torque and road reaction force estimation of an electric power assisted steering using sliding mode observer with unknown inputs[END_REF] the two angles are measured. In [START_REF] Tordesillas Illán | Oscillation annealing and driver/tire load torque estimation in Electric Power Steering systems[END_REF], an observer computing the driver and load torque is designed by pole placement. In this case, the measured inputs of the observer are the torque sensor signal and the steering column speed. Another estimator is introduced in [3] and [START_REF] Chabaan | Control of Electrical Power Assist Systems : H∞ Design Torque Estimation and Structural Stability[END_REF] using the torque sensor signal and the assist motor input current. The estimator is deduced from the EPS transfer function model. However, in order to get a stable transfer function in the specified frequency range, approximations have been made in computing the estimator. In [START_REF] Cholakkal | Fault tolerant control of electric power steering using robust filter-simulation study[END_REF] an indirect estimation is used to compute the driver torque, since it is deduced from estimation of the road reaction torque and the motor shaft torque. These torques are computed using a model of the vehicle dynamics and a Direct Current (DC) motor model. In [START_REF] Parmar | Modeling and sensorless optimal controller design for an electric power assist steering system[END_REF] a Kalman filter is implemented to estimate all the states of a dual-pinion EPS system. Motor angle and voltage are the required measurements.

These previous studies show that driver torque estimation is still an open issue. The challenge is to avoid the need of a torque sensor to estimate the driver torque. Moreover the estimation algorithm should be implementable under some software/hardware constraints e.g ECU storage space. Moreover, the implementation of the observer on a massproduced vehicle requires developing new safety concepts with respect to embedded complexity, which is made easier by using a linear observer.

C. Contribution

The main contribution of this paper is the applicationoriented design of a PI observer to estimate the driver torque on an EPS system. The steering wheel angle and the assist motor angle are required inputs measurements for the PI observer. Then, the observer gains are chosen to get a fast estimation and to minimize the effects of road disturbances (resp. sensor noise) on the driver torque estimation error, in an H ∞ framework (resp. H 2 framework).

This paper is organized as follows. Section II presents the EPS model limited to its mechanical relations. Section III describes the multi-objectives design of the PI observer, and section IV shows synthetic results based on numerical simulation. Then, performance of the PI observer are illustrated by simulation results using experimental data in section V. Finally, section VI leads to the conclusions and future works.

II. EPS SYSTEM MODEL

Mechanical model of EPS

The EPS system can be decomposed into four main parts: the steering wheel, the assist motor, the pinion and the rack. The pinion type EPS (P-EPS) simplified mechanical model structure is illustrated in Figure 1 similarly as in [START_REF] Badawy | Modeling and Analysis of an Electric Power Steering System[END_REF].

J c θc = τ d -D tb ( θc -θp ) + K tb (θ c -θ p ) -B c θc (1) J m θm = τ m -D g ( θm -G θp ) + K g (θ m -Gθ p ) -B m θm (2)
J p θp = D tb ( θc -θp ) + K tb (θ c -θ p ) + G D g ( θm -G θp ) + K g (θ m -Gθ p ) -R p D pr (R p θp -Ẋr ) + K pr (R p θ p -X r ) (3)
J r Ẍr = τ road + D pr (R p θp -Ẋr ) + K pr (R p θ p -X r ) -D r Ẋr + K r X r -B r Ẋr (4)
where the variables are θ c the steering wheel angle (rad), θ m the motor angle (rad), θ p the pinion angle (rad), X r the rack position (m) and the mechanical parameters are described in Table I.

Hence the state-space representation of the P-EPS system In this paper the model identification and validation are not detailed. However, the given model was validated upon experimental data on the relevant frequency domain. Figure 2 shows the typical experimental and simulated frequency responses (computed as described in [START_REF] Poussot-Vassal | A new semi-active suspension control strategy through LPV technique[END_REF]) used to validate the model. In this case, the signal is applied as the motor torque input (control signal). Using a sinus of amplitude 10N.m with frequency from 0.1 to 60Hz, the obtained transfer function from control signal to motor angle T θmu is shown below in Figure 2.

ẋ = Ax + Bu + Ed + W w y = Cx + N n (

III. PI OBSERVER DESIGN

The main objective is to estimate the unmeasured driver torque input d, considering the unknown input w as a disturbance. Such a problem is similar to the unknown input observer design problem as described in [START_REF] Hui | Observer Design for Systems With Unknown Inputs[END_REF], [START_REF] Xiong | Unknown disturbance inputs estimation based on a state functional observer design[END_REF]. However, the system (5) does not fulfill the exact disturbance decoupling condition. Therefore, a H ∞ /H 2 PI observer is proposed to minimize the effect of sensor noise and road disturbances while providing a fast convergence of the estimation. Moreover, it should be ensured that driver torque acts in low frequency, since the PI observer dynamic allows to estimate constant and slowly time varying input.

In order to specify the frequency range for disturbance attenuation, a weighting function W f (s) is introduced as:

w = W f (s) w ⇔ ẋw = A w x w + B w w w = C w x w + D w w (6)
Indeed, the road reaction force can reach up to 10kN for parking manoeuvres. This means even an attenuation of -60dB implies an error on the estimation up to 10N.m, whereas the driver torque doesn't exceed 30N.m. The weighting function is therefore used to specify a frequency domain where the disturbance w should be attenuated.

The augmented representation combining ( 5) and ( 6) is:

ẋa = A a x a + B a u + E a d + W a w y = C a x a + N n (7)
where x a = x T x T w T are the augmented states and the associated matrices are

A a = A W C w 0 A w , W a = W D w B w , E a = E 0 , B a = B 0 , C a = C 0 .
The PI observer for the augmented representation ( 7) is written as follows :

ẋa = A a xa + E a d + L p (y -C a xa ) + B a u ḋ = L i (y -C a xa ) (8)
The condition for the existence of the observer is that the pair and assuming ḋ ≈ 0 (however this might be a restrictive hypothesis for evasive manoeuvres that imply fast driver dynamic), the dynamic of the estimation error deduced from [START_REF] Chabaan | Control of Electrical Power Assist Systems : H∞ Design Torque Estimation and Structural Stability[END_REF] and ( 8) is: ėad = (A ad -L a C ad )e ad + W ad w + L a N n z = D ad e ad [START_REF] Badawy | Modeling and Analysis of an Electric Power Steering System[END_REF] where L a = L p L i T is the observer gain to be deter-mined, z = d -d is the observer output and the associated

A a E a , C a 0 is detectable i.e rank   pI -A a -E a 0 pI C 0   = n xa + n d (9
matrices are A ad = A a E a 0 0 , W ad = W a 0 , C ad =
C a 0 , D ad = 0 1 . From [START_REF] Badawy | Modeling and Analysis of an Electric Power Steering System[END_REF], A ad -L a C ad is Hurwitz if and only if the pair (A ad , C ad ) is detectable, or equivalently:

rank     pI -A -W C w -E 0 pI -A w 0 0 0 pI C 0 0     = n x + n w + n d (11)
∀p ∈ C such that Re(p) ≥ 0. Here n x and n w denotes the number of states x in (5) and x w in [START_REF] Tordesillas Illán | Oscillation annealing and driver/tire load torque estimation in Electric Power Steering systems[END_REF]. Since a stable weighting function is chosen for ( 6), [START_REF] Poussot-Vassal | A new semi-active suspension control strategy through LPV technique[END_REF] comes to:

rank   pI -A -E 0 pI C 0   = n x + n d (12)
for all p such that Re(p) ≥ 0. For (5) the condition ( 12) is satisfied thus the observer described in (8) exists.

Figure 3 shows the augmented representation and the observer. Plant inputs are the unknown input w specified on a frequency range, the driver torque d to estimate and the control signal u. The PI observer inputs are the control signal u, the measured plant output y affected by sensor noise n. Its output is the estimated driver torque d. As explained below, the first design objective is to ensure the performances of the estimation through the specifications of the observer poles. On the other hand, even if the observer estimates all the internal states xa , the performance objective will focus only on the minimization of the driver torque estimation error d -d which explains the choice of z. Then, the design of the observer aims at minimizing the effect of road disturbances (resp. sensor noise) subject to H ∞ -norm (resp. H 2 -norm) on the driver torque estimation error. This needs to handle the sensitivity functions deduced from [START_REF] Badawy | Modeling and Analysis of an Electric Power Steering System[END_REF] :

T z w = D ad (sI -(A ad -L a C ad )) -1 W ad (13)
T zn = D ad (sI -(A ad -L a C ad )) -1 L a N (14)

A. Pole placement

As the observer dynamic response is defined by the eigenvalues of (A ad -L a C ad ), setting an upper bound λ min on the real part of the observer poles ensures a fast convergence of the estimate. Indeed, the other poles are in that case greater than λ min . Such a constraint is formulated as the following Linear Matrix Inequality (LMI) condition developed in [START_REF] Dugard | Full vertical car observer design methodology for suspension control applications[END_REF], [START_REF] Chilali | H∞ Design with pole placement constraints: an LMI Approach[END_REF]:

A T ad P + P A ad -C T ad Y -Y C ad + 2λ min P < 0 (15)
where Y = P L a with P a positive definite matrix, the observer gain L a is deduced as L a = P -1 Y after solving the LMI equations. Using L a as in [START_REF] Chilali | H∞ Design with pole placement constraints: an LMI Approach[END_REF] for the observer gain ensures that all the poles are in the left-half plane region Re(λ) < λ min .

B. H ∞ performance

As stated in the beginning of section III, the road reaction disturbance w should be highly attenuated. In practice, it is sufficient to minimize the disturbance effects on the frequency domain of interest for the considered application. The weighting function described in ( 6) is a filter, where the design parameters are the minimal bandwith ω b , the steady state error k 1 and the maximum of the sensibility function K 1 . In practice, the road disturbance has to be rejected in low frequencies (up to 30Hz), as it is the range of frequencies for driver's torque inputs. The minimization of the disturbances w on the driver torque estimation error z is handled using H ∞ -norm, determined by bounding the transfer function z to w. The problem formulation is given as: minimize γ ∞ such that ||T z w|| ∞ ≤ γ 2 ∞ . Thus the LMI below is satisified:

A T ad P + P A ad -C T ad Y -Y C ad + D T ad D ad P W ad * -γ 2 ∞ < 0 (

C. H 2 performance

The minimization of the effect of sensor noise n on the driver torque estimation error z is handled using the generalized H 2 -norm, of the transfer function from z to n (sse [START_REF] Scherer | Linear matrix inequalities in control[END_REF]). The problem formulation is given as: minimize γ 2 such that ||T zn || ∞ ≤ γ 2 2 . Thus the LMIs below are satisfied:

A T ad P + P A ad -C T ad Y -Y C ad -Y N * -I < 0 P D T ad * γ 2 2 > 0 (17)
It is worth noting that, as in subsection B, a weighting function could be introduced to specify the frequency range on which sensor noise should be attenuated. Besides, sensor noise itself is a high frequency signal.

D. Design procedure in the

H ∞ /H 2 framework
Proposition 1: Consider the system model ( 5) under the assumption ḋ = 0 and the PI observer described in [START_REF] Cholakkal | Fault tolerant control of electric power steering using robust filter-simulation study[END_REF].

If there exists P = P T > 0 and Y satisfying: min αγ ∞ + (1-α)γ 2 with α ∈ [0; 1] with positive scalars γ ∞ , γ 2 subject to the LMIs (15)- [START_REF] Michelis | Method for determining a torque set value for a steering wheel for a power steering system of a motor vehicle[END_REF]. Then, the observer ( 8) is an H ∞ /H 2 PI observer with the gain L a = P -1 Y .

Note that a common Lyapunov matrix P has been used to ensure convexity to the multi-objective problem [START_REF] Chilali | H∞ Design with pole placement constraints: an LMI Approach[END_REF].

IV. DESIGN ANALYSIS

In this part, numerical synthesis results computed from the design described in Section III are presented. First, the results obatined with an observer satisfying the hypothesis of P roposition 1 are shown. Then, a comparison is made between design performances of H ∞ , H 2 and H ∞ /H 2 .

Finding L a such as described in P roposition 1, the following results are deduced: taking α = 0.5, the attenuation level is given by γ ∞ = 0.2365 and γ 2 = 0.645. Moreover, Figure 4 shows the region of the resulting pole placement. To see the efficiency of H ∞ /H 2 , a comparison with H ∞ and H 2 performance has been carried out. From [START_REF] Badawy | Modeling and Analysis of an Electric Power Steering System[END_REF] the effective attenuation subject to minimization of H ∞ -norm is shown by the Bode diagram of the transfer function from the disturbance w to the driver torque estimation error z:

|T zw | < |W -1 f |γ 2 ∞ ( 18
)
Figure 5 shows the resulting attenuation of the disturbance on the estimation error subject to the minimization problems H ∞ , H 2 and H ∞ /H 2 . Note that all methods meet the requirements (18). 5) and ( 8): In Table II

T dd = D ad (sI -(A ad -L a C ad )) -1 L a .C(sI -A) -1 E (19)

V. RESULTS

In this part the experimental results obtained with the designed observer described in Section V are shown.

The PI observer has been tested using experimental data gathered on a development vehicle equipped with a P-EPS system. Besides the usual motor resolver and torque sensor used in P-EPS system, a steering sensor and driver torque sensor have been added to the vehicle. Figure 8 In the following simulation results, the measured steering wheel angle, motor angle and driver torque are sampled every 10ms; the implemented observer is also discretized at the same rate. From figure 9, the frequency response of the continuous and discretized observer are very close, ensuring a good implementation. The experimental data comes from a rolling condition without assist torque (u = 0N.m) on a test track.

A. Test 1: rolling on the track

Figure 10 shows the comparison between the measured driver torque and the estimated torque, when the driver goes around the track. Therefore, vehicle speed and direction are time-varying.

B. Test 2: sinusoidal driving

Figure 11 shows the comparison between the measured driver torque and the estimated torque, in a particular driving condition. It involves following an "8" shaped pattern on the track at a constant speed of 20kph. This test allows to simulate a sinusoidal driver torque. The root mean squared where m is the number of data points. The resulting RMSE for test 1 is 3.21N.m and for test 2 is 3.63N.m. Moreover the normalized RMSE on the whole range is 7% for test 1 and 8.32% for test 2. According to the above simulation results, it can be seen that the observer has a good dynamic transition and an acceptable convergence to the estimated value. Thus performance of the PI observer are validated on experimental data.

VI. CONCLUSIONS

In this paper, an estimation of a P-EPS system driver torque has been proposed using a PI observer. The observer has been designed subject to pole placement and H ∞ /H 2 minimization objectives. Then, the performance of the resulting observer has been evaluated by simulation using experimental measurements.

An industrial interest could be, at first, to develop a limphome operation mode. This means that even if a failure occurs on the torque sensor, an assistance torque (computed from the PI observer) could be provided to the driver. Future works will concentrate on further improvements to the ob-server performance. The road reaction force, e.g represented with Lugre model, should be included in the EPS model [START_REF] Marouf | Driver torque and road reaction force estimation of an electric power assisted steering using sliding mode observer with unknown inputs[END_REF].

Even though the observer can still be improved, it is satisfactory enough to be tested on an actual vehicle using an existing control scheme [START_REF] Michelis | Method for determining a torque set value for a steering wheel for a power steering system of a motor vehicle[END_REF]. Unlike previous field testings, the observer will be evaluated in the presence of an assist torque (u = 0). Depending on the results, a motor model might be added in [START_REF] Marouf | Driver torque and road reaction force estimation of an electric power assisted steering using sliding mode observer with unknown inputs[END_REF] or a new controller based on the PI observer might be designed to get a more efficient observer. Furthermore a robustness analysis of the closed-loop system including the observer and controller could be carried out.

Fig. 1 .

 1 Fig. 1. P-EPS mechanical model

  5) with x = θc θm θp Ẋr θ c θ m θ p X r T ∈ R nx are the internal states, d ∈ R is the driver torque to be estimated, w ∈ R the road reaction torque is the unknown input, and u ∈ R the assist motor torque is the control signal. The available measurements are y = θ c θ m T ∈ R ny which are affected by white gaussian noise n ∈ R. Moreover N is the weighting matrix associated to n in the form N = βI with β ∈ R + .

10

 10

Fig. 2 .

 2 Fig. 2. Frequency domain model validation of T θmu

  ) ∀p ∈ C such that Re(p) ≥ 0. Here n xa = n x + n w and n d denote the number of states x a in (7) and n d = 1 refers to the state d of the driver torque. Defining the extended state x T a d T T and the estimation error e ad = (x a -xa ) T (d -d) T T

Fig. 3 .

 3 Fig. 3. Block diagram of the augmented system and observer

  [START_REF] Scherer | Linear matrix inequalities in control[END_REF]) let * denotes the symmetric element.

Fig. 4 .

 4 Fig. 4. Observer's poles

Fig. 5 .Figure 6

 56 Fig. 5. Comparison between H∞, H 2 and H∞/H 2 design performance on T zw Figure 6 shows the resulting attenuation of the sensor noise on the estimation error subject to the minimization problems H ∞ , H 2 and H ∞ /H 2 . The left side illustrates the attenuation of the driver torque estimation error z depending on the sensor noise on the steering wheel angle y 1 = θ c +βn and the right one depending on the sensor noise on the motor angle y 2 = θ m + βn.

2 Fig. 6 .

 26 Fig. 6. Comparison H∞, H 2 and H∞/H 2 design performance on T zn

2 Fig. 7 .

 27 Fig. 7. Comparison between H∞, H 2 and H∞/H 2 design performance on T dd

Fig. 8 .

 8 Fig. 8. Block diagram of the implemented observer

Fig. 9 .

 9 Fig. 9. Comparison between continuous PI observer and discretized PI observer at 10ms

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Driver torque estimation -Test 1

  the norms of ||T zw || ∞ and ||T zn || 2 obtained from the resulting observers following the three minimization problems H ∞ , H 2 and H ∞ /H 2 are shown.

		H 2	H∞	H∞/H 2
	T zw ∞ T zθc 2 T zθm 2	6.5382E-3 2.8892E-7 4.5015E-3 1.974 9.6728 0.6459 0.5917 0.3028 7.2380E-2

TABLE II COMPUTED

 II NORM OF T zn AND T zn SUBJECT TO STRATEGIES H∞, H 2 Hence the mixed H ∞ /H 2 problem is the best compromise between attenuation of disturbance and sensor noise.

	AND H∞/H 2

ACKNOWLEDGMENT

This work has been supported by JTEKT corporation in Japan: Miyama Shinpukujicho Okazaki-shi, Aichi-ken, 444-2106 Japan.


Wenmeng Xiong

José Picheral

Sylvie Marcos

ANTENNA ARRAY DESIGN FOR MUSIC IN THE PRESENCE OF SPATIALLY DISTRIBUTED SOURCES

Keywords: array signal processing, distributed sources, array geometry, performance, MUSIC

In this paper, the impact of array geometry on the direction of arrival (DOA) estimation of spatially distributed sources impinging on a sensor array is considered. Taking into account the coherently distributed source model proposed in [1], we establish closed-form expressions of the MUSIC-based DOA estimation error as functions of the positions of the array sensors in the presence of model errors due to the angular dispersion of the signal sources. The impact of the array geometry is studied and particular array designs are proposed to make DOA estimation more robust to source dispersion. The analytical results are validated by numerical simulations.

INTRODUCTION

The direction of arrival (DOA) estimation based on snapshots received on a sensor array has been widely studied with plenty of methods [START_REF] Krim | Two decades of array signal processing research: the parametric approach[END_REF]. Among these methods, the multiple signal classification (MUSIC) [START_REF] Stoica | Music, maximum likelihood, and cramer-rao bound[END_REF] is famous for its high resolution in the case of punctual sources. However, in many applications, such as acoustic source imaging [START_REF] Chu | Robust Bayesian super-resolution approach via sparsity enforcing a priori for near-field aeroacoustic source imaging[END_REF], the physical sources can no longer be considered as point sources and a spatially distributed model of the sources could be more appropriate, which degrades the performance of the DOA estimation obtained by MUSIC.

The models for spatially distributed sources have been classified into two types, namely incoherently distributed (ID) sources and coherently distributed (CD) sources. On one hand, for ID sources, signals coming from different points of the same distributed source can be considered uncorrelated. On the other hand, in the scenario of CD sources, the received signal components are delayed and scaled replicas from different points of the same one [START_REF] Valaee | Parametric localization of distributed sources[END_REF]. For CD sources, the performances of MUSIC with discretely distributed sources and continuously distributed sources have been investigated in [START_REF] Astély | The effects of local scattering on direction of arrival estimation with music[END_REF], and [START_REF] Xiong | Performance analysis of MUSIC in the presence of model errors due to the spatial distributions of sources[END_REF], respectively. As expected, the mismatch due to the angular dispersion between the steering vector model of MUSIC and the actual steering vectors of the sources causes estimation errors. MUSIC-based methods for the joint esti-mation of both the DOA and the angular dispersion parameter can be used to avoid the angular dispersion mismatch problem [START_REF] Valaee | Parametric localization of distributed sources[END_REF]. More recently, another MUSIC-based method [START_REF] Abou Chaaya | Localization of spatially distributed near-field sources with unknown angular spread shape[END_REF] makes it possible to have a good performance regardless of the angular distribution shape of the signal sources.

The array geometry on the DOA estimation has been studied in plenty of publications. Optimal array geometries have been designed to reach isotropic and/or optimal performance based on the Cramér-Rao bound (CRB) criterion (eg : [START_REF] Gazzah | Optimum ambiguity-free directional and omnidirectional planar antenna arrays for doa estimation[END_REF][START_REF] Gazzah | Cramer-rao bounds for antenna array design[END_REF]). More recently, based on the spatial aliasing phenomenon, a class of non-uniform array geometries composed of two or more uniform linear arrays (ULAs) with different interelement spacing has been used to reduce the computational burden of Maximum Likelihood (ML) estimator [START_REF] Shin | Reduced-complexity maximum likelihood direction-of-arrival estimation based on spatial aliasing[END_REF].

In this paper, we focus on the impact of the array geometry on the performance of the MUSIC estimator with the CD source model, in the presence of model errors due to the angular dispersion of the source signals. Based on the work in [START_REF] Xiong | Performance analysis of MUSIC in the presence of model errors due to the spatial distributions of sources[END_REF], we first propose an analytical expression of the DOA estimation error as an explicit function of the sensor positions for a single source. Subsequently, the array design of special geometries to ensure robustness to the angular dispersion of the signal sources will be investigated. In the case of two sources, the DOA estimation error with a uniform circular array (UCA) will be proved to be an explicit function of the DOA separation between the two sources. Simulation results illustrate the validation of the theoretical results.

The organization of this paper is as follows. The signal model and a brief recall of MUSIC are given in section 2. The particular array geometry design is investigated in section 3. Numerical simulations are presented in section 4 to validate the analytical results of the previous section. Finally, conclusions are given in section 5.

SIGNAL MODEL AND MUSIC ESTIMATOR

Let us consider q spatially CD far-field narrow-band sources impinging on an array of M sensors. The sources arrive from the DOA θ 1 , ..., θ q , and the position of the m -th sensor is given by the polar coordinate ρ m and α m . Without loss of generality, the signals and the sensors are assumed to be in the same plane, as shown in Figure 1a. The q source signals and the M signals received by the array at moment t are denoted by s(t) = [s 1 (t), . . . , s q (t)] T and y(t) = [y 1 (t), . . . , y M (t)] T , respectively, with:

y(t) = C(θ)s(t) + n(t), (1)
where n(t) ∈ C M ×1 represents the complex Gaussian distributed additive noise, C(θ) = [c h1 (θ 1 ), . . . , c hq (θ q )] ∈ C M ×q is the array steering matrix composed of q steering vectors c hi (θ) that can be written as proposed in [START_REF] Valaee | Parametric localization of distributed sources[END_REF] by:

c hi (θ i ) = π 2 -π 2 a(θ i + φ)h i (φ)dφ, (2)
where i = 1 . . . q, and a(θ) is the steering vector for a point source, which can be given by:

a(θ i ) = e -j2π ρ 1 λ cos(θi-α1) , . . . , e -j2π ρ M λ cos(θi-α M ) T , (3)
where λ is the wavelength, and [•] T is the transpose operation. The function h(φ) is introduced to describe the angular spread distribution (for instance, Uniform and Gaussian distributions) and it can be parameterized by an angular dispersion ∆ which is omitted in the notation. The source signals and the additive noise are considered to be complex centered Gaussian independent random variables. Assuming that signals and noises are uncorrelated and the sources are uncorrelated with each other, the correlation matrix is given by:

R = E[yy H ] = CR s C H + σ 2 b I, (4)
where E[.] is the expectation operator, R s and σ 2 b are the source covariance matrix and the noise variance, respectively.

Under the hypothesis that q < M and R s and C are not rank deficient, it is well known that the decomposition of R into eigenvalues λ m and eigenvectors e m is as follows :

R = M m=1 λ m e m e m = UΛ s U H + σ 2 b VV H , (5)
where U = [e 1 , . . . , e q ] spans the signal subspace defined by the columns of C and V = [e q+1 , . . . , e M ] spans the noise subspace defined as the orthogonal complement of U,

Λ s = diag{λ 1 , . . . , λ q }.
The MUSIC method makes use of the orthogonal property of subspaces spanned by C(θ) and V to estimate the DOA θ. In practice it is difficult to know exactly the angular dispersion of the actual sources, consequently, the steering vector model of the punctual source a(θ) is used here instead of c(θ) to estimate the value of θ so as to have:

θi = argmax θ 1 a H (θ)V 2 . ( 6
)
We note Π = VV H the noise subspace projector corresponding to the actual sources, and Π = I -A(A H A) -1 A H the noise subspace projector corresponding to the model in the estimator, with I the M × M identity matrix, A = [a(θ 1 ), . . . , a(θ q )]. Due to the mismatch of the angular dispersion parameter between the model of the MUSIC estimator a(θ) and the actual source c hi (θ i ), and the fact that the number of snapshots can not be infinite in practice, an estimation error then arises [START_REF] Xiong | Performance analysis of MUSIC in the presence of model errors due to the spatial distributions of sources[END_REF]. In this paper, we focus on the influence of the antenna array geometry on the DOA estimation error due to the model error, the covariance matrix is thus assumed to be perfectly known. In the following, for simplicity we note c i for c hi (θ i ), and a i for a(θ i ), respectively. For source i, from eq(13) in [START_REF] Xiong | Performance analysis of MUSIC in the presence of model errors due to the spatial distributions of sources[END_REF], the DOA estimation error can be given by:

∆θ i = θi -θ i = Re{ ȧH i Πc i } ȧH i Π ȧi . (7)
where ȧi = ∂a ∂θ .

INFLUENCE OF THE ANTENNA ARRAY GEOMETRY

In the following, we show that ∆θ i can be expressed as an explicit function of the antenna positions, the effects of special arrays can thus be investigated.

Single source

In this subsection, we study the particular scenario where there is only one source arriving from the DOA θ. Notice that in this case Π = I -1 M aa H , the DOA estimation error can then be given by:

∆θ = Re{ ȧH c -1 M ȧH aa H c} ȧH ȧ -1 M | ȧH a| 2 . ( 8
)
For the m -th sensor, let us define:

ϕ m = 2π ρ m λ cos(θ -α m ), φm = -2π ρ m λ sin(θ -α m ), u m = cos(ϕ m (θ + φ) -ϕ m (θ))h(φ)dφ, v m = sin(ϕ m (θ + φ) -ϕ m (θ))h(φ)dφ. (9)
Introducing the notations defined in ( 9) in ( 8), ∆θ can be given by:

∆θ = φ • v -φ • v φ2 -( φ) 2 , (10)
where

x = 1 M M m=1
x m is the mean value operator for the variable x. Eq (10) makes it possible to design the geometry in order to minimize ∆θ, Assuming that the angular dispersion of the distributed source is small enough with a symmetrical distribution, we introduce a third order Taylor approximation in φ of ϕ m (θ + φ) -ϕ m (θ) so as:

v m ≈ 1 2 φm σ 2 , ( 11
)
where φm (θ) = -2π ρm λ cos(θ -α m ), and σ 2 φ 2 h(φ)dφ is the parameter related to the angular dispersion of the source signal.

Let us introduce (11) in [START_REF] Shin | Reduced-complexity maximum likelihood direction-of-arrival estimation based on spatial aliasing[END_REF], ∆θ can be given by:

∆θ = 1 2 σ 2 φ • φ -φ • φ φ2 -( φ) 2 . ( 12
)
We can observe that this approximative expression of the DOA estimation error is proportional to the angular dispersion of the signal source.

Uniform circular array

Inspired by (12), we find that for a UCA with an even number of symmetrical elements, where ρ m = ρ, and α m = 2π (m-1)

M

with M even, the terms in the numerator can be null:

φ = - 1 M 2π ρ λ M m=1 sin(θ -α m ) = - 1 M 2π ρ λ   M 2 m=1 sin(θ -α m ) + M 2 m=1 sin(θ -α m + π)   = 0, (13)
and similarly :

φ • φ = (2π ρ λ ) 2 M M m=1 sin(θ -α m ) cos(θ -α m ) 1 2 sin(2θ-2αm) = 0, (14)
which leads ∆θ to be 0, The result illustrates that a symmetrical even elements UCA can be robust to the model error due to the source dispersion, in the scenario of one source.

V-shape array

The V-shape array (VA) is another type of particular geometries which has been widely studied, thanks to its easiness to parameterize. It has been proven in [START_REF] Gazzah | Optimum ambiguity-free directional and omnidirectional planar antenna arrays for doa estimation[END_REF] that with the array angle fixed to a specific value, a VA with any size can have an isotropic behavior, and outperforms the UCA in terms of CRB.

Let us assume that, without loss of generality, the two identical branches are symmetrical with respect to the y axis, with no sensor placed at the origin, as depicted in Figure 1b, and that, the distance between two sensors is fixed to d. The coordinates of the sensors can be expressed as:

ρ m = md, 1 m M 2 (m -M 2 )d, M 2 m M , (15)
and:

α m = α, 1 m M 2 π -α, M 2 m M , (16)
respectively.

Considering the described V-shape geometry, (12) can be derived as:

∆θ = 1 2 σ 2 sin θ cos θ (4 + 4M ) cos 2 α + (2 -M ) sin 2 α (4 + 4M ) cos 2 α sin 2 θ + (2 -M ) sin 2 α cos 2 θ , (17)
we can see that the DOA estimation error is independent of d.

To make ∆θ equal to 0, the numerator of (17) should be 0, which yields:

tan 2 α = 4M + 4 M -2 , (18)
the expression reveals that the array angular α can be fixed to a special value to be robust to the model error due to the source dispersion, which depends roughly on the array size, and is robust to other parameters.

Linear array

Notice that the scenario of α m = 0 in (12) corresponds to the linear array (LA), after some straightforward calculus from (12) with α m = 0, the DOA estimation error is derived as:

∆θ = 1 2 σ 2 cot θ, (19)
which reveals that in the case of LA, the DOA estimation error does not depend on the sensors positions but only on the DOA and the model error. Taking into account that the model error can not be modified, the sensor array should be rotated to reduce the value of cot θ, so as to get a smaller estimation error, if θ is roughly known.

Two sources in the case of a UCA

In the case of multiple sources, we take the example of two sources arriving from θ i , where i = 1, 2 and assuming that θ 2 > θ 1 . For the m -th sensor, let us introduce the notation:

∆ϕ m = ϕ m (θ 2 ) -ϕ m (θ 1 ) = ϕ m,2 -ϕ m,1 . (20)
We here focus on the above-mentioned case of a UCA with an even number of symmetrically positioned sensors, which is easier to manipulate mathematically comparing to other geometries. Notice that with such a geometry:

Re{ ȧH i a i } = - M m=1 2π ρ λ sin(θ -α m ) = 0, (21)
and:

Re{ ȧH i c i } = M m=1 φm,i u m,i =   M/2 m=1 Φ(θ i ) + M/2 m=1 Φ(θ i + π)   =0, ( 22
) where Φ(θ i ) = φm (θ i ) [cos(ϕ m (θ i + φ)) cos(ϕ m (θ i ))+ sin(ϕ m (θ i + φ)) sin(ϕ m (θ i ))] h(φ)dφ.
Introducing ( 21), ( 22), and the expression A = [a 1 , a 2 ] in [START_REF] Abou Chaaya | Localization of spatially distributed near-field sources with unknown angular spread shape[END_REF] gives:

∆θ 1 = Re{-M ȧH 1 a 2 a H 2 c 1 + ȧH 1 a 2 a H 1 c 1 a H 2 a 1 } ȧH 1 ȧ1 (M 2 -|a H 2 a 1 | 2 ) -M | ȧH 1 a 2 | 2 . ( 23
)
Introducing the notations of ( 9) and ( 20) in ( 23). The DOA estimation error as an explicit function of the sensor positions can be given by:

∆θ 1 = φ1 • sin(∆ϕ) cos(∆ϕ) • u -cos(∆ϕ) • u φ2 1 -φ2 1 • cos 2 (∆ϕ) -φ1 sin(∆ϕ) 2 . (24)
Let us note δ = θ 2 -θ 1 , and introduce a third order approximation in δ in (24), after some straightforward calculus, the expression of the DOA estimation error for θ 1 can be derived as:

∆θ 1 = 18 φ2 -3δ 2 φ2 + φ4 φ2 1 • u 1 -φ2 1 • u 1 δ -9( φ2 1 ) 3 + 12 φ2 • φ4 1 + φ2 1 -δ 2 ( φ2 1 + φ4 1 ) 2 .
(25) Let us consider the sign of ∆θ 1 : firstly, Π is Hermitian, so the denominator ȧH i Π ȧi in ( 7) is non-negative, it is enough to focus on the numerator in [START_REF] Abou Chaaya | Localization of spatially distributed near-field sources with unknown angular spread shape[END_REF] for θ 1 which yields δ 3 (18 φ2 -

3δ 2 φ2 + φ4 )( φ2 1 • u 1 -φ2 1 • u 1
) after some straightforward calculus; secondly with a small δ, 18 φ2 -3δ 2 φ2 + φ4 can be approximated by 18 φ2 which is positive; thirdly, by the CauchySchwarz inequality, φ2

1 • u 1 -φ2 1 • u 1 is negative. Therefore, ∆θ 1 is opposite to the sign of δ.
For θ 2 , after some similar calculus, we have:

∆θ 2 = -∆θ 1 . (26)
It is interesting to see that, the value of the estimation error for both the two sources are inversely proportional to the value of δ for small δ, which illustrates that the estimation error decreases as the sources move away from each other. Besides, ∆θ 1 is negative and ∆θ 2 is positive, which implies that when the two sources are close to each other, MUSIC is always able to give the estimation of two DOAs separately, despite that the value of the estimation error increases.

NUMERICAL SIMULATIONS

In this section, we validate the theoretical results through numerical experiments. In all simulations, uniformly distributed sources with the angular dispersion ∆ 0 = 10 • are considered, the theoretical covariance matrix is assumed to be known, with SN R = 10dB. All the arrays are composed of M = 20 sensors, the sensors are uniformly distributed on a circle with the diameter equal to 5λ for the UCA, for the VA, 10 sensors are identically uniformly placed at the two branches, spaced by d = λ/4, and the VA with α = 0 is considered as the case of LA. Such array configurations ensure that the UCA and the LA have the same horizontal aperture. The snapshots number N = 1000 and 100 experiments are realized for the simulation results. In Figure 2, we compare the performance of MUSIC with the UCA and that with the LA versus the DOA separation, in the scenario of two sources. The DOA separation varies with θ c = 1 2 (θ 1 + θ 2 ) fixed to 60 • , assuming that θ 2 > θ 1 . We can see that when the two sources are too close to each other, the estimation error can not be ignored despite the array geometry, due to the interference between the sources. For LA, ∆θ 1 decreases as the DOA separation increases and as θ 1 tends to zero, by contrast, ∆θ 2 decreases as the interference phenomenon fades, but increases again along with the DOA separation, which coincides with the well-known conclusion that the ULA has a better performance for broadside sources. For UCA, both ∆θ 1 and ∆θ 2 tend to zero quickly as the DOA separation increases, which is in adequacy with the scenario of one source where the DOA estimation error is null and outperforms its LA partner. In addition, it is interesting to see the validity of the approximated expression (25) and (26) for describing the trend and the sign of the estimation error when the sources are close to each other.

In Figure 3, the absolute value in logarithm of the estimation error of one uniform distributed source with ∆ 0 = 10 • is plotted with respect to the DOA and the array aperture angle β = π -2α, to validate the approximated expression (17). We note β min for the array aperture angle corresponding to the minimum value of |∆θ|, and β min is plotted versus θ in the dashed line. According to (18), to have a null estimation error, β ≈ 49.7 • regardless of the value of DOA, which corresponds to the simulation result illustrated in Figure 3b.

In Figure 4, the DOA estimation error of MUSIC versus the angular dispersion of the signal source with θ = 45 • are presented, in the case of LA and VA with β = 100 • . We can observe the validation of expression (17) and (19) which makes it possible to consider the DOA estimation error as an explicit convex quadratic function of the angular dispersion of the signal source. Furthermore, regardless of a more complicated geometry, VA outperforms LA, which reveals the interest to optimize the array geometry to have a better performance.

Fig. 1 :

 1 Fig. 1: Illustration of planar arrays

Fig. 2 :Fig. 3 :Fig. 4 :

 234 Fig. 2: DOA estimation error vs. DOA separation (UCA with M = 20, LA with M = 20 and d = λ/4, two sources with uniform angular dispersion, ∆ 0 = 10 • , θ c = 1 2 (θ 1 + θ 2 ) = 60 • )

CONCLUSION

In this paper, we have investigated the impact of the array geometry on the performance of the MUSIC-based DOA estimation in the presence of the spatially distributed sources. We have found that the DOA estimation error can be reduced and even canceled for particular array geometries in the case of one source or in the case where the DOA separation between two sources is large enough. Simulations which are carried out are in adequacy with the proposed theoretical results.


Vincent Strubel

Nicolas Fillot

Fabrice Ville

Philippe Vergne

Alexandre Mondelin

Yves Maheo

Rolling Hybrid

Hybrid Rolling Element Bearings: a New Answer Facing the Problem of Debris Entrapment

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Twin-disc machine

Numerical tool

Context/Objectives Hybrid Rolling Element Bearings: a New Answer Facing the Problem of Debris Entrapment

Vincent STRUBEL [START_REF] Ville | On the Two-Disc Machine: A Polyvalent and Powerful Tool to Study Fundamental and Industrial Problems Related to Elastohydrodynamic Lubrication[END_REF] , Nicolas FILLOT [START_REF] Ville | On the Two-Disc Machine: A Polyvalent and Powerful Tool to Study Fundamental and Industrial Problems Related to Elastohydrodynamic Lubrication[END_REF] , Fabrice VILLE [START_REF] Ville | On the Two-Disc Machine: A Polyvalent and Powerful Tool to Study Fundamental and Industrial Problems Related to Elastohydrodynamic Lubrication[END_REF] , Philippe VERGNE 1 , Alexandre MONDELIN 2 , Yves MAHEO [START_REF] Dwyer-Joyce | The Entrainment of Solid Particles into Rolling Elastohydrodynamic Contacts[END_REF] 1 Université de Lyon, INSA-Lyon, CNRS, LaMCoS, UMR5259, F-69621, Villeurbanne, France 2 SKF Aerospace, F-26300, Châteauneuf-sur-Isère, France

In the aeronautical domain, due to high speeds and weightsaving, hybrid rolling element bearings (REBs) are foreseen to replace classical steel REBs. They offer indeed many advantages including a lower density.

The higher hardness appears also as a preponderant feature to warranty integrity in contaminated lubrication.

Results

Conclusions

Using numerical simulations

Mixing both experimental and numerical works

As particles are dispersed in the lubricant, they can approach the contact region from different positions, named release positions (RP).

Key parameters governing particle entrapment were found : The type of contaminants found in lubricated mechanisms can be very diversified (debris, pollution particles) and from different sources (external or internal).

•
The twin-disc machine [START_REF] Ville | On the Two-Disc Machine: A Polyvalent and Powerful Tool to Study Fundamental and Industrial Problems Related to Elastohydrodynamic Lubrication[END_REF] allows to model an elastohydrodynamic (EHD) point contact.

Considering a contaminated lubricant supposes to take into account both the lubricant behavior and the particles dynamics. Firstly by analyzing streamlines and secondly by determining the theoretical particle trajectories.

The release positions normal to the rolling direction (RP y ) and across the confinement (RP z ) influence the entrapment probability.

The closer to the contact center line and the nearer from the walls the particles are released, the more likely they will be entrapped.

Using experiments

Assumption: each particle is responsible for a single dent  quantitative and qualitative comparison of entrapment phenomena.

Tests with different couples of materials  discs made of silicon nitride withstand indentation by ductile particles (M50 steel).

Spherical

  M50 steel particles (left) are commonly used to study entrapment phenomena. Surface spalling on REBs rings (right) can appear when rolling fatigue occurs in presence of particle dents. LaMCoS, Université de Lyon, CNRS, INSA-Lyon UMR5259, 18-20 rue des Sciences -F69621 Villeurbanne Cedex vincent.strubel@insa-lyon.fr Recurrent problem: debris contamination with destructive behavior

  Dependence of particle entrapment with • Release position • Particle size and mass • Changing contacting materials has a minor influence on particle entrapment phenomena • Silicon nitride discs are able to withstand indentation by ductile particles Simulation of the fluid flow and determination of particle trajectories STLE 70 th Annual Meeting Dallas, May 17-21, 2015

From

  quantitative comparisons: -larger particles are more likely to be entrapped -contacting materials have a little influence on the entrapment ratios … to hybrid REBs "Faster, Higher, Stronger" Larger particles are more likely to be entrapped, in accordance with 2 , due here to a larger entrapment width. Dents Count Tests with the twin-disc machine and the contamination test bench: Controlled level of contamination (nature, size, quantities) Sieved M50-steel powder Silicon nitride and steel discs Pure rolling and fully flooded conditions Funded by the industrial chair: SKF -Fondation INSA de Lyon Numerical simulation : (Contaminants range: 32-40 µm) (Tests performed within steel-steel contacts) Particle size: 40 µm