Didier Lesueur
email: didier.lesueur@lhoist.com

Pierre Métais

Patrick Pibis

Samyr El Bedoui

Hervé Ruat

Sté- Phane Bouron

Ferhat Hammoum

Improving the durability of asphalt mixtures with hydrated lime: Field results from recent French sections

Keywords: Hydrated Lime, Asphalt Mixture, Durability

Although known for a long time, hydrated lime (HL) still attracts a strong technical interest as an asphalt additive. This presentation describes recent field results obtained in France in two different contexts:

-First, we describe the 7-yrs experience obtained on 20 sections with and without HL in the mountainous Vosges department. This eastern part of France is characterized by cold winters and significant elevation, making it a pretty challenging environment for asphalt mixtures. The types of mixtures were essentially semidense Asphalt Concrete (French BBSG). The results after 7-yrs of monitoring are described.

-Second, we describe the 5-yrs experience obtained on 5 sections with and without HL in the oceanic climate of Highway A84 in Normandy. The types of mixtures were also semi-dense Asphalt Concrete (French BBSG) and a thorough testing program was established from the beginning. The results after 5-yrs of monitoring are described. As a result, the results allow quantifying the impact of HL in terms of durability improvement.

Introduction

The use of hydrated lime in asphalt surfacings is a tried and tested practice which developed on a large scale in the late 1970s in the United States, as a response to the widespread and premature appearance of damage such as raveling and aggregate stripping [1]. In Europe the practice is growing, and the Netherlands has the greatest experience, since in this country, Porous Asphalt (PA), the main asphalt type for local motorways, always contains 2% of hydrated lime [1]. Recently, Switzerland and Austria have also been using this additive on a large scale for all types of surfacings [START_REF] Lesueur | Proc. 6th E&E Congress[END_REF].

In France, interest was renewed in 2004 due to the generalization of the use of hydrated lime in surface courses by private highway company Société des Autoroutes du Nord Est de la France (SANEF) [1]. This decision followed the observation that the RILEM International Symposium on Bituminous Materials, ISBM Lyon 2020, Lyon, France, 14 au 16 décembre 2020 first PA laid on its network with hydrated lime in it in 1984, then a first in France, had an exceptional lifetime of over 16 years [START_REF] Sainton | Presented to Club d'échanges "Route et Transport" Île de France[END_REF].

In a similar way, many project owners in various countries across the world now associate the presence of hydrated lime in asphalt mixtures with extended durability of the order of 25%, or greater [1]. This effect is due to the multifunctional effect of hydrated lime in asphalt mixtures, which materializes in an improved resistance to moisture damage, a reduced age-hardening and improved mechanical properties [START_REF] Lesueur | [END_REF].

Vosges department sections

The winters of recent years have caused extensive damage to French road surfaces, in particular with the premature appearance of raveling and aggregate stripping. Damage is generally accentuated by repeated freeze-thaw cycles in the presence of water and de-icing salts. This degradation is aggravated by the phenomenon of expansion of water which, when it freezes, can lead to the delamination of the surface layers [4].

In the recent years, the Vosges Département Road Agency has observed a reduction in the service lives of surface courses, mainly Semi-Coarse Asphaltic Concretes (BBSG -AC surf according to NF EN 13108-1), with disorders such as raveling and aggregate stripping, associated with deficient resistance to moisture damage. They thus decided to put hydrated lime to the test.

Materials

A preliminary study was performed observing that the shorter service lives observed with BBSGs were not associated with a decrease in performance as measured in the Duriez test (NF EN 12697-12 method B) where the loss of compressive strength after total immersion (7 days in water at 18°C) is recorded. It was therefore proposed to make the test more severe by extending it to 14 days, increasing the temperature to 40°C and adding 2% NaCl in the water. With this new set-up, the retained resistance of a reference BBSG dropped from 86% (standard) down to 52% [START_REF] Lesueur | Proc. 6th E&E Congress[END_REF][START_REF] Pibis | [END_REF]. Current specifications being >80%, the mixture was classified as having sufficient resistance with standard procedure but this was not true anymore with the severe one.

In parallel, mixtures with hydrated lime had more than 90% retained resistance for all testing conditions [START_REF] Lesueur | Proc. 6th E&E Congress[END_REF][START_REF] Pibis | [END_REF]. Thus, it was decided to construct 20 test sections in 2011-13 at elevations 300-1090 m, with traffic intensities 35-500 heavy vehicles/day [START_REF] Pibis | [END_REF]. Sections consisted in BBSG mixtures with 5.7-5.8% binder (50/70 pen grade bitumen except for 2 sections in 35/50 and 70/100) with or without up to 1.8% hydrated lime [10]. All mixtures were produced in an industrial batch plant, based on two local aggregate sources (granite and crushed Meuse river gravel).

Monitoring Methodology

For each section, the Mean Texture Depth (MTD -NF EN 13036-1) was measured by the Vosges Department laboratory over time every 20 meters with "axis" (center line)

and "wheel-path" values. 20 points were measured for each lane and the average value was reported.

Results

Average MTD for each section was measured as a function of age with and without lime (Fig. 1). For all new BBSGs, the MTD after construction was equal on average to 0.59 mm. The standard deviation for the 20 values for all the sections was 0.04 mm, showing great uniformity. No significant difference on construction appeared between the formulae with and without lime (Fig. 1).

The differences started to appear from the 3 rd year on (Fig. 1). At this age, the average MTD of the mixtures without lime reached 0.68 mm, whereas that of surfacings with lime stayed close to 0.61 mm. At the same time the standard deviation of the nonmodified surfacings doubled compared to its initial value, indicating a greater heterogeneity of the surface conditions. The field observations enabled this increase to be associated with initial, and scattered, raveling/stripping. The increased standard deviation (0.07 mm) thus illustrated that surface damage was starting to develop.

After the 4 th year, these differences became even more accentuated (Fig. 1) with an average MTD of 0.79 mm for BBSGs without lime, compared to 0.66 mm for the treated-ones. These differences remained similar at the end of the 5 th year (Fig. 1), with a number of sections without lime exceeding 0.8 mm. Degradations was now so marked that some sections required preventive maintenance (microsurfacing or surface dressing). The scatter in MTD was also highly pronounced without lime, with a standard deviation reaching 0.12 mm, well above the repeatability of the method (12%), showing that the degradations were heterogeneously distributed.

After 6 years the sections with lime still had a satisfactory average value of 0.69 mm and a standard deviation of 0.09 mm, i.e. the value observed for BBSGs without lime after 3 years.

3

Highway A84 sections

Materials

37 000 t of AC surf 0/14 (BBSG) were laid in 7-cm-thick layers on highway A84 between the cities of Caen (Normandy) and Rennes (Brittany) in 2012-13. The mixtures were made with siliceous sandstone and 4.4% 35/50 bitumen. 15% Reclaimed Asphalt Pavement (RAP) was introduced in all productions amounting to 0.7% additional binder in the mixture. Hydrated lime was added at 0, 0.5 or 1.5 %. All mixtures were produced in an industrial batch plant. The sections were 1.6 to 6 km-long and where built on both directions of the highway.

Monitoring Methodology

Monitoring was organized since the beginning of the construction works by Cerema Rouen and Ifsttar Nantes. Raw materials and mixtures were tested at the initial stage. Slabs were extracted after 3 and 5 years. They were cut to various dimensions with widths 400-600 mm and lengths 600-700 mm. Specimens from the wearing course could be extracted by further sawing and/or coring. The 35/50 bitumen was tested as such and after short-term (Rolling Thin Film Oven Test RTFOT EN 12607-1) and long term (Pressure Aging Vessel PAV EN 14769) aging. Binders were also recovered (EN 12697-3) from the all mixtures. All binders were tested for standard properties (penetration EN 1426, Ring & ball softening temperature EN1427) and advanced rheological ones (complex shear modulus EN 14770 and Bending Beam Rheometer BBR EN 14771). In parallel, chemical indicators were also measured on the binders such as the asphaltenes content (based on a test method similar to ASTM D6560) and FTIR indices [7].

The mixtures were tested for moisture damage using a test adapted from EN 12697-12, performed on the specimens cored from the slabs, then broken in Indirect Tension in dry conditions and after a conditioning of 14 days at 40°C. Mean Texture Depth (MTD -NF EN 13036-1) were also measured as described in section 2.

Results

The properties measured on the binders are summarized in Fig. 2. First, binder agehardening in the absence of hydrated lime could be observed thru the increase in softening point from an initial value of 54.4°C to 68.4°C after 5 years (Fig. 2). This hardening was also perceived using other indicators, as illustrated by a similar increase in the temperature at which the complex modulus G*/sin delta reaches 1 kPa at 10 rad/s or a decrease in penetration value (Fig. 2).

On the low temperature side, the temperature at which the BBR modulus S(t) after 60s reached 300 MPa (TS=300) or the slope m(t) reached 0.3 S(t) (Tm=0.3), are seen to increase form an initial value of -15.3°C for both parameters to -12.4 and -9.4°C respectively, showing the well-known observation that the cracking temperature increases with binder aging, hence increasing the risk of observing thermal cracking.

Note that the initial binder after RTFOT and RTFOT+PAV reached properties that were quite similar to the ones from the binder extracted right after construction and after 5 yrs respectively (Fig. 2). Fig. 2. Rheological properties for the neat initial bitumen, the same binder after laboratory aging (RTFOT and RTFOT+PAV) and the binder extracted from A84 sections after 5 yrs.

In the presence of hydrated lime and after 5 yrs, the binder showed a reduced aging with a softening temperature of 58.8 or 58°C for respectively 0.5 and 1.5% hydrated lime. These values were very close to that of the binder right after construction, confirming the reduction in age-hardening imparted by hydrated lime [START_REF] Pibis | [END_REF]. The other indicators showed the same fact, including for the low temperature fracture-criteria that remained essentially unchanged. Note that the chemical indicators, based on both asphaltenes and FTIR, confirmed that the reduced aging in the presence of lime was due to a reduction in the asphaltenes creation. Asphaltenes remained close to 16% after 5 yrs for the lime sectors when it was 19after 5 years in the reference one. The slower formation of oxidized compounds such as carbonyls and sulfoxides was also recorded thru FTIR testing.

The results for Indirect Tensile Strength of the mixtures are shown on Erreur ! Source du renvoi introuvable.. After construction, the reference sector showed a pretty high value of 93.8% for the ITSR, that dropped to 49.5% after 5 years (Erreur ! Source du renvoi introuvable.). In the presence of lime, the ITSR could be maintained above 85%. Note that air voids for the 5 years-old specimens with lime were also lower (6-7.4% vs 9.3%), thus participating to the improved resistance to moisture damage.

Conclusions

Mixtures with and without hydrated lime were compared in two different French contexts.

First, 20 sections with semi-coarse ACs (BBSG) in the Vosges department were followed for 6 years based on visual observations and measurements of Mean Texture Depth (MTD). This clearly showed that the mixtures with hyrdated lime expereicned less raveling/stripping and pavements with lime after 6 years were in a similar damagestate as the ones without lime after 3 years. The monitoring will be pursued over time in order to see the performance over a longer time frame, trying also to extract cores and test them in a similar way to what was done on highway A84.

Second, the sections on highway A84, also based on BBSG-type mixtures, showed that hydrated lime slowed-down bitumen aging and allowed the mixture to maintain a good resistance to moisture damage even after 5 years. However, no surface damage was yet visible on the test sections regardless of their lime content, even if some very limited raveling/stripping could be observed on the untreated section, as maybe captured by a slightly higher MTD (Erreur ! Source du renvoi introuvable.). Clearly, monitoring must be continued in order to see if the changes recorded in the indicators materialize into premature damage for some of the sections.

Fig. 1 .

 1 Fig. 1. Mean Texture Depth as a function of pavement age for the 20 sections built in 2011-13 based on BBSG with or without Asphacal® hydrated lime in the Vosges Département.

Fig. 3 .

 3 Fig. 3. Indirect Tensile Strength (ITS) for the initial mixtures and after 5 yrs. The ratio (ITSR) between dry and wet ITS is also given. Void content and Mean Texture Depth are also given.


domain: spi phys

text:
Raphael Levy
email: raphael.levy@onera.fr

Maxime Duquesnoy

Guillaume Aoust Mirsense

Jean Michel

Melkonian Dphy

Myriam Raybaut Dphy

New signal processing for fast and precise QEPAS measurements

Keywords: QEPAS, quartz, resonator, photoacoustic, spectroscopy I

Quartz Enhanced Photoacoustic Spectroscopy (QEPAS) gas sensors have been widely developed over the last decade. This technique takes advantage of a high quality factor tuning fork to enable high sensitivity and high selectivity miniature gas sensors. Lock-in detection is classically used to measure the resonator amplitude which is proportional to the gas concentration, but this technique is slow and leads to measurement drifts as it does not follow the resonator frequency drifts over temperature and pressure. This paper presents a new QEPAS signal processing technique that allows faster and more stable measurements that will enable accurate and fast multi gas sensors.

INTRODUCTION

Photo-acoustic spectroscopy is now a well-established technique in which a gas is detected through the sound it generates by non-radiative relaxation consecutive to absorption of modulated light at a wavelength matching a vibrational transition [START_REF] Bell | Upon the production of sound by radiant energy[END_REF] or a rotational transition [START_REF] Spagnolo | THz Quartz-enhanced photoacoustic sensor for H2S trace gas detection[END_REF]. In traditional photoacoustic spectroscopy the sound pressure is amplified by an acoustic cavity and detected by a microphone. In Quartz-Enhanced Photo-Acoustic Spectroscopy (QEPAS), the sound pressure is generated between the prongs of a quartz tuning fork which actuates it and generates an electrical signal through quartz piezo-electric transduction.

QEPAS sensors have been widely studied since the first publication in 2002 [START_REF] Kosterev | Quartz-enhanced photoacoustic spectroscopy[END_REF], analytical models of the sensor have been developed [START_REF] Petra | Theoretical analysis of a quartz-enhanced photoacoustic spectroscopy sensor[END_REF][START_REF] Aoust | Theoretical analysis of a resonant quartz-enhanced photoacoustic spectroscopy sensor[END_REF] as well as optimized resonant detectors coupled to micro tubes acoustic resonators in different configurations [START_REF] Patimisco | Quartzenhanced photoacoustic spectroscopy: a review[END_REF][START_REF] Aoust | Optimal tuning fork for quartz enhanced photoacoustic spectroscopy[END_REF], but signal processing has not been studied. Usually, lock-in detection, as shown in figure 1, is used to extract the signal at a given frequency from the surrounding noise within a reduced bandwidth. However this detection scheme shows a major drawback that is slow measurement, preventing QEPAS to be used for rapid and multi-gas sensing (measurement time < 100ms).

XXX-X-XXXX-XXXX-X/XX/$XX.00 ©20XX IEEE Indeed, in the frame of lock-in detection, the resonator is not actively controlled at resonance, and as the resonance shifts over temperature, pressure or aging, a sweep of the laser modulation frequency must be performed to detect the frequency that corresponds to the resonator maximum displacement. For each frequency point, the measurement time corresponds to:

𝜏 = 𝑄 𝜋 𝑓 𝑟 , ( 1
)
where Q is the resonator quality factor, and fr its eigen frequency.

This measurement time can be as long as a few tens of seconds as Q must be high (>10 000) and fr low (<35 kHz) to achieve high sensitivity [START_REF] Kosterev | Quartz-enhanced photoacoustic spectroscopy[END_REF]. This long measurement time also decreases accuracy as measurements integrate the signal drifts. Thus it hinders high speed measurement of gas with no idle time.

Only a few other detection schemes were presented to our best knowledge. A technique called electrical modulation cancellation method [START_REF] Zheng | Ppb-level QEPAS NO2 sensor by use of electrical modulation cancellation method with a high power blue LED[END_REF] was used to outbalance the background noise induced by the illumination of the QTF prongs by the laser source. Another method consists in actuating the QTF with optical actuation as well as electrical actuation at the same time [START_REF] Mordmüller | An electrically and optically cooperated QEPAS device for highly integrated gas sensors[END_REF], while the resonator displacement is still measured with a lock-in detector: this enables to keep the resonator at resonance or to measure periodically the resonator Q factor. Yet another technique consists in using a beat frequency to obtain fast measurements, down to the millisecond, of the sensor quality factor, resonance frequency as well as gas concentration [START_REF] Wu | Beat frequency quartz-enhanced photoacoustic spectroscopy for fast and calibration-free continuous tracegas monitoring[END_REF].

In this paper, we propose a closed loop detection scheme enabling fast and precise measurements as well as background noise cancellation. Its theoretical analysis is presented in the following sections as well as first experimental results.

II. THEORETICAL ANALYSIS OF QEPAS CLOSED LOOP DETECTION

A. Principle of operation

The principle of operation shown in Fig. 2 is as follows: the tuning fork is actuated at resonance by means of an oscillator circuit, the oscillator output signal is phase shifted and used to modulate the laser intensity. The laser beam passes in between the tuning fork prongs, generating a photoacoustic signal at the resonator eigen frequency. A phase quadrature is set between the resonator piezoelectric actuation force from the oscillator circuit and the photoacoustic actuation force from the laser beam. This induces a shift of the oscillator frequency that is proportional to the photoacoustic force and thus to the gas concentration. This frequency shift is instantaneous allowing fast measurements.

B. Theoretical model

The closed loop detection scheme is presented in Fig. 2, its equivalent block scheme is shown in Fig. 3. Two feedback branches contribute to actuate the resonator: the oscillator circuit branch that leads to the piezoelectric actuation force Fpz, and the photoacoustic branch that leads to the photoacoustic force Fpa. In order to obtain an oscillator frequency shift depending on the photoacoustic force, these two forces must be in phase quadrature. This is performed by adjusting the phase shift Φ 𝑒𝑙𝑒𝑐 through the phase shifter circuit. Thus, the two modulated forces are in phase quadrature if:

π 2 = Φ pa + Φ propagation + Φ elec .
 Φ pa is the photoacoustic phase shift between the laser modulation and the force on the resonator prongs [START_REF] Kosterev | Quartz-enhanced photoacoustic spectroscopy[END_REF]:

Φ pa = π 2 -𝑎𝑡𝑎𝑛 (2𝜋𝑓 𝑟 𝜏 𝑉𝑇 ), 
where fr is the modulation frequency, τvt the vibrationaltranslational relaxation time of the gas, depending on the wavelength and the gas mixture. Φ propagation is the acoustic propagation delay from the laser beam to the resonator prongs.

Assuming there is no acoustic radiation from the resonator, it can be expressed as:

Φ propagation = 2πL λ 𝑎 , 
where L is the distance from the beam to the resonator branch, and λa the acoustic modulation wavelength. For our custom tuning fork described in section III, this phase shift is: Φ propagation = 0.4 rad.

If the two actuation forces act in phase quadrature, the total force actuating the resonator is:

𝐹 tot = √|𝐹 𝑝𝑧 |² + |𝐹 𝑝𝑎 |² sin (2𝜋𝑓 𝑟 𝑡 + 𝜃)  𝜃 = 𝑎𝑡𝑎𝑛( 𝐹 pa |𝐹 pz | ), 
where  is the additional phase shift in the oscillator loop induced by the two parallel feedback loops (see Fig. 3). This additional phase shift in the oscillator loop induces an oscillator frequency shift:

∆𝑓 𝑟 = 𝜃 2𝑄 𝑓 𝑟 . 
Considering that the photoacoustic force is order of magnitudes smaller than the piezoelectric force, we have:

∆𝑓 𝑟 = 1 2𝑄 𝑓 𝑟 𝐹 pa |𝐹 pz |
. 

C. Differential Measurements enabled

The active detection scheme described above also makes it possible to remove the influence of resonator eigen frequency variations over temperature by periodically adding a π phase shift through the phase shifter circuit. Indeed a π additional phase shift inverts the sign of the acoustic force Fpa applied to the resonator, and thus inverts the sign of the oscillator frequency variation (see Eqs. 6 & 7). Let us write the frequency of the oscillator without supplementary phase shift is:

𝑓 𝑟 | 0 = 𝑓 𝑟0 + ∆𝑓 𝑟 + ∆𝑓 𝑟 (𝑇),
  where fr0 is the initial eigen frequency of the resonator, ∆𝑓 𝑟 is the photoacoustic induced frequency variation, ∆𝑓 𝑟 (𝑇) is the temperature induced frequency variation. The frequency of the oscillator with a π supplementary phase shift is:

  𝑓 𝑟 | 𝜋 = 𝑓 𝑟0 -∆𝑓 𝑟 + ∆𝑓 𝑟 (𝑇) 
The frequency variation over temperature is then cancelled by switching the additional phase shift between 0 and π and subtracting the two measured frequencies:

 𝑓 𝑟 | 0 -𝑓 𝑟 | 𝜋 = 2∆𝑓 𝑟   

D. Limitation of the technique

The addition of a photoacoustic actuation force in phase quadrature with the piezoelectric force causes a phase shift in the oscillator loop. This phase shift must remain low enough for the oscillation conditions to remain valid, i.e. it must remain in the resonator bandwidth:

∆𝑓 𝑟 ≪ 𝑓 𝑟 𝑄 .
   This is the case in the experiment detailed in section III; considering 1% concentration of C2H2 diluted in air and the experimental set-up described; ∆𝑓 𝑟 = 8 mHz, and

𝑓 𝑟 𝑄 = 2.6 Hz
as shown in Fig. 9 and Table 1. Considering [START_REF] Zheng | Ppb-level QEPAS NO2 sensor by use of electrical modulation cancellation method with a high power blue LED[END_REF], this implies that the photoacoustic force is more than two orders of magnitude lower than the piezoelectric force.

E. Noise calculation

In order to compare the resolution of lock-in open loop detection and closed loop detection, the force equivalent noise is calculated for both detection techniques.

1) Lock-in Detection Noise

It has been shown that the ultimate QEPAS noise source originates from the resonator Brownian noise, provided that the system is used near resonance and that the electronical scheme has been correctly designed [START_REF] Petra | Theoretical analysis of a quartz-enhanced photoacoustic spectroscopy sensor[END_REF]. The white Brownian noise power spectral density can be expressed as:

𝐹 open loop = 4𝑘𝑇 𝑀2𝜋𝑓 𝑟 𝑄 (N²/Hz), ( 13
)
where k is the Boltzmann constant, T the temperature, M the equivalent mass of the resonator, fr its eigen frequency, and Q its quality factor.

2) Closed loop detection noise

The main noise source for the closed loop detection originates in the oscillator phase noise. Its power spectral density can be derived from the Leeson equation [START_REF] Rubiola | Phase noise and frequency stability in oscillators[END_REF]:

𝑆 𝜙 (𝛿𝑓) = 𝐹𝑘𝑇 𝑃 0 (1 + 𝑓 𝑟 2 4𝑄 2 𝛿𝑓 2 )(1 + 𝑓 𝑐 𝛿𝑓 ) (rad²/Hz), ()
where 𝛿𝑓 = 𝑓 -𝑓 𝑟 is the carrier offset frequency, F is the noise factor of the oscillator feedback circuit, fc is the flicker phase noise corner of the oscillator feedback circuit, P0 the power through the resonator. As detailed in [START_REF] Levy | Phase noise analysis and performance of the vibrating beam accelerometer[END_REF], the corresponding frequency noise power spectral density can be expressed as:

𝑆 𝐹 (𝛿𝑓) = 𝛿𝑓. 𝑆 𝜙 (𝛿𝑓) (15).
As explained in section II.B, the photoacoustic force is calculated from the eigen frequency variation of the oscillator. In the same way, the force noise power spectral density of the closed loop detection oscillator can be calculated from frequency noise, taking into account [START_REF] Zheng | Ppb-level QEPAS NO2 sensor by use of electrical modulation cancellation method with a high power blue LED[END_REF]:

𝐹 𝑐𝑙𝑜𝑠𝑒𝑑_𝑙𝑜𝑜𝑝 = (𝑆 𝐹 × 𝐹 𝑃𝐴 ∆𝑓 𝑟 ) 2 (16)
Finally with equation ( 16), ( 14) and (15) one can obtain:

𝐹 closed loop (𝛿𝑓) = ( 2𝑄𝐹 pz 𝑓 𝑟 )²𝛿𝑓²(1 + 𝑓 𝑟 2 4𝑄 2 𝛿𝑓 2 )(1 + 𝑓 𝑐 𝛿𝑓 ) 𝐹𝑘𝑇 𝑃 0
(N²/Hz) (17)

3) Noise comparison

In order to compare noise levels of the two detection schemes, the noise power spectral density is computed from ( 13) and ( 16) and shown in Fig. 4. The two configurations are compared: the open loop detection scheme (red line) and the closed loop detection scheme (blue line) working with the same resonator used for experimental results. The parameters of the resonator and the oscillator circuit taken into account are listed in Table 1.

Considering this resonator, the noise of this new technique is worse than the noise of the open loop scheme. This drawback can be overcome by decreasing the piezoelectric actuation force as can be deduced from ( 16). This can be done by adjusting electrodes configuration on the resonator to decrease actuation force and increase piezoelectric detection efficiency.

III. EXPERIMENTAL RESULTS

A. Experimental set-up

The experimental scheme is shown in Fig. 5. The resonator is placed in a hermetic chamber filled with C2H2 diluted at 1% in air at atmospheric pressure and room temperature. The laser source is a telecom laser diode, its wavenumber is set at 6,490.05 cm -1 corresponding to a rovibrational transition of C2H2 as shown in Fig. 6. There is no interference from residual water vapor at this wavelength.

The oscillator circuit is an analog self-sustained oscillator, its output signal controls a waveform generator that enables to control the phase shift of the photoacoustic branch, as well as the waveform and the amplitude of the Mach-Zehnder input for amplitude modulation of the laser radiation. The oscillator frequency shift, proportional to the C2H2 concentration, is measured with a frequency counter (Agilent 53132A). The tuning fork resonator used for the experiment is shown in Fig. 7. It has been designed and realized at ONERA for QEPAS operation [START_REF] Aoust | Ultra high quality factor resonators operated in fluids[END_REF]. In order to increase the eigen frequency variations over photoacoustic actuation force, the piezoelectric force must be kept as small as possible, as shown in [START_REF] Zheng | Ppb-level QEPAS NO2 sensor by use of electrical modulation cancellation method with a high power blue LED[END_REF]. For this purpose only one third of the actuation electrodes have been used for piezoelectric actuation. The parameters of the resonator and the oscillator circuit used for the experiment are shown in Table 1.

B. Experimental results 1) Setting the phase shift

In order to maximize the eigen frequency variation over the photoacoustic actuation force, the phase shift Φ elec must be adjusted with the phase shifter circuit to set a phase quadrature between the photoacoustic force and the piezoelectric force, as shown in [START_REF] Spagnolo | THz Quartz-enhanced photoacoustic sensor for H2S trace gas detection[END_REF]. For that purpose the oscillator frequency variation ∆𝑓 𝑟 is measured over the phase variation of the phase shifter circuit Φ elec as shown in Fig. 8. The phase values corresponding to the minimum and maximum frequency variations are measured:

Φ elec | min = 60° (18)

Φ elec | max = Φ elec | min + 180°= 240° (19)
These values are measured at the beginning of the experiment, and do not need to be measured again, provided that the sound velocity in the gas does not change. 2) Frequency shift measurement Fig. 9 shows the oscillator frequency shift over laser operation. When the laser is on, the frequency increases when the phase shift is set to Φ elec | min (red line), and decreases when the phase shift is set to Φ elec | max (blue line), as described in Eqs. 9 and 10. An additional frequency decrease over time is common to the two oscillator frequencies: it corresponds to the frequency drift with temperature. This temperature sensitivity can be cancelled by computing the difference of the two frequencies, as shown in Fig. 9 and described in [START_REF] Rubiola | Phase noise and frequency stability in oscillators[END_REF].

The measurement time used in Fig. 9 and 10 is 1 s, the measurements being taken alternatively at 0° and 180° phase shifts. 15 measurements are taken alternatively while the laser is off and on. Measurements have also been performed down to 50 ms in order to validate the fast measurements ability of the closed loop detection.

Figure 9: Resonator eigen frequency shift vs. laser on/off. The frequency shift is positive when the additional phase shift is null (red line), whereas it is negative when additional phase shift is 180° (blue line). In addition to the photoacoustic signal frequency sensitivity, the frequency also shifts over time due to temperature variations.

3) Noise equivalent relative absorption measurement

The relative absorption noise can be related to the Allan deviation of the oscillator frequency over measurement time:

𝜎 𝑟𝑎 (𝜏) = 𝐾 𝛼 𝜎 𝑓 (𝜏), (20)
where 𝜎 𝑓 (𝜏) is the oscillator frequency Allan deviation measured with different gate times 𝜏 , α is the absorption coefficient at the line center calculated from HITRAN database and shown in Fig. 5 (α = 1.9 10 -3 cm -1 ), and K is the ratio between the gas absorption and the resulting differential oscillator frequency variation. K has been measured experimentally in the case of the experiment described in section A/ for 1% C2H2 and the oscillator parameters of Table 1, K = 0.24 cm -1 /Hz. Relative absorption noise is plotted in Fig. 11, and measurements with integration time of 100 ms is shown in Fig. 12. The sensitivity of a photoacoustic detector can be qualified by a figure of merit, called normalized noise equivalent absorption (NNEA). It represents the minimum detectable absorption, independently from the laser power, the electrical bandwidth and the absorption coefficient of the gas. It is defined as follows:

𝑁𝑁𝐸𝐴 = 𝑃 𝐿 𝛼 𝑆𝑁𝑅 √𝐵𝑊 , ( 21
)
where PL is the laser power (17 mW), SNR is the signal to noise ratio calculated from relative absorption Allan deviation, and BW is the measurement bandwidth, ie the inverse of the Allan deviation integration time. The calculated NNEA is: NNEA = 3 10 -7 W.cm -1 /√Hz. This value, obtained by this first experiment of QEPAS closed loop detection, is worse than the NNEA obtained with open loop lock-in detection. This is in agreement with theoretical model of noise shown in Fig. 4. The SNR can be improved by improving the electrodes design. Indeed, from [START_REF] Zheng | Ppb-level QEPAS NO2 sensor by use of electrical modulation cancellation method with a high power blue LED[END_REF] it is clear that the piezoelectric actuation force must be as low as possible to increase the measured frequency shift. One way to do this is to make the actuation electrodes very small, enabling to use the resonator surface available to design longer detection electrodes, and then to obtain an improved signal to noise ratio. 1.

IV. CONCLUSION

A new active detection scheme for QEPAS has been described; a theoretical model has been shown and validated by experimental measurements. This new detection scheme enables fast and precise measurements. Fast measurements down to 50 ms integration time have been performed and confirm the precise and thermal drift free measurements. Future work will focus on the signal to noise ratio improvement by optimizing the resonator and its electrodes scheme, and by implementing the detection scheme in synchronous digital electronics on FPGA.

V. ACKNOWLEDGMENT

Figure 1 :

 1 Figure 1: QEPAS lock-in detection scheme.

Fig. 2 .

 2 Fig. 2. QEPAS active detection scheme. The resonator oscillations are sustained by means of an oscillator circuit. The oscillator output signal is phase shifted and controls the modulation of the laser intensity. The measured frequency shift is proportional to the resonator vibration amplitude.

Fig. 3 .

 3 Fig. 3. Closed loop detection scheme. The oscillator feedback circuit branch is presented as well as the photoacoustic feedback. The oscillator frequency shift, proportional to the resonator displacement, is measured by means of the frequency counter.

Figure 5 :

 5 Figure 5: Experimental set-up to test the closed loop detection scheme.

Fig. 6 :

 6 Fig. 6: Calculated absorption of C2H2 at 1 percent concentration at standard atmospheric conditions, using a Voigt broadening model and the HITRAN database

Figure 7 :

 7 Figure 7: Tuning fork resonator used for the experiment.

Figure 8 :

 8 Figure 8: Resonator eigen frequency shift over phase of the phase shifter circuit.

Figure 10 :

 10 Figure 10: Differential eigen frequency shift vs. laser on/off. The common mode frequency variation over temperature has been cancelled by the differential measurement.

Figure 11 :

 11 Figure 11: Measured relative absorption Allan deviation over integration time.

Figure 12 :

 12 Figure12: Differential frequency shift vs. laser on/off for integration time of 1 s (top) and 100 ms (bottom) for 1% C2H2 and the oscillator parameters of Table1.

Table 1 :

 1 Parameters of the resonator and oscillator circuit used for calculations and experimental results. Force equivalent noise power spectral density for open loop detection and closed loop detection, calculated with the parameters of Table1.

	Parameter	Unit	Value
	Q factor		8000
	Motional resistance	MΩ	3
	Resonator eigen frequency (fr)	Hz	21000
	Resonator equivalent mass	kg	1.4 10 -4
	(M)
	Oscillator Noise Factor (F)		1
	Oscillator flicker phase noise	Hz	2
	corner frequency (fc)
	Oscillator power (P0)	W	8.3 10 -10
	Piezoelectric force (Fpz)	N	10 -6
	Photoacoustic force (Fpa)	N	10 -13

The project was partially funded by the FCS Campus Paris-Saclay under grant 2014-0455I (MICSI).


Dr Wadih Ghattas
email: wadih.ghattas@u-psud.fr

Dr Virginie Dubosclard

Dr Sybille Tachon

Dr Morane Beaumet

Prof. Dr Régis Guillot

Dr Marius Réglier

Dr A Jalila Simaan

Jean-Pierre Mahy

] A Jalila Simaan

Jean-Pierre Mahy

J.-P Mahy

CuII -Containing 1-Aminocyclopropane Carboxylic Acid Oxidase Is an Efficient Stereospecific Diels-Alderase

Keywords: Artificial enzymes, Catalysis, Docking, Diels-Alder, ACCO

come

Artificial metalloenzymes (AMEs) can be prepared by incorporating transition metals or their complexes into biomacromolecules. [1] The use of various transition metals should grant AMEs catalytic versatility, whereas biomacromolecules should enable a stereoselective ecofriendly catalysis. Accordingly, AMEs are considered first-class catalysts, particularly in the context of ecofriendly chemistry. [2] Their preparation strategies have evolved to produce in vivo examples, [3] but the "metal substitution strategy" remains the simplest one. This strategy requires stripping a metalloenzyme from its naturally coordinated metal center, followed by the incorporation of a different metal that may provide an artificial reactivity. [4] The substitution of the natural Fe(II) (KD = 4.4 10 -6 M) of 1aminocyclopropane carboxylic acid (ACC) oxidase (Fe(II)-ACCO) by other metal cations such as Ni(II), Co(II), Zn(II) or Cu(II) has been reported. [5] In most cases the exogenous metal cations have been reported to occupy the active site and coordinate the same amino acids as those that originally coordinated Fe(II), namely two histidine and one aspartic acid residues. [5] In particular, Cu(II) binding has been investigated and shown to bind at the active site of the enzyme with an affinity similar to that of the native Fe(II) ion Scheme 1. Proposed mechanism for the catalysis of the oxidation of ACC by natural Fe(II)-ACCO leading to the formation of the phytohormone ethylene. The electrons are provided by ascorbate.

(KD = 5.5 10 -6 M). [5c] While the natural Fe(II)-ACCO catalyzes the oxidation of ACC by dioxygen to form the phytohormone ethylene, all its non-iron reconstituted versions have appeared to be inactive in catalyzing ethylene production. [5b-c] The proposed mechanism of the catalysis proceeds by the coordination of the substrate ACC on Fe(II) in a bidentate fashion, via its amino and carboxylate groups. [6] Subsequently, dioxygen binds to the Fe(II), which leads to the oxidation of ACC and the formation of ethylene, water, cyanide and carbon dioxide (Scheme 1). [7] Inspired by the active site of Fe(II)-ACCO and by the way ACC binds to Fe(II), we decided to investigate Cu(II)-ACCO as an artificial Diels-Alderase. Diels-Alder cycloadditions are important reactions for chemical synthesis. [8] They often take place in organic solvents and require the use of polluting catalysts. Natural enzymes that catalyze these reactions are rare, [9a-f] which has often incited scientists to prepare artificial alternatives that can function under ecofriendly conditions. [10] Known natural Diels-Alderases are metal-free proteins that catalyze the reaction by positioning substrates into the corresponding transition state conformation in what is known as the template effect. [9a-f] Aided by computational design, artificial metal-free Diels-Alderases relying on the template effect have been prepared and genetically evolved to become highly efficient. [9g] Catalytic antibodies that rely on this approach were also developed to catalyze Diels-Alder cycloaddtions. [9h-k] By contrast, Mg(II)-containing macrophomate synthase was suspected to be a metallo-Diels-Alderase that coordinates its dienophile substrate making it electron-deficient and thus more prone to react with the diene substrate. [9l-m] It was further demonstrated, that macrophomate synthase functioned instead by catalyzing a Michael addition and a subsequent aldol addition. [9m-n] Nonetheless, artificial metallo-Diels-Alderases catalyzing Diels-Alder cycloadditions are developed based on the coordination of dienophiles on transition metals using benchmark substrates 2-azachalcon (1a) and cyclopentadiene (2) that provide up to four product isomers (3). It has been shown that Cu(II) is the most efficient water compatible transition metal catalyst for this cycloaddition [11] and therefore most artificial metallo-Diels-Alderases have been prepared by incorporating Cu(II) complexes into biomacromolecules. [12] The mechanism of the reaction has been described to proceed by the bidentate coordination of dienophile 1a on Cu(II) via the nitrogen and the oxygen of its amino and carbonyl groups, respectively. Diene 2 has been shown to be more prone to react with coordinated, hence electron-deficient 1a, as opposed to non-coordinated 1a. [11a-b, 13] Thus, the bidentate binding mode of 1a on copper centers is similar to that of ACC on the active-site metal ion in ACCO. [6a-b, 7] Cu(II)-reconstituted ACCO thus contains the desired transition metal catalyst for the Diels-Alder cycloaddition of 1a and 2 and the required coordination sites to accommodate bidentate 1a are available, along with an easily accessible active site as demonstrated via the binding of diverse substrate analogues and inhibitors (Scheme 2A). [5d] Accordingly, the mechanistic requirements for the catalysis of this Diels-Alder cycloaddition are seemingly fulfilled by Cu(II)-ACCO.

Assays were performed to catalyze the cycloaddition of 1ae and 2, using Cu(II)-ACCO at 2.5 mol% (Scheme 2 A and B, see SI). The obtained results (Table 1) revealed that with substrate 1a, the reaction was quantitative, the endo/exo ratio of 3a isomers was 99/1, and the ee was 95 % in favor of the major endo product. Similar results were obtained with 1b though with a slightly lower yield (95 %), most probably due to its electron donating methoxy substituent as previously observed. [11a-b, 14] Nitro substituted 1c on the other hand, provided the lowest yield (85 %) and ee (54 %) and an endo/exo ratio of 83/17, but was fully converted which indicated that side reactions took place as previously noted with this reactive substrate. [11a-b, 14] It is noteworthy that 1d-e, which are only capable of monodentate coordination on Cu(II) instead of the bidentate coordination that is allowed with reactive 1a-c, did not react. Consequently, it can be assumed that here also the Diels-Alder reaction requires the bidentate coordination of the dienophiles 1a-c on Cu(II). Thereafter, diene 2 reacted with the bound and electron deficient dienophiles. To our knowledge, no other artificial Diels-Alderase was ever prepared by following the "metal substitution strategy", instead, only the "covalent anchoring strategy" and supramolecular insertion of metal cofactors into proteins using either the "hostguest strategy" or the "Trojan-horse strategy" were followed. [3b, 15] The most efficient of these was prepared by the covalent anchoring of 1,10-phenanthroline-Cu(II) complex into the mutant M89C of the lactococcal multidrug resistance regulator. [15d] It provided 3a in up to 93 % yield and 97 % ee with an endo/exo ratio of 95/5. In light of the herein reported results, Cu(II)-ACCO stands out as an artificial Diels-Alderase that is simple to prepare and provides similarly highly stereoselective catalysis but with quantitative yield associated with an endo/exo ratio of 99/1. The crystal structure of natural Fe(II)-ACCO, [5b] revealed a tetrameric form with the C-terminal part of one enzyme monomer interlocking with the C-terminal part of an adjacent monomer. This C-terminal part is thought to be involved in the catalysis as it has been found essential to the natural activity, [16] thus the crystallographic tetrameric form is probably not representing an active conformation of the enzyme. In order to investigate the importance of the C-terminal part in the Diels-Alder activity, we produced an enzyme variant truncated from residue 290, so called ACCOt [16b] and prepared the corresponding Cu(II)-ACCOt AME.We found that unlike Fe(II)-ACCOt that was incapable of ethylene production, both Cu(II)-ACCO and Cu(II)-ACCOt AMEs showed similar catalytic properties (Table 1). These results suggest that the C-terminal loop is neither involved in the coordination of Cu(II) nor in that of substrates 1a-c. Scheme 2. A) Proposed mechanism for the catalysis of the Diels-Alder cycloaddition reaction of substrates 1a-c to 2 by Cu(II)-ACCO leading to the formation of products 3a-c. B) Monodentate substrate analogues.

Table 1. Results of the catalysis of the Diels-Alder cycloaddition of the corresponding substrates 1a-e 2 mM and 2 20 mM after 72 h at 4 ºC.

Catalyst [a] Yield (%) ±3 [b] endo/exo [b] ee (%) ±1 [b]

N Cu II O N O N R R 1a H 1b OMe 1c NO 2 O N O 1d 1e A B R N O N Cu II O N R N O N Cu II N O si face re face O N O N O N O N 3a-c 1a-c 2 re-exo re-endo si-exo si-endo R R R R
ACCOt [16b] and prepared the corresponding Cu(II)-ACCOt AME.We found that unlike Fe(II)-ACCOt that was incapable of ethylene production, both Cu(II)-ACCO and Cu(II)-ACCOt AMEs showed similar catalytic properties (Table 1). These results suggest that the C-terminal loop is neither involved in the coordination of Cu(II) nor in that of substrates 1a-c.

The observed selectivities of AMEs have been investigated with the help of molecular modelling by different groups, including ours. [15a-c] Schrödinger Biologics Suite was used to perform the molecular docking of the substrate 1a in Cu(II)-ACCO (see SI). The results showed that only one pose of 1a in the protein was possible (Figure 1). The Cu(II)-bound carbonyl group and pyridine of substrate 1a occupied a small cavity inside the catalytic pocket. The phenyl group on the other hand, dangled in a larger cavity. The double bond of 1a that is involved in the Diels-Alder reaction was sandwiched between isoleucine 184 from the si-face and lysine 158 from the re-face. Additionally, phenylalanine 250 covered the double bond from the non-coordinated side of 1a. Thus, the si-face of the double bond of 1a stacked on the protein at the inter cavity and was therefore protected from both the siendo and the si-exo attack of 2. The re-face of the double bond was however only protected from the re-exo attack of 2 as it remained accessible from the coordinated side of 1a. According to these results, the 3a si-endo, si-exo and re-exo isomers were all disfavored and the AME should produce selectively the 3a reendo isomer i.e. the (1S,2R,3R,4R) isomer. This is in good agreement with the experimental results as the product (95 % ee) was indeed in the re-endo configuration as determined by circular dichroism, following an already described procedure (see SI). [15c] The minor formation of si-face attack products could possibly be assigned to the reaction taking place in solution i.e. outside the protein catalytic pocket, due to the high concentration of both substrates (5 % yield in buffer in the absence of ACCO and Cu(II)). To verify this hypothesis and in prospect of improving the reaction conditions, we performed kinetic studies.

The kinetics of the catalysis by the AME was explored using 1a and 2 as substrates (See SI). To our knowledge, kinetic parameters have not been determined for other protein-based artificial metallo-Diels-Alderases, most probably because of the relatively high quantities of AMEs that are required for such studies. The products formation was monitored at 4°C as a function of time at different concentrations of 1a while the concentration of the AME was maintained at 50 µM and that of 2 at 20 mM (Figure 2), and initial velocities were extracted (Figure 2 inset). The results revealed that Cu(II)-ACCO follows a pseudo first order Michaelis-Menten kinetics with apparent parameters KM = 230 µM and kapp = 3 h -1 . At Vmax, the reaction catalyzed by the AME was about 100 times faster than that catalyzed by Cu(II) in the absence of protein. [3b, 15b] In comparison, Fe(II)-ACCO was reported to catalyze the oxidation of ACC at 29°C with KM = 150 µM and kcat = 11 min -1 . [16a] Note that under the initially used conditions, 1a and 2 resulted in 5 % yield of all four 3a isomers in the absence of ACCO and Cu(II).

With the kinetic parameters in hand, we attempted to further increase the selectivity of the catalysis by avoiding the formation of products outside the catalytic pocket of the AME. To limit this non-selective reaction of substrates in excess with respect to the AME (2.5 % AME), the addition of both 1a and 2 was fractioned in time. With 1 equiv. AME, the reaction was performed by addition of 5 equiv. 1a and 50 equiv. 2 every 8 h, for up to 40 times followed by an additional 48-50 h then extraction and products analysis. By following this procedure, the excess of substrates was limited. These conditions were very effective as adding up to 90 equiv. of 1a formed exclusively the 3a re-endo product isomer in a quantitative yield (TON = 90, ee > 99 %). Furthermore, the full conversion of up to 185 equiv. of 1a was possible but 3a re-Table 2. Results of the catalysis using 50 mM AME upon fractioned addition of substrates. equiv. of 1a [a] Reaction time for full conversion (h) [b] endo/exo [c] ee (%) ±1 [c] [d] 97/3 92

[a] Fractioned in 5 equiv. and added every 8h with 50 equiv. of 2. [b] Time necessary to achieve full conversion of 1a. [c] All measurements were performed in triplicates (see SI). [d] Full conversion not reached, 96 % yield after 400 h. endo was obtained with 95 % ee. Upon further addition of substrate, full conversion could not be reached most probably because the enzyme was denaturing after more than 16 days at 4 ªC (Table 2, see SI).

In summary, ACCO reconstituted with Cu(II) following the metal substitution strategy was discovered to be an efficient and stereoselective Diels-Alderase. Computational studies were performed and revealed the origin of the high selectivity of the reaction. Additionally, thanks to the simplicity of the preparation of the artificial enzyme, kinetic parameters could be measured, and reaction conditions adapted to perform quantitative and stereospecific Diels-Alder cycloadditions of up to 90 equiv. substrate.

Figure 1 .

 1 Figure 1. substrate 1a is shown as sticks docked in the active site of Cu(II)-ACCO which is shown as a surface colored by heteroatoms. The Cu(II) center is shown as an orange sphere and the relevant amino acid residues as sticks.

Figure 2 .

 2 Figure2. Kinetic of the formation of 3a. The concentration of 2 was 400 times that of the enzyme. Insert, the extracted initial velocities (Vi) as function of substrate concentration. All measurements were repeated at least three times (see SI).

Acknowledgements

S. Tachon was paid by an A*MIDEX grant (n° ANR-11-IDEX-0001-02) funded by the « Investissements d'Avenir » French Government program, managed by the French National Research Agency (ANR). V. Dubosclard and W. Ghattas were financially supported by the grant REBAR (n° ANR-12-PDOC-0007-1) funded by the ANR « Programme Retour Post-Doctorant ».

Entry for the Table of Contents


Razvan Chitic

Frédéric Giraud
email: frederic.giraud@univ-lille1.frwww.univ

Betty Lemaire-Semail

Richard Béarée
email: richard.bearee@ensam.euwww.ensam.fr

Monique Favre
email: mfavre@sinaptec.frwww.sinaptec.fr

Jacky Jehanno
email: jacky.jehanno@melox.frwww.areva.com

Modeling of an Ultrasonic Powder Transportation System

Keywords: «Friction control», «powder transportation», «ultrasound», «vibrations»

This paper presents a new powder transportation system that uses a high frequency flexural stationary wave coupled with a low frequency horizontal displacement of a beam to produce the transport of the powder. The ultrasonic wave is produced with the help of piezoelectric cells glued under the beam and is used to decrease the friction coefficient between the powder and the beam surface.

Introduction

Powder transportation systems can convey powder amount over long distances, but add no value to the product in an industrial process. These systems often use conveyor belts or vibratory conveyors which are often noisy mechanical systems, requiring a massive energy consumption and costly maintenance due to their moving parts. The Ultrasonic Powder Transport Systems use ultrasonic vibrations to transport powder. For example, in [START_REF] Mracek | A system for powder transport based on piezoelectrically excited ultrasonic progressive waves[END_REF] and [START_REF] Kühne | Dosage of highly disperse powders by ultrasound agitated tube modules[END_REF], an ultrasonic progressive wave propagates into a specific pipe material with appropriate absorption properties. So a simple device is designed with no moving part, which can transport powder with good performances. However, the pipe material needs to absorb the reflective wave, dissipating the electrical power into heat. In this paper, we use a high frequency flexural stationary wave coupled with a low frequency horizontal displacement. This principle already demonstrated its capability to convey powder at reasonable speed. The aim of the paper is to focus on the modeling of this device.

Presentation of the system

The system is built up with a steel ruler animated with a low frequency horizontal movement (LFM). If alone, this movement doesn't produce any powder transportation. This is why, piezoelectric cells are glued on the ruler in order to produce a High Frequency Vibration (HFV), by synchronizing HFV on LFM, we can produce the powder transportation. The figure 1 shows the experimental test setup, and the connections between each component. The figure 2 shows the high frequency deformation of the ruler when HFV is switched on. The system works as follows. Powder is deposited on the ruler. Then the linear motor (M) is switched on, producing the LFM. Position of the ruler is measured by a linear position sensor PS. Typical movement produced are oscillations at 1mm of amplitude and 25Hz. According to ruler's position, the High Frequency Amplifier (HFA) is switched ON or OFF. It supplies piezo-cells with an alternative voltage at fixed frequency (around 25kHz) and variable amplitude (0 when HFA is OFF, 200V when HFA is ON). When energized, the piezo-cells produce a stationary wave in the ruler. Its amplitude is about 1µm, and the wavelength is 25mm. Under those high frequency vibrations, the powder is shacked, and lifted up, producing a friction reduction, [3] [4] between the ruler and the powder. Hence, in such contact conditions, the powder can move relatively to the ruler. Performances of the powder transportation depend on the time when HFV is switched on or off. It's in order to determine this moment, that a modeling is introduced, and detailed in next section of this paper.

Modeling of the Ultrasonic Transportation System Effect of the ultrasonic vibrations

Ultrasonic vibrations have several effects on powder, leading to several applications in industry. In this paper, we focus on the friction reduction between the vibrating plate and the powder. To illustrate this point, we inclined the plate, and we put powder on it. Without vibration, the powder stays on the plate, due to friction between powder and the plate. When ultrasonic vibrations are switched on, then the friction decreases, and the powder slips on the vibrating plate almost without inclination of the beam. We then measure the flowing rate of the powder. In this try, there is no LFM, powder is transported only by gravity. As it can be seen, the larger the vibration amplitude is, the lower friction is obtained. However, the behavior depends on powder properties. For example, Alumina has lower friction reduction than Corundum for the same vibration amplitude.

In this paper, we make the assumption that friction has a static behavior, and also depends on powder's speed. We will write:

dt dx f F F v f + = 0 (1)
where F f is the global friction force, F 0 is a static friction force, f v is a viscous friction coefficient, and dx/dt is the powder's speed. Given the results of figure 3, F 0 and f v both depend on W, the amplitude of the HFV.

Modeling behavior

In the final paper, we describe the dynamic behavior modeling of the powder. For that purpose, we assume that the powder stays compact, and can be considered as a solid material [5] [6]. We then consider two separate states:

• The powder slides on the plate is state S1, and the force F f acts on it,

• The powder is static on the plate is state S2. Coulomb's friction law is then laid down on this system: powder changes from state S1 to state S2 when its speed relatively to the plate is null; plate changes from state S2 to S1 when the force of the plate on the powder is above a threshold namely F s . This behavior can be summarized by the chart of figure 4.

Fig. 4: The two states of the powder On a dynamical point of view, there are only 2 forces applied to the powder, which are the force of the plate on the powder and gravity. In our case, we are only interested on the tangential behavior of the powder. Moreover, since dx/dt represents the speed of the powder in the tangential axis, and measured in a moving frame attached to the vibrating plate, we obtained equation (2) to calculate the powder's speed:

) sin( P F ) t sin( MA dt x d M α - - ω ω = 2 2 2 (2)
where M is the powder's mass, A is the LFM 's amplitude, ω is the LFM's pulsation and α is the inclination of the beam. In state S2, the powder doesn't move, so dx/dt=0: F is calculated consequently.

Finally, to take into account HFV, we have two sets of parameters (F 0 , f v and F s ) which are used according to HFV state. To improve the quality of the simulation results, we have tried to measure these parameters by using our test bench. We have measured only the static friction coefficient for the corundum powder in the two cases with HFV ON and with HFV OFF. The coefficients were determined by measuring the angle of beam inclination needed for starting the slipping of the powder on the ruler. The values are presented in the table I.

Table I: Static friction coefficients for HFV ON and OFF

Corundum HFV OFF

HFV ON 0.84 0.14

We can see that, by applying the ultrasonic vibration, the friction coefficient diminish quickly.

Knowing the static friction coefficients, noted µ s , and the gravitational acceleration g, we can calculate the threshold F s using the equation:

Mg F s s µ = (3)
Due to the difficulty of measurement, the viscous friction coefficient was imposed in an arbitrary manner. The results of the figure 5 show the dynamic of the powder calculated on a LFM period. It shows that the speed of the powder measured in the moving reference frame affixed to the plate is not constant. It equals 0 when the powder is in state S2. When HFV are switched on, the powder starts sliding on the plate (at t=0.012sec). When HFV are switched off (at t=0.038), the powder breaks until dx/dt=0. Hence, ultrasonic vibrations produce an asymmetry in contact conditions, which lead to a powder transportation. The dynamic modeling we developed helps to highlight and understand this phenomenon. The use of ultrasound is visible if a comparison between the curve of the absolute velocity of the powder and the speed of the bar is made. It is observed that, at the beginning of the transport cycle, whereas the HFV is OFF, the powder follows the movement of the ruler; then, when HFV is ON, the velocity of the powder follows a different path from that of the bar, which takes up about 0.04s. But, at 0.025s the ruler changes the direction of movement, which does not affect the movement of the powder, thanks to the ultrasonic vibration and also to the kinetic energy that it has accumulated. This simulation allows us to better analyze the influence of the HFV on transport and also to justify the moments when the HFV is activated relative to LFM movement of the ruler.

We can then calculate the average speed of the powder relatively to the plate, and compare the simulation to the measurements. As many parameters are unknown, the approach needs an identification step so as to adapt the value of the parameters of equation ( 1) until simulation and experiment matches as we can see in figure 6. Therefore, to match the two models, we have modified the values of the viscous friction coefficient for the two states: with and without HFV.

Conclusion

In this work, a system of powder transportation is described. It uses a low frequency oscillation movement combined with ultrasonic vibration synchronized on the plate's position. We propose a modeling which allows us to calculate the average speed of the powder in the plate's reference frame. Simulation results are compared to measurements and both are consistent. Future work will focus on the comparison of the system with a more conventional powder transportation system which uses low frequency movements only, in terms of efficiency, reliability and cost.

Fig. 1 : 2 :

 12 Fig. 1: Experimental test setup Fig. 2: The High Frequency Deformation of the Ruler

FlowFig. 3 :

 3 Fig. 3: Powder flow rate for an inclination of 15°

TimeFig. 5 :

 5 Fig. 5: Dynamic behavior of the corundum powder

Fig. 6 :

 6 Fig. 6: Comparison simulation-experiment


Lucie Heim
email: lucie.heim@ensam.eu

Kevin Candelier

Eric Badel

Louis Denaud

Lydie Dufour

Remy Marchal

Rémy Marchal

Hybrid walnut wood quantity

Keywords: Enhancing ecosystem services provision by agroforestry systems Agroforestry, wood quality, growth, wood productivity, Juglans Regia x Nigra

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

This work started in January 2020 and is divided into two main axes: on the one hand, we aim to understand how agroforestry trees adjust to their specific growing conditions, and on the other hand, we aim to study how these growing conditions affects the quality of the wood. As well as the annual wood productivity level, the notion of wood quality, which is very important because it influences the potential valorisation ways of the woods, may allowing to improve their economical values. In our study, the wood quality will be addressed by considering the two following aspects: the mechanical performances and the biochemical composition. Two species are studied: hybrid walnut (Juglans x Regia x Nigra) for the high aesthetic value in the wood market and poplar (Poplar sp) for the packaging and building market. The work is conducted by comparing, for each species, the growth and quality of wood from agroforestry trees with that from trees grown under more conventional conditions in terms of stand densities.

In a first stage, we focused on hybrid walnut species (Juglans x Regia x Nigra) located into the Restinclières Agroforestry Platform (RAP), near Montpellier, France, managed by lab ABsys INRAE (https://umrsystem.cirad.fr/en/the-unit/research-and-training-platform-in-partnership/restinclieres-agroforesteryplatform-rap ). Two plots with 25 years old walnuts are studied: an Agroforestry plot (AF, with 140 walnuts) and a Forestry Control plot (FC, with 235 walnuts).

Agroforestry walnut trees were spaced 4, 8, 12 or 16 m on a same planting line and the density in the plot was 100 trees/ha. Between each line, cereal crops were cultivated. In FC plot, the walnut tree density was almost 200 trees/ha and were mixed with other trees, i.e. alder (Alnus spp) and ash trees (Fraxinus excelsior). Both plots were north-south oriented. In the agroforestry stand, the annual crop was fertilized with approximately 150 kg N ha -1 yr -1 and the soil in the inter-row was usually ploughed to 20 cm every year before the winter crop was sown. The forest control plot did not receive any fertilization.

Dendrometric measurements have been recorded by UMR System since 1995 for the height and 2002 for the circumference and allow the comparison between wood production of agroforestry walnut and forestry walnuts. Our first results show that radial growth of the stem in AF plot was on average 50% greater than in FC plot, with mean values of circumference of 95 cm and 63 cm, respectively (Figure 1 -left). The height of walnut trees in AF plot was also more important than those of trees from FC plot with average values of 12.8 m and 11.5 m high, respectively (Figure 1 -right). Concerning radial growth, our results are in agreement with observations made on the impact of stand density on growth. In forest stands, density has significant effects on radial tree growth and crown size.

Low stand density leads to the formation of trees with larger diameters than trees planted at higher densities [START_REF] Jiang | Effect of stand and tree attributes on growth and wood quality characteristics from a spacing trial with Populus xiaohei[END_REF]). In addition, [START_REF] Cabanettes | Diameter and height growth curves for widely-spaced trees in European agroforestry[END_REF] showed that tree spacing in agroforestry systems favoured radial growth. This could be explained by the higher availability of water, light and enhanced by a higher mechanical stimulation due to the wind [START_REF] Bonnesoeur | Forest trees filter chronic windsignals to acclimate to high winds[END_REF][START_REF] Niez | Water-stressed or not, the mechanical acclimation is a priority requirement for trees[END_REF].

Concerning height growth, our results show the AF walnuts are taller than FC walnut, which differs from other studies indicating that agroforestry trees are usually smaller than forest tree [START_REF] Cabanettes | Diameter and height growth curves for widely-spaced trees in European agroforestry[END_REF]) and that axial growth could be reduced by higher mechanical stimulations due to wind [START_REF] Niez | Water-stressed or not, the mechanical acclimation is a priority requirement for trees[END_REF].

After this first quantitative approach on agroforestry walnut wood, all the other experiments aim at investigating the qualitative aspects of these specific woods; in particular the duraminisation kinetic reactions, the mechanical properties and the peeling aptitude for an industrial target. These experiments are currently in progress.

The continuation of the study on the growth of agroforestry trees will focus on the agroforestry poplar, its mechanical properties and its aptitude for peeling. We will also focused on the impact of wind on the primary and secondary growth dynamics of the trees and on the juvenile/adult transition of the wood whose age of establishment influences the mechanical quality of the wood produced by the tree. In parallel to these two main studies, the effect of wind on tree growth (thigmomorphogenesis) is studied on young Robinia trees (Robinia pseudoacacia) in an experimental agroforestry system in South of France, in order to distinguish the effects of water and light resources from the impact of higher exposure to wind.

Figure 1 .

 1 Figure 1. Left: Circumference at 1,30 m height of AF walnuts and FC walnuts, p-value = 2.2e-16 / Right: Height of AF walnuts and FC walnuts, p-value = 0.00012 -(2017)


K Beddek

S Clénet

O Moreau

Y Le Menach

Solution of Large Stochastic Finite Element Problems -Application to ECT-NDT

Keywords: Finite Element Method, Stochastic Model, Quasi Static Fields, Iterative Solver, Uncertainty Quantification

This paper describes an efficient bloc iterative solver for the Spectral Stochastic Finite Element Method (SSFEM). The SSFEM was widely used to quantify the effect of input data uncertainties on the outputs of finite element models. The bloc iterative solver allows reducing computational cost of the SSFEM. The method is applied on an industrial Non Destructive Testing (NDT) problem. The numerical performances of the method are compared with those of the Non-Intrusive Spectral Projection (NISP).

I. INTRODUCTION

ccounting for inherent uncertainties of numerical model parameters emerges nowadays as an important step to perform robust analysis and reliability assessment. In order to quantify the effect of input data uncertainties on the quantities of interest, numerous techniques were developed. The most popular method is the Monte Carlo Simulation (MCS) method. Nevertheless, as the MCS method requires a large number of deterministic simulations to achieve high accuracy, it becomes unfeasible for very costly numerical simulations. The past decade has seen the rapid development of probabilistic spectral approaches that are based on polynomial chaos representation. The polynomial chaos expansion (PCE) methods have the advantage to exhibit a functional expression of the random quantity of interest. The Galerkin method (spectral stochastic finite element method (SSFEM)) can be applied to solve the stochastic Maxwell equations [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF]. Taking advantage of the Kronecker structure of the linear systems, the storage requirement can be significantly reduced but the computation time of the solution of the linear system is still very high [START_REF] Herzog | Intrusive versus Non-intrusive methods for stochastic finite elements[END_REF]. The non intrusive approaches like the Non Intrusive Spectral Projection (NISP) method are often preferred to SSFEM because they are easier to implement [START_REF] Beddek | Adaptive method for Non-Intrusive Spectral Projection -Application on eddy current NDT[END_REF]. Moreover, even though the linear systems to solve are numerous, they are smaller yielding to a natural parallelization. However, the convergence proof and the error bounds are not so well established as with the SSFEM. In [START_REF] Beddek | Spectral stochastic finite element method for solving 3D stochastic eddy current problems[END_REF], an iterative bloc solver taking advantage of the Kronecker form of the matrix is presented. It enables the reduction of the SSFEM memory requirements. In this paper, we proposed some improvements of the iterative solver. It will be applied to solve a non destructive testing (NDT) problem. The computational costs and results obtained with the SSFEM will be compared to those given by the adaptive NISP presented in [START_REF] Beddek | Adaptive method for Non-Intrusive Spectral Projection -Application on eddy current NDT[END_REF].

II. SPECTRAL STOCHASTIC FINITE ELEMENT METHOD

Consider a conducting domain D with a boundary Γ subdivided into M disjoint subdomains D i on which the permeability and conductivity are supposed to be random but constant and equal to µ i (θ) and σ i (θ) respectively (θ the outcome belonging to the event space Θ ). We suppose that the permeability µ(x,θ) (x denotes the spatial coordinates) and the conductivity σ(x,θ) on D can be written:

( ) ( ) ( ) { } σ µ ≡ κ θ κ = θ κ ∑ = , , , 1 M i D i x I x i (1)
with I Di (x) a function equal to 1 on D i and 0 elsewhere. The magnetic field H(x,θ), the magnetic flux density B(x,θ), the eddy current density J(x,θ) and the electric field E(x,θ) are then random fields. In the frequency domain, these fields verify the following equations:

) , x ( j ) , x ( curl ) x ( ) , x ( ) , x ( θ ω
- = θ + θ = θ B E J J curlH 0 (2)
where ω is the angular velocity and J 0 (x) a deterministic source term. Boundary conditions are added to (2) to insure uniqueness of the solution. As the differential operators (curl, div, grad) operate only on the spatial dimension x, potentials can be introduced, as in the deterministic case, to express the electric and magnetic fields. We consider in this paper the electric potential formulation where the stochastic magnetic potentials A(x,θ) and the stochastic electric scalar potential φ(x,θ) are introduced such that the stochastic magnetoharmonic problem (2) reads:

( ) ( ) ) x ( ) , x ( ) , x ( j ) , x ( ) , x ( ) , x ( 0 J grad A curlA curl = θ ϕ + θ ω θ σ + θ θ µ -1 (3)
We define the functional space H curl (D) and H grad (D) such that:

{ } { } 0 0 2 2 2 2 = ∈ ∈ = = × ∈ ∈ = Γ Γ E E v ), D ( v / ) D ( L v H ), D ( / ) D ( H grad curl L grad n v L curlv L v (4)
We denote by L 2 (Θ) the space of random variables with finite variance. The unknown vector potential A(x,θ) and scalar potential φ(x,θ) are respectively sought in H curl (D)⊗L 2 (Θ) and in H grad (D)⊗L 2 (Θ). A weak form of (3) can be derived by applying the weighted residual for all u belonging to H curl (D)⊗L 2 (Θ):

( )

] ) , ( u ) ( J [ ] ) , ( u ) , ( grad ) , ( A ) , ( [ ] ) , ( curlu ) , ( curlA ) , ( [ 0 γ θ ⋅ = γ θ ⋅ θ ϕ + θ ω θ σ + γ θ ⋅ θ θ µ ∫ ∫ ∫ - d x x E d x x x j x E d x x x E D D D 1 (5)
A where E[X(θ)] denotes the expectation of the random variable X(θ). To solve numerically the problem (5), a discrete form has to be determined on finite dimensional spaces which are subspaces of H curl (D), H grad (D) and L 2 (Θ).

A. Spatial dimension discretization

Usually, with the finite element method, the approximation functions used in computational electromagnetics are the Whitney shape functions [START_REF] Bossavit | Whitney forms: a class of finite elements for threedimensional computations in electromagnetism[END_REF]. We consider a mesh of D with n 0 nodes, n 1 edges, n 2 facets and n 3 elements. The space W 0 of nodal shape functions and the space W 1 of edge shape functions will be used as approximation subspaces of H grad (D) and H curl (D) respectively. In the following, the nodal and edge shape functions are respectively denoted w 0 and w 1 .

B. Random dimension discretisation

In the present work, the stochastic dimension is discretized using the Polynomial Chaos Expansion (PCE). Wiener [START_REF] Wiener | The Homogeneous Chaos[END_REF] was the first to suggest the use of the Hermite chaos to represent Gaussian processes. In [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF], Xiu and al. generalize the concept to more general processes. Let consider a vector ξ ξ ξ ξ(θ)=(ξ 1 (θ),…,ξ M (θ)) of M independent random variables. We denote f ξi the probability density function of ξ 1 (θ). Any output Y of a given model, having as input the random vector ξ ξ ξ ξ(θ), can be written as mapping Y:Θ M →R. The quantity Y is a random variable and we have Y(θ)=Y(ξ ξ ξ ξ(θ)). Let consider ψ j (ξ i (θ)) a univariate orthogonal of order j with respect to the probability measure f [START_REF] Wiener | The Homogeneous Chaos[END_REF] are orthogonal with respect to the joined probability measure

ξi . A multivariate polynomials Ψ α α α α (ξ ξ ξ ξ(θ)) defined as ( ) ( ) ( ) M M i i ,..., with ) ( )) ( ( j α α = θ ξ ψ = θ Ψ = θ Ψ ∏ = α 1 1 α ξ α α
∏ ≤ ≤ = m i i f f 1 ξ ξ That is to say: [ ] β α β α ≠ = θ Ψ θ Ψ if ) ( ) ( E 0 (8)
If Y(θ) has a finite variance, the PCE refers to the representation of the random variable Y(θ) as a linear combination of multivariate polynomials Ψ α α α α (θ):

∑ ∈ θ Ψ = θ m N ) ( y ) ( Y α α α (9)
In practice, the expansion ( 9) is truncated up to the polynomials of order p. If we denote Z M P the space of the Mtuples α α α α which satisfy ||α α α α|| L1 ≤p, then the total number of polynomials in the PCE basis is equal to:

! p ! M )! p M ( P + = (11)
Later, we denote by C P M the space of multivariate polynomials Ψ α α α α (θ) such that α α α α∈ Z M P .

C. Discrete stochastic problem

In the following, we will consider the finite dimensional spaces W 1 ⊗C P M and W 0 ⊗C P M where A(x,θ) and ϕ(x,θ) will be approximated, that is to say:

∑ ∑ ∑ ∑ = ∈ = ∈ θ Ψ ϕ ≈ θ ϕ θ Ψ ≈ θ 0 1 1 0 1 1 n i Z i i n i Z i i M p M p ) ( ) x ( w ) , x ( ) ( ) x ( A ) , x ( α α α α α α w A (12)
Substituting (12) in the weak form [START_REF] Bossavit | Whitney forms: a class of finite elements for threedimensional computations in electromagnetism[END_REF] and applying the n 1 P test functions w 1 i (x)Ψ α α α α (θ) and the n 0 P test functions gradw 0 i (x)Ψ α α α α (θ) leads to a linear equation system:

A s U = B s (11)
The size of the square matrix A s is equal to (NP)x(NP), with N=n 0 +n 1 , the vector U of the unknowns A ij and ϕ ij is of size NP The matrix A s has the following Kronecker structure:

) ( ) ( i i M i i i P s σ σ = µ µ ⊗ + ⊗ + ⊗ = ∑ H S H S A Ι A 1 0 (12)
With I p the identity matrix of size PxP, A 0 is the mean matrix obtained by considering the deterministic magnetoharmonic problem taking the means of the reluctivity and of the conductivity (denoted

     θ Ψ θ Ψ σ - θ σ = θ Ψ θ Ψ µ - θ µ = σ - - µ )] ( ) ( ) ) ( [( E ) ( )] ( ) ( ) ) ( [( E ) ( i i i i i i β α αβ β α αβ S S 1 1 (13)
The matrices H i µ and H i σ are symmetric and of size N×N. The coefficients (H i µ ) jk are given by: [ ]

else if 0 1 2 1 1 1 = ∈ γ ⋅ = ∫ µ n , ) k , j ( if d ) ( k D j jk i i curlw curlw H (14)
And the coefficients (H i σ ) jk by:

[ ] ] ] [ ] ] ] N , n * n , ) k , j ( if d w ) ( N , n ) k , j ( if d w w j ) ( n , ) k , j ( if d j ) ( k D j jk i k D j jk i k D j jk i i i i 1 1 1 0 2 1 0 0 2 1 1 1 1 , , 1 1 ,
∈ γ ⋅ = ∈ γ ⋅ ω = ∈ γ ⋅ ω = ∫ ∫ ∫ σ σ σ w grad H grad grad H w w H (15)
We can notice that the matrix A s is never stored entirely but only the 4M matrices S i µ , S i σ , H i µ and H i σ .

III. ITERATIVE BLOC SOLVER

One way to solve the linear system (11) is to use classical iterative solver and do particular matrix vector product. Another way is to use an iterative bloc solver. Let U 0 an initial solution of the linear system (11). Based on [START_REF] Cung | Efficient numerical strategies for spectral stochastic finite element models[END_REF], a bloc Jacobi iterative solver has been proposed in [START_REF] Beddek | Spectral stochastic finite element method for solving 3D stochastic eddy current problems[END_REF] to solve the linear system (11). The relationship between the solution obtained at the iterations i and i-1 can be written:

1 1 1 0 - - σ σ = µ µ = +         ⊗ + ⊗ - = ⊗ ∑ i s i j j M j j j i P ) ( ) ( T B U H S H S U A Ι (16)
As the matrix (I P ⊗A 0 ) is bloc diagonal of size P, for each iteration i we solve P independent linear systems of size NxN. We can notice that the left hand side matrix is always the same and equal to A 0 . At each iteration i, we solve one linear system of size NxN with multiple and independent right hand sides:

A 0 X=V (17)
where the matrices X and V are of size NxP and defined by:

P 1 1 1 ≤ ≤ ∀ + - = j ), N * j : N * ) j (( ) j (:, U X (18) P 1 1 1 ≤ ≤ ∀ + - = j ), N * j : N * ) j (( ) j (:, T V
The first advantage of this approach lies in a possible parallel computing. At each iteration step, the P linear system resolutions can be performed simultaneously. In our case, a BiCG solver has been used to solve (17). Moreover, since the matrix A 0 is the same for all the linear systems to solve, a unique preconditionner can be applied to reduce the calculation time. The Ainv preconditioner was used whereas in [START_REF] Beddek | Spectral stochastic finite element method for solving 3D stochastic eddy current problems[END_REF] it was a diagonal preconditioner. The residual r i for each iteration i is given by: ( )

1 1 1 0 - = σ σ µ µ = σ σ µ µ -         ⊗ + ⊗ = -         ⊗ + ⊗ + ⊗ = - = ∑ ∑ i i M i i i i i s i M i i i i i i P s i s i ) ( U U H S H S B U H S H S U A I B U A r (19)
By denoting ||X|| 2 , the norm of the vector X, The stopping criterion is then defined by:

2 2 ) ( / s i i B r = ε (20)
With respect to [START_REF] Beddek | Spectral stochastic finite element method for solving 3D stochastic eddy current problems[END_REF], to accelerate the convergence of the iterative procedure, we define a local stopping criterion ε i j for each column of the matrix X such that :

2 0 2 2 1 2 * i i j i i i j X X r r - = ε (21)
With X j i the j-th column of the matrix X (also called j-th mode). This local criterion is such that the modes with high norms (contribute more to the global solution) are well computed and those with low norms, which contribute weakly to the global solution, are less constrained.

IV. APPLICATION AND RESULTS

Our industrial application concerns the detection investigation of clogging in tube support plate in steam generator of nuclear power plant, which may produce performances reduction and safety concern. The clogging detection is performed by an Eddy current NDT technique, called the SAX ratio, using differential bobbin coil signal. In Fig. 1, we give the variation of the flux ∆Φ in function of the position of the probe in the tube. We can see that the variation of ∆Φ is not the same at the upper edge (∆Φ in ) and the lower edge (∆Φ out ). To evaluate the level of clogging, the idea consist in correlating the difference of magnitudes of ∆Φ in and ∆Φ out to the amount of deposit at the inlet foil (Fig. 1). From the control signal magnitude, a clogging SAX ratio indicator is usually extracted as:

in out SAX out Ratio ∆Φ ∆Φ ∆Φ - = (22)
Deterministic simulations are achieved using a finite element model with 1.789.946 spatial unknowns. In the probabilistic model, four parameters are assumed to be random uniform and independent variables: the conductivities and relative permeabilities of the magnetite (clogging product) and the support plate of steam generator tubes. Intervals of variation of the four considered random variables are given in Table 1. These intervals have been determined from expertise and measurements made on samples of magnetite. We fix the polynomial chaos order to p=4, leading to P=70 multidimensional polynomials in the random dimension. The quantities of interest are the flux difference between the two coils of the probe (50 positions of the probe inside the tube) and the Sax ratio. We carry out the simulation using the SSFEM method with the dedicated iterative bloc solver presented in section III and with the Non-Intrusive Spectral Projection (NISP) where a sparse adaptive grid is used to compute the multidimensional integrals [START_REF] Herzog | Intrusive versus Non-intrusive methods for stochastic finite elements[END_REF]. First, we compare performances of the bloc iterative solver (BLOC) to the classical BiCG iterative solver with particular matrix vector product for one position of probe inside the tube (we chose the position 16). With the BiCG solver, the whole linear system (11) is considered during the iterative proces. In the case of the BLOC solver, at each iteration, P linear systems of size NxN are solved (see ( 16)). Two preconditioners are retained for the comparison: Diag for the Jacobi preconditioner and Ainv for the preconditioner based on the approximation of inverse of the matrix A 0 . As reported in Fig. 2, the number of matrix vector products of size N (N is equal to spatial degrees of freedom) of the bloc iterative solver increases linearly with the number of random variables (number of stochastic dimensions), whereas, this number with the BiCG increases exponentially. Notice that the Ainv preconditioner contributes slightly to the acceleration of the two solvers due to singularity and bad conditioning of the linear system resulting from A-ϕ potential formulation. In the second step, the SSFEM and the NISP methods are applied for all probe positions inside the tube. The chosen solver for the SSFEM is the bloc iterative one with global stopping criterion (20) equal to 10 -9 and a local stopping criterion (21) equal to 10 -6 . The stopping criterion of the adaptive anisotropic procedure in the NISP has been chosen equal to 10 -7 and 10 -9 . The computation cost of the SSFEM, reported in Fig. 3, is almost the same for the whole positions of the probe. By contrast, computational costs of the NISP method are higher between the positions 10 and 25. It depends on the position probe, where the variability of the control signal is primarily located [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]. Furthermore, it can be observed that performances of the NISP are closely related to the prior stopping criterion. Both method leads to very close statistical moments on the variation of the flux. The statistical moments of the random SAX ratio obtained with the SSFEM and with a NISP (ε=10 -9 ) are given in Table II. In that case, relative errors of 10%, 30% and 28% are obtained respectively on the standard deviation, the skewness and the kurtosis of the imaginary part of the random SAX ratio. According the size of the problem and the high sensitivity of the ratio SAX, the error between the two methods is acceptable.

V. CONCLUSION

In this paper, the SSFEM is developed in order to propagate uncertainties on the input data through an eddy current problem. A bloc iterative method based solver is proposed to solve the resulting linear system. The SSFEM computational time and the statistical moments of the solution are compared to those given by the anisotropic adaptive NISP method. The large-scale application confirms the success of the proposed iterative solver to reduce the SSFEM computational complexity. The solver can be certainly improved by applying an efficient preconditioner to solve the multiple right hand side linear system.

Fig. 1 .

 1 Fig. 1. Geometry of the studied system and the clogging detection principle.

Fig. 2 .

 2 Fig. 2. Number of matrix vector product of size N (same as deterministic system) of the preconditioned BiCG and Bloc iterative solver for different numbers of random variables

Fig. 3 .

 3 Fig. 3. Number of matrix vector product of size N of the SSFEM and the NISP for each probe position. The SSFEM is equipped with bloc iterative solver (BLOC) and the NISP is equipped with adaptive anisotropic sparse grid with two stopping criteria.TABLE II STATISTICAL MOMENTS OF THE IMAGINARY PART OF THE SAX RATIO,

TABLE II STATISTICAL

 II MOMENTS OF THE IMAGINARY PART OF THE SAX RATIO, COEFFICIENT OF VARIATION AND THE NUMBER OF MATRIX-VECTOR PRODUCT OF THE SSFEM AND THE NISP

		Mean	S-dev.	Skewness Kurtosis CV	Mat-vec
	NISP	0.598	0.093	-0.871	2.761	15%	5.900.000
	SSFEM 0.595	0.103	-1.185	3.539	17%	7.300.000


Guillaume Sahuquet

Ronny Gueguen

Samuel Mer

Adrien Toutant

Françoise Bataille

Gilles Flamant
email: gilles.flamant@promes.cnrs.fr

Particle flow stability in tubular fluidized bed solar receivers

Introduction and experimental setup

Concentrating solar power (CSP) is a promising technology currently undergoing brisk development [START_REF] Barlev | Innovation in concentrated solar power[END_REF]. Among CSP systems, solar towers have the potential to generate power efficiently [START_REF] Behar | A review of studies on central receiver solar thermal power plants[END_REF]. Molten salt is the commonly used heat transfer fluid (HTF) [START_REF] Tyner | eSolar's Modular, Scalable Molten Salt Power Tower Reference Plant Design[END_REF], but solar salt has limited operating temperatures which decrease the efficiency of the power block. To overcome this issue and reduce costs, new particle receivers are under development worldwide as the falling particle receiver of Sandia National Laboratories, the centrifugal receiver of the German Aerospace Center and the fluidized particle-in-tube of the french National Center for Scientific Research [START_REF] Ho | A review of high-temperature particle receivers for concentrating solar power[END_REF]. The European project Next-CSP aims to use solid-gas fluidized bed as HTF in tubular solar receivers. The material chosen is olivine [START_REF] Kang | Particles in a circulation loop for solar energy capture and storage[END_REF] with a Sauter diameter of 30.0 µm (group A of the Geldart classification) that permits to reach high temperatures, up to 800°C and more. This temperature level allows the use of a combined Brayton-Rankine thermodynamic and s-CO2 cycles with a 50 % efficiency. A pilot-scale test unit is currently under construction at Thémis, France, while a smaller cold mock-up was constructed to study the hydrodynamics of the particle suspension flow at PROMES laboratory. Figure 1 shows a schematic (a) and a 3D (b) representation of the cold mock-up. The particles are fluidized in a dispenser thanks to a homogeneously distributed air flowrate. By applying pressure into the dispenser with a pressure-control valve, the fluidized bed passes into a 3m-long transparent glass tube of 4.5 cm internal diameter and is collected in a fluidized bed collector, creating a particle flowrate GP measured by a weighing device. A secondary air flowrate qae is applied at an injection nozzle of diameter dae at the bottom of the tube. Depending on the values of qae and GP, the particle flow is either stable or unstable. When an instability occurs the total pressure in the dispenser decreases and there are significant variations of the particle flowrate. This paper is focused on the conditions of occurrence of this instability that can result, in the solar application, in an emergency stop of the power plant and aims to define a stability diagram of the particle flow in order to define the best operating conditions for the full-scale unit.

Results and discussion

The stability of the particle flow as a function of the secondary air flowrate qae (in Nm 3 /s) and of the solid flowrate GP (in kg/(m².s)) is studied. Firstly, the particle flow is always unstable without aeration. Secondly, two specific values of qae and one caracteristic value of GP delimit two main domains (Fig. 2):

• For values of qae under a threshold air flowrate called qunstable , the mass flow is unconditionally unstable whatever the value of GP . This domain represents the unstable area.

• For values of qae above a threshold air flowrate called qstable , the particle flow is unconditionally stable whatever the value of GP . This domain represents the stable area. For values of qae between qunstable and qstable , the mass flow is either stable or unstable depending on the value of GP. The boundary between these two states follows a sigmoid curve. This domain is called transition area.

The unconditionally stable area is splitted into two sub-areas :

• For small values of GP (under a certain value Gp transition ), where the fluidized bed gets denser with the height in the tube (i.e. the partial derivative of the solid volume fraction with respect to the height is positive);

• For high particle flowrate (GP above Gp transition ), where the fluidized bed dilutes with the height (i.e. the partial derivative of the solid volume fraction with respect to the height is negative). Figure 2.a represents these four domains for a nozzle's diameter of the secondary aeration dae of 0.5 mm. These experimental data show that a secondary aeration at the bottom of the tube is a necessary condition to get a stable particle flow and that the aeration flowrate must overcome a given threshold to reach unconditional particle flow stability. The latter caracteristic is essential for the solar application. In addition, the boundaries of these domains depend on the nozzle's diameter dae of the secondary aeration, with an increase of qstable and qunstable with dae (Fig. 2.b), but with no influence on Gp transition . For our tested diameters, the thinner the nozzle's diameter the smaller the air flowrate necessary to get unconditional stability. In fact, the sigmoid tends to be left shifted with the reduction of the nozzle's diameter. In the final paper, more experimental data would be presented and discussed.

Fig. 1 :

 1 Fig. 1: Schematic representation of the cold mock-up (a), and 3D view (b).

Fig. 2 :

 2 Fig. 2: a) Stability diagram in function of qae and GP. b) Stability diagram in function of dae and qae .

Acknowledgements

Next-CSP is a project that has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 727762.