Thomas H Kauffmann
email: thomas.kauffmann@metz.supelec.fr

Marc D Fontana
email: fontana@metz.supelec.fr

Optical sensor of salt concentration : uncertainty evaluation

Keywords: Raman sensor, salted solution, error analysis, measurement uncertainty

We report on the uncertainty upon the concentration of NaCl in an aqueous solution, as determined from the change in intensity of the OH band within a new Raman optical sensor. The various sources of errors are considered and the standard uncertainty is calculated. The accuracy of the sensor is discussed according to the application.

A c c e p t e d M a n u s c r i p t

Introduction

Raman spectroscopy is a well known optical technique generally used to study the vibrational properties of different media (solid, liquid or gas) in order to get some information on the structure and phase transformations of the sample [START_REF] Colomban | Raman imaging of materials and heterogeneous devices[END_REF][START_REF] Lyon | Raman spectroscopy[END_REF]. In addition, a Raman spectrum can indirectly provide several information about this sample. Thus, the composition of a mixture can be derived from the positions of the different peaks, the mechanical stress from the shift of a particular line, and the partial order or disorder of the structure from the FWHM (Full Width at Half Maximum) of some lines [START_REF] Colomban | Raman imaging of materials and heterogeneous devices[END_REF][START_REF] Lyon | Raman spectroscopy[END_REF].

The Raman spectroscopy is a contactless technique allowing, in situ and fast measurements (sometimes less than one second), so that it can be integrated in sensors to probe some physical properties of a substance [START_REF] Durickovic | Experimental study of NaCl aqueous solutions by Raman spectroscopy : Towards a new optical sensor[END_REF][START_REF] Claverie | Optical Sensor for Characterizing the Phase Transition in Salted Solutions[END_REF].

A c c e p t e d M a n u s c r i p t

The spectroscopic instrument used as a sensor is often composed by many elements (laser source, optical components, spectrometer, camera...) which are all causes of measurement uncertainties. Moreover, the applied method and the signal processing used to extract the information of the signal are also sources of uncertainties.

We have recently conceived a new Raman sensor to detect the sodium chloride (NaCl) in an aqueous solution which is able to determine its concentration [START_REF] Durickovic | Experimental study of NaCl aqueous solutions by Raman spectroscopy : Towards a new optical sensor[END_REF] or its phase [START_REF] Claverie | Optical Sensor for Characterizing the Phase Transition in Salted Solutions[END_REF][START_REF] Durickovic | Water-ice phase transition probed by Raman spectroscopy[END_REF] (solid or liquid). The probe is based on the measurement of the integrated intensity of the Raman stretching OH band in salted solutions. The ratio between two parts of the spectrum gives a parameter, S D , which provides the concentration at a given temperature from an appropriate calibration [START_REF] Claverie | Optical Sensor for Characterizing the Phase Transition in Salted Solutions[END_REF]. The main application devoted to this sensor is the control of the use of NaCl as de-icer by the road winter service, so that these last investigations had specially covered the large concentrations (above 50 g/l). The sensor could be interesting as well to detect low concentrations of NaCl in basins or rivers.

Here, the various sources of the measurement uncertainties in the case of wide (0 up to 200 g/l) and low (0 up to 15 g/l) concentration range are evaluated and discussed.

Sensor description and method

NaCl concentration determination

Our experimental method is based on the detection by Raman spectroscopy of the OH-stretching band of water [START_REF] Carey | Measurement of the Raman spectrum of liquid water[END_REF]. This broad band lying between 2800 and 3860 cm -1 was shown to be affected by the NaCl incorporation in the water [START_REF] Terpstra | Effects of salts on dynamics of water: A Raman spectroscopy study[END_REF]. This change is reflected by an increase of the intensity of the upper part of the band and a simultaneous decrease of the lower part, as shown in Figure 1.

As a consequence, our method to determine the concentration of the NaCl consists into the calculation of the ratio between the integrated intensity of these two areas, J 1 and J 2 , which are affected by the concentration. This parameter S D parameter is therefore defined by:

S D = J 1 J 2 = � 3860 3330 I(w) • dw � 3330 2800 I(w) • dw (1)
where w is the wavenumber and I the intensity. The main advantages of this method considering the ratio of integrated intensities are to avoid the 2 A c c e p t e d M a n u s c r i p t deconvolution of the spectrum, to discard the effect of the spectrometer drift and to reduce the influence of the measurement noise. The signal treatment is trivial and consists in recording directly the Raman spectrum from the CCD camera and in applying several treatments on the useful area (OH-stretching band): baseline correction, intensities normalization and smoothing. Then the two areas, J 1 and J 2 , are calculated in order to get the S D ratio. The last part of the treatment consists in using a calibration S D = f (C) to determine the concentration value C.

Instrumentation and calibration

The main part of our optical sensor (see Figure 2) is composed by a Raman spectrometer including a laser diode (Laser Quantum Ventus) with an exciting line at the wavelength of 532 nm and an output power of 300 mW, an optical head with a Notch filter and a Si CCD camera (Andor technology) of 1024x256 pixels (26 µm 2 ) cooled at -70 • C. The Raman spectrometer (Horiba scientific HE) has a grating of 920 g/mm and a spectral range from 975 to 3870 cm -1 with a spectral resolution of 3 cm -1 . The CCD is linked to the computer by an USB plug and the two optical fibers for the excitation and for the collection have a core diameter of 100 µm. The laser light is focused on the sample through a 50X longworking-distance objective (Olympus) at about 8 mm.

The calibration of the sensor was obtained from the study of the dependence of the parameter S D on the solution concentrations which was determined by the titration method with silver nitrate AgNO 3 .

To achieve the calibration of our probe several measurements were done on NaCl solutions with different concentrations at various temperatures. Indeed the Raman OH-stretching band of the water spectrum is not only affected by the NaCl concentration but also by the temperature, as shown in Figure 3. Two different calibrations of the sensor were done according to the aimed applications. The first calibration concerns the wide concentration range from 0 to 200 g/l with a step of 10 g/l with a temperature excursion from 0 to 10 • C with a step of 0.5 • C. A temperature stage (Linkam THMS600 ) controlled by computer was used to realize this excursion. The second series of measurements were carried out for low NaCl concentrations from 0 to 15 g/l with a step of 1 g/l, from 0 to 30 • C with a step of 0.5 • C.

Each value of S D was obtained from the measurement of the Raman integrated intensities on a salted solution, within constant experimental conditions and the same integrated time (10 accumulations of 1 second in order A c c e p t e d M a n u s c r i p t to reduce the random spectral noise). 21 samples of salted solutions with varying concentration from 0 to 200 g/l were used to get the calibration for the wide range and 16 samples from 0 to 15 g/l were used to achieve the calibration for the low concentration range.

In Figures 4a and4b are plotted the calibrations recorded for different temperatures for the wide and low range of concentration respectively. According to a linear fit applying the least square method, the calibration equation can be expressed for each range by ŜD = a • C + b where ŜD is the predicted value of S D , a is the slope and b the y-intercept. It is observed in Table 1 that the slope a between ŜD and the concentration C is nearly independent of the temperature within each concentration range. By contrast, b strongly varies with the temperature in both concentration ranges.

Measurement uncertainty calculation

At first it is reminded that according to the ISO GUM [START_REF]Guide to the Expression of Uncertainty in Measurement[END_REF], the standard measurement uncertainty can be evaluated via two kinds of methods. The type A evaluation consists into a statistical study from several measurements of the same quantity obtained in assumed constant experimental conditions. The result is chosen as the average of the data and the standard uncertainty is given by the standard deviation of the mean. Uncertainty evaluation by means other than statistical analysis provides the type B contributions.

In our case of the probe of NaCl concentration the different sources of uncertainty are as follows. The calibration gives to rise to a origin of uncertainty linked to the linear regression which is evaluated by a type A method. Then the error in the titration of the sample solution and the error on the measurement of the integrated intensities have to be taken into account. The effect on the determination of the concentration is thus calculated for each source of error. Finally the influence of the temperature on the accuracy of our sensor is considered and evaluated.

In the calculation of the whole uncertainty on the concentration, the squared of standard uncertainties obtained on the various sources of error have to be added. To deduce the standard uncertainty from the experimental uncertainty of type B, a rectangular distribution is conventionally assumed since no indication of the confidence, with which the error was estimated, was supplied by the constructor of the different used instruments.

A c c e p t e d M a n u s c r i p t

Uncertainty related to the calibration

The confidence on the calibration depends on the reliability of the results derived from the fit of the plot S D = f (C) with a regression linear law ŜDi = a • C i + b. For this we apply the type A method and calculate the standard deviation between experimental data of S D and the values as predicted by the law [START_REF] Rosengard | Probabilités et statistiques en recherche scientifique[END_REF][START_REF] Feinberg | Validation interne des méthodes d'analyse[END_REF]:

s(S D ) = 1 √ n -2 • � � � � n � i=1 (S Di -ŜDi ) 2 (2)
where S Di and ŜDi are the experimental and predicted values respectively and n is the number of data used for the calibration. Then we calculate the standard uncertainty on the concentration C coming from the linear regression by the inverse calibration method with the formula [START_REF] Feinberg | Validation interne des méthodes d'analyse[END_REF][START_REF] Lu | Uncertainty evaluation of humidity sensors calibrated by saturated salt solutions[END_REF]:

u calib (C) = s(S D ) a • � 1 N + 1 n + (S D -S D ) 2 a 2 • � n i=1 (C i -C) 2 (3)
where s(S D ) is given by equation ( 2), S D is the value coming from N replicates, S D and C are the averages of the calibration values and C i the different concentration values used in the calibration.

Uncertainty on the titration

The titration of the NaCl solutions used as samples in the calibration procedure is determined with an uncertainty of ΔC = ±0.1 g/l. By applying a rectangular distribution for this type B evaluation, the concentration uncertainty related to the titration is given by u titr (C) = ΔC/ √ 3.

Instrumental uncertainty

There are concerned the uncertainties u instr (S D ) related to the measurement of the S D ratio which depends on the characteristics of the instruments (spectrometer, laser power, ...). The S D ratio is rewritten from equation (1) as:

S D = J 2 J 1 = I 2 • dw 2 I 1 • dw 1 (4)
where I 1 and I 2 are the intensity average of the two areas, respectively J 1 and J 2 , and dw 1 and dw 2 the two Raman ranges used for the areas calculation.

5

A c c e p t e d M a n u s c r i p t

We can therefore deduce the relative standard uncertainty of the S D ratio by:

( u(S D ) S D ) 2 = ( u(I 2 ) I 2 ) 2 + ( u(dw 2 ) dw 2 ) 2 + ( u(I 1 ) I 1 ) 2 + ( u(dw 1 ) dw 1 ) 2 (5)
where u(I) and u(dw) are respectively the standard uncertainty on the intensity I and the wavenumber range dw. The two chosen areas have the same width so that dw 1 = dw 2 = dw and assuming that we have

I 1 ≈ I 2 ≈ I, it comes: ( u(S D ) S D ) 2 = 2 • ( u(I) I ) 2 + 2 • ( u(dw) dw ) 2 (6)
The uncertainty of the intensity I, which is in fact the spectral noise, is estimated at ±0.5% of the spectrum intensity in the optimal condition measurement and the uncertainty of the Raman shift is given by the constructor and is equal to ±1/2 pixel. This value represents for our spectrometer an uncertainty of ±1.5 cm -1 for a band half width of 530 cm -1 , i.e. a relative uncertainty of ±0.3%. The standard uncertainties u(I) of the intensity and u(dw) of the Raman shift are estimated with a rectangular distribution for a type B evaluation.

It is to be mentioned that the drift of the spectrometer causing a shift of the Raman peaks positions is here assumed to be negligible with respect to other errors entering. In fact to diminish this drift a calibration of spectrometer by means of the position peak of Si was made between each series of measurements.

Finally we can find the standard concentration uncertainty related to the instruments by dividing the uncertainty on S D by the slope a of the calibration:

u instr (C) = S D • √ 2 a • � ( 0.005 √ 3 ) 2 + ( 0.003 √ 3 ) 2 (7)

Whole standard uncertainty

The combined standard uncertainty of the concentration can be now obtained by adding the squared values of each uncertainty coming from different error sources, i.e. from the calibration (type A evaluation) and from the titration and the instrumentation (type B evaluations):

u c (C) = � u 2 calib (C) + u 2 titr (C) + u 2 instr (C) (8)
A c c e p t e d M a n u s c r i p t

Temperature influence

We have observed above a strong influence of the temperature on the Raman spectrum of the OH-stretching band so that the defined parameter S D is affected by the NaCl solution temperature as well (see Figures 4a and4b). The temperature therefore can play the role, as frequently in sensors, of an influence external parameter which can impede a correct determination of the quantity to be found. The temperature T is therefore an additional source of error for the concentration determination. For this we establish the relationship between a temperature variation and the measured parameter S D , in our probe when the concentration is assumed to be constant. This leads to perform an additional calibration S D vs T .

Calibration S D = f (T )

In Figure 5 is plotted the S D = f (T ) for three different NaCl concentrations. It is pointed out that even a relatively small temperature change can give rise to a large variation of the parameter S D , used in the determination of the NaCl content with our sensor. Thus a variation of T of about 5 • C causes nearly the same shift of S D as a change of concentration of 10 g/l. This has as first consequence that S D has to be recorded considering both the salt concentration and the temperature. The second effect concerns the additional error on the predicted value of S D , and to the salt content.

The temperature dependences reported in Figure 5 are linear, for each concentration so that Ŝ′ D = a ′ • T + b ′ where Ŝ′ D is the predicted value of S D , a ′ is the slope and b ′ the y-intercept. The linear regression of these plots was obtained by applying the least square method. We calculate the standard deviation s ′ (S D ) between measured and predicted values of S D for p different temperatures by:

s ′ (S D ) = 1 √ p -2 • � � � � p � i=1 (S Di -Ŝ′ Di ) 2 (9)
Thus we can obtain the standard uncertainty u T calib (S D ) of the S D ratio due to the temperature calibration by applying the formula:

u T calib (S D ) = s ′ (S D ) • � 1 N ′ + 1 p + (S D -S D ) 2 a ′2 • � p i=1 (T i -T ) 2 (10)
A c c e p t e d M a n u s c r i p t

where s ′ (S D ) is given by equation ( 9), S D is the value coming from N ′ replicates, S D and T are the averages of the temperature calibration values and T i the different temperature values used in the calibration.

Temperature stage uncertainty

Then the uncertainty of the temperature T is given by the temperature stage constructor as equal to ΔT = ±0.1 • C. We apply therefore a type B evaluation with a rectangular distribution to find u(T ) = ΔT / √ 3. In order to determine the uncertainty on S D caused by the error on the temperature, the standard temperature uncertainty is then multiplied by the slope a ′ of the temperature calibration u T stage (S D ) = a ′ • u(T ).

Concentration uncertainty due to the temperature influence u T (C)

The standard uncertainty on the concentration due to the temperature influence can be evaluated by associating the two above uncertainties of the S D ratio and by dividing the value by the slope a of the concentration calibration:

u T (C) = � u 2 T calib (S D ) + u 2 T stage (S D ) a 2 (11)

Whole concentration uncertainty

At last we combine the uncertainty due to the temperature influence to the uncertainty calculated above arising from the intrinsic sources of errors of type A and B associated to the proposed method:

u c,T (C) = � u 2 c (C) + u 2 T (C) (12)
Finally the expanded uncertainty U c,T (C) of the concentration for a confidence level of 95% is given by U c,T (C) = k • u c,T (C) with a coverage factor of k = 2.

Numerical applications and discussion

Results for the wide range of concentration

Calculations in wide range of concentration (0-200 g/l) are illustrated for a concentration of 50 g/l at a temperature of 10 • C. At this working point S D is equal to 1.3511 and the calibration equations are S D = 0.004380•C +1.1412 for the concentration and S D = 0.009205 • T + 1.2666 for the temperature.

A c c e p t e d M a n u s c r i p t

We calculate first the concentration uncertainty regardless the temperature influence. We get a standard deviation between experimental and predicted data in the concentration calibration equal to s(S D ) = 1.1 • 10 -2 and with N = 10 replicates, n = 21 samples, S D = 1.5792 and C = 100 g/l we find an uncertainty u calib (C) = 1.1 g/l. Then with a titration uncertainty for the NaCl solutions of ΔC = 0.1 g/l we find u titr (C) = 0.1/ √ 3 = 0.06 g/l. The standard uncertainty related to the instruments is found to be u instr (S D ) = 6.4 • 10 -3 so that we get u instr (C) = 1.5 g/l. Finally the combined standard uncertainty of the concentration is found equal to u c (C) = √ 1.1 2 + 0.06 2 + 1.5 2 = 1.9 g/l and with a coverage factor of k = 2 we find an expanded uncertainty of U c (C) = 3.8 g/l.

Then we evaluate the concentration uncertainty due to the temperature influence. We find s ′ (S D ) = 3.6 • 10 -3 as standard deviation between experimental and predicted data. Then with N ′ = 10 replicates, S D = 1.3126 and T = 5 • C we obtain u T calib (S D ) = 1.4 • 10 -3 . From the uncertainty on the temperature stage given by the constructor (±0.1 • C) we get u(T ) = 0.1/ √ 3 = 0.06 • C and then by multiplying by the slope a ′ we obtain u T stage (S D ) = 5.5 • 10 -4 . Finally we calculate the standard uncertainty on the concentration due to the temperature by associating the two above uncertainties and by dividing the value by the slope a of the concentration calibration:

u T (C) = � (1.4 • 10 -3 ) 2 + (5.5 • 10 -4 ) 2 0.004380 2 = 0.3g/l. ( 13
)
At last we combine the uncertainty due to the temperature influence with the uncertainty coming from the other sources (calibration, titration, instruments) to get the total concentration uncertainty u c,T (C) = √ 1.9 2 + 0.3 2 = 1.9 g/l. Finally the expanded uncertainty of the concentration with a confidence level of 95% is U c,T (C) = 2 * 1.9 = 3.8 g/l for C = 50 g/l and T = 10 • C.

Results for the low range of concentration

To illustrate the case of the low range of concentration (0-15 g/l), calculations of the uncertainty are given for a concentration of 10 g/l at a temperature of 20 • C. At this working point S D is equal to 1.3042 and the calibration equations are S D = 0.003602 • C + 1.2682 for the concentration and S D = 0.008402 • T + 1.1401 for the temperature.

A c c e p t e d M a n u s c r i p t

We proceed with the same way as for the wide concentration range. Thus we find a standard deviation between experimental and predicted data coming from the procedure equal to s(S D ) = 2.3 • 10 -3 and with N = 10 replicates, n = 16 samples, S D = 1.2952 and C = 7.5 g/l we get an uncertainty u calib (C) = 0.3 g/l. The titration uncertainty is unchanged and equal to u titr (C) = 0.06 g/l. The standard uncertainty related to the instruments is u instr (S D ) = 6.2 • 10 -3 giving rise to u instr (C) = 1.7 g/l. Finally the combined standard uncertainty of the concentration is found equal to u c (C) = √ 0.3 2 + 0.06 2 + 1.7 2 = 1.7 g/l and with a coverage factor of k = 2 we find an expanded uncertainty of U c (C) = 3.4 g/l.

The temperature calibration was made for the low range between 0 • C and 30 • C by a step of 0.5 • C i.e. p = 61 samples. We find s(S D ) = 5.4 • 10 -3 as standard deviation between experimental and predicted data within the temperature calibration and with N ′ = 10 replicates, S D = 1.2661 and T = 15 • C we obtain u T calib (S D ) = 1.9 • 10 -3 . Then with the temperature stage uncertainty given by the constructor we find u T stage (S D ) = 0.008402 * 0.06 = 5 • 10 -4 . Finally we calculate the standard uncertainty on the concentration due to the temperature by:

u T (C) = � (1.9 • 10 -3 ) 2 + (5 • 10 -4 ) 2 0.003602 2 = 0.5g/l. ( 14
)
We note that the contribution related to the temperature stage uncertainty is negligible (u T stage (C) = 0.1 g/l) compared to the contribution coming from the calibration (u T calib (C) = 0.5 g/l). At last after associating the various sources of uncertainty we get the total concentration uncertainty u c,T (C) = √ 1.7 2 + 0.5 2 = 1.8 g/l and U c,T (C) = 3.6 g/l for the low concentration 10 g/l at 20 • C.

Discussion

All the identified sources of uncertainty on the concentration as well as their relative contribution are summarized in Table 2 for the wide and low concentration ranges. They are defined by the ratio of the squared of considered standard uncertainties.

At first we note that the expanded uncertainty has nearly the same value for both concentration ranges. This gives rise to a good relative accuracy (U/C < 10%) when probing large NaCl concentration (50 g/l) but a poor accuracy if our sensor is devoted to measure low concentration (U/C ∼ 25% if C = 15 g/l).

A c c e p t e d M a n u s c r i p t

If we look at the intrinsic origins of the uncertainty, it can be pointed out from Table 2 that the instrumental error has nearly the same value and provides the largest contribution for both concentration ranges. On the contrary, the uncertainty coming from the calibration procedure is large in the wide concentration range but negligible for low concentration measurements. This discrepancy arises from the much higher density of data (number of data within the concentration range) used in the last case.

The second remark concerns the temperature as an external parameter influencing the measurement and the sensor accuracy. Table 2 shows that the whole uncertainty has the same or nearly the same value if regarding or not the temperature effect for both concentration ranges. This observation does not mean that the temperature does not play an important role in the collecting data procedure with our sensor. In fact the reliability of the measurements depends on the knowledge of the calibration at the correct temperature. In other terms the non-accounting of the temperature should lead from the measured value of S D to an erroneous determination of the concentration.

Our estimation of the different sources of errors to be considered within the sensor proposed to determine the salt concentration, provides a too large uncertainty if probing small values of concentration (36% for a concentration of 10 g/l). By contrast, the results should be more reliable if measuring concentration values larger than 50 g/l. In this case the relative accuracy should be reduced to 10% or less (3.8% for a concentration of 100 g/l). This better accuracy comes from the fact that the main error contribution is independent of the concentration values. This renders our optical sensor more efficient for applications for which large values of concentration are purchased [START_REF] Marchetti | Détermination de la quantité résiduelle de fondants routiers par une technique spectroscopique sans contact[END_REF]. In fact even for low concentrations, the accuracy of our sensor can be improved, regarding and analyzing results reported in Table 2.

The only source, on which we can efficiently act to reduce the total uncertainty of our sensor, is related to the instruments (spectrometer and detector). It is indeed possible to reduce the spectral noise by increasing the number of accumulations and collecting a larger signal in using a more powerful laser and increasing the exposure time.

Additionally CCD detector used in the calibration can be replaced by another one with more pixels (2048 instead of 1024 for example). In this case the resolution will be twice better with the same spectrometer grating (±0.75 cm -1 for the same band width of 530 cm -1 ) and the uncertainty due to the instruments can be recalculated to find u instr (C) = 1.5 g/l instead of A c c e p t e d M a n u s c r i p t 1.7 g/l and a total concentration uncertainty u c,T (C) = 1.6 g/l instead of 1.8 g/l for the low concentration range.

The useful Raman signal covers a spectral range of about 1060 cm -1 from 2800 cm -1 to 3860 cm -1 while the full spectral range with the current grating starts at 975 cm -1 . So a huge part of the range is useless in our method and the grating can therefore be replaced by another one with more grooves per millimeter in order to use all the CCD dimension and thus improve the spectral resolution.

Conclusion

In this study, was reported and discussed an estimation of uncertainty to be considered in a new Raman sensor to probe the NaCl concentration in a salted solution. This method provides quick, in situ and remote measurements. The various sources of uncertainties were taken into account and calculated. In particular, the influence of the temperature change during probing was pointed out and discussed. Based on the analysis of estimation of expanded uncertainty at 95% confidence level with a coverage factor of k = 2, the values of concentration measurement uncertainty were evaluated at 3.8 g/l for the wide range of concentration (0-200 g/l) and at 3.6 g/l for the low range (0-15 g/l). The concentration uncertainty is therefore nearly the same for the two ranges so that the sensor is better suitable for the high concentration applications.

The main application of this sensor is the road winter service with the measurement of high concentrations but other applications using small concentrations could be exploited by our sensor. In these cases a better accuracy is required and can be provided by a more appropriate choice and performances of instruments. Table 1. Slope and intercept of the calibration linear fit for several temperatures and for the two concentration ranges. Table 2. Uncertainties calculated for 50 g/l at 10

Figure 1 .

 1 Figure 1. Evolution of the Raman OH-stretching band with the NaCl concentration. In the insert is reported the whole OH band spectrum.Figure 2. Experimental setup description.Figure 3. Influence of the temperature on the Raman OH-stretching band of a NaCl solution at 10 g/l.Figure 4a. Concentration calibration lines from 0 to 200 g/l (wide range) at different temperatures.Figure 4b. Concentration calibration lines from 0 to 15 g/l (low range) at different temperatures. Figure 5. Temperature dependence of the parameter S D for different concentrations.Table1. Slope and intercept of the calibration linear fit for several temperatures and for the two concentration ranges. Table2. Uncertainties calculated for 50 g/l at 10 • C (wide concentration

Figure 2 .

 2 Figure 1. Evolution of the Raman OH-stretching band with the NaCl concentration. In the insert is reported the whole OH band spectrum.Figure 2. Experimental setup description.Figure 3. Influence of the temperature on the Raman OH-stretching band of a NaCl solution at 10 g/l.Figure 4a. Concentration calibration lines from 0 to 200 g/l (wide range) at different temperatures.Figure 4b. Concentration calibration lines from 0 to 15 g/l (low range) at different temperatures. Figure 5. Temperature dependence of the parameter S D for different concentrations.Table1. Slope and intercept of the calibration linear fit for several temperatures and for the two concentration ranges. Table2. Uncertainties calculated for 50 g/l at 10 • C (wide concentration

Figure 3 .

 3 Figure 1. Evolution of the Raman OH-stretching band with the NaCl concentration. In the insert is reported the whole OH band spectrum.Figure 2. Experimental setup description.Figure 3. Influence of the temperature on the Raman OH-stretching band of a NaCl solution at 10 g/l.Figure 4a. Concentration calibration lines from 0 to 200 g/l (wide range) at different temperatures.Figure 4b. Concentration calibration lines from 0 to 15 g/l (low range) at different temperatures. Figure 5. Temperature dependence of the parameter S D for different concentrations.Table1. Slope and intercept of the calibration linear fit for several temperatures and for the two concentration ranges. Table2. Uncertainties calculated for 50 g/l at 10 • C (wide concentration

Figure 4a .

 4a Figure 1. Evolution of the Raman OH-stretching band with the NaCl concentration. In the insert is reported the whole OH band spectrum.Figure 2. Experimental setup description.Figure 3. Influence of the temperature on the Raman OH-stretching band of a NaCl solution at 10 g/l.Figure 4a. Concentration calibration lines from 0 to 200 g/l (wide range) at different temperatures.Figure 4b. Concentration calibration lines from 0 to 15 g/l (low range) at different temperatures. Figure 5. Temperature dependence of the parameter S D for different concentrations.Table1. Slope and intercept of the calibration linear fit for several temperatures and for the two concentration ranges. Table2. Uncertainties calculated for 50 g/l at 10 • C (wide concentration

Figure 4b .

 4b Figure 1. Evolution of the Raman OH-stretching band with the NaCl concentration. In the insert is reported the whole OH band spectrum.Figure 2. Experimental setup description.Figure 3. Influence of the temperature on the Raman OH-stretching band of a NaCl solution at 10 g/l.Figure 4a. Concentration calibration lines from 0 to 200 g/l (wide range) at different temperatures.Figure 4b. Concentration calibration lines from 0 to 15 g/l (low range) at different temperatures. Figure 5. Temperature dependence of the parameter S D for different concentrations.Table1. Slope and intercept of the calibration linear fit for several temperatures and for the two concentration ranges. Table2. Uncertainties calculated for 50 g/l at 10 • C (wide concentration

Figure 5 .

 5 Figure 1. Evolution of the Raman OH-stretching band with the NaCl concentration. In the insert is reported the whole OH band spectrum.Figure 2. Experimental setup description.Figure 3. Influence of the temperature on the Raman OH-stretching band of a NaCl solution at 10 g/l.Figure 4a. Concentration calibration lines from 0 to 200 g/l (wide range) at different temperatures.Figure 4b. Concentration calibration lines from 0 to 15 g/l (low range) at different temperatures. Figure 5. Temperature dependence of the parameter S D for different concentrations.Table1. Slope and intercept of the calibration linear fit for several temperatures and for the two concentration ranges. Table2. Uncertainties calculated for 50 g/l at 10 • C (wide concentration

  Figure 1. Evolution of the Raman OH-stretching band with the NaCl concentration. In the insert is reported the whole OH band spectrum.Figure 2. Experimental setup description.Figure 3. Influence of the temperature on the Raman OH-stretching band of a NaCl solution at 10 g/l.Figure 4a. Concentration calibration lines from 0 to 200 g/l (wide range) at different temperatures.Figure 4b. Concentration calibration lines from 0 to 15 g/l (low range) at different temperatures. Figure 5. Temperature dependence of the parameter S D for different concentrations.Table1. Slope and intercept of the calibration linear fit for several temperatures and for the two concentration ranges. Table2. Uncertainties calculated for 50 g/l at 10 • C (wide concentration

Figure1A

  Figure1

Figure3A

  Figure3

Figure4b

  Figure4b

A c c e p t e d M a n u s c r i p t range) and for 10 g/l at 20 • C (low concentration range).


Christian Vessaz
email: christian.vessaz@epfl.ch

Loïc Andolfatto
email: loic.andolfatto@epfl.ch

François Avellan
email: francois.avellan@epfl.ch

Christophe Tournier
email: christophe.tournier@lurpa.ens-cachan.fr

Christophe Tournier Lurpa

Pelton Runner Pelton

Toward design optimization of a Pelton turbine runner

Keywords: Pelton turbine, Bucket shape parameterization, Design optimization, High-dimension, Finite volume particle method

The objective of the present paper is to propose a strategy to optimize the performance of a Pelton runner based on a parametric model of the bucket geometry, massive particle based numerical simulations and advanced optimization strategies to reduce the dimension of the design problem. The parametric model of the Pelton bucket is based on four bicubic Bézier patches and the number of free parameters is reduced to 21. The numerical simulations are performed using the finite volume particle method, which benefits from a conservative, consistent, arbitrary Lagrangian Eulerian formulation. The resulting design problem is of High-dimension with Expensive Black-box (HEB) performance function. In order to tackle the HEB problem, a preliminary exploration is performed using 2'000 sampled runners geometry provided by a Halton sequence. A cubic multivariate adaptive regression spline surro-

Introduction

Over the past decades, the production of renewable energy has been constantly growing. This expansion is very likely to accelerate considering many countries reinforced their renewable energy policies. This growth includes the hydro power production at a similar pace as the other renewable energy sources.

In this context, the exploitation of hydro power potential becomes one vector of this expansion towards more renewable. The Pelton turbine is the most popular machine type for the exploitation of high head and low discharge power plants. Since the early water wheel concept featured with several double half-cylindrical buckets patented by Lester Pelton in 1880 [START_REF] Pelton | Water-wheel[END_REF], the geometry of Pelton runners illustrated in Fig. 1 has been considerably improved. Most of the progress made take their roots in practical experience, know-how and extensive experimental tests.

In the new context where harvesting small hydro potentials can become economically viable, there is a need to provide solutions to reduce the design cycle time and the design cost for Pelton runners. Such objectives commonly relies on the use of numerical simulation tools and optimization techniques to solve the runner design problem.

The design problem addressed in this paper can be informally described as finding the runner geometry providing the best performance. The implementation of a design methodology to solve it therefore relies on three underlying main pillars [START_REF] Falcidieno | Processing free form objects within a Product Development Process framework[END_REF]. The runner geometry can be defined by a set of design parameters instantiating a parametric model. The performance of a runners can be evaluated according to a performance indicator, generally computed through numerical simulations. Finally, finding the design parameters leading to the best runner with respect to this performance indicator relies on the implementation of an appropriate optimization strategy.

The work presented in this paper aims at highlighting and circumventing the hurdles toward the optimal design of Pelton runners thanks to numerical simulation. The literature review of Section 2 reveals the obstacles related to he high-dimension of the problem and the computational expensiveness of the numerical simulation together with the consequences on the associated optimization problem.

A geometric model is proposed and presented in Section 3. It is based on 4 bicubic Bézier patches modeling the inner surface of the bucket and on a definition of the outer surface thanks to a thickness map. A major effort is made to keep enough degrees of freedom to avoid suboptimal geometries while reducing the number of free parameters to z P = 21.

The performance of the runners is evaluated according to the simulated torque under specified operating conditions. The FVPM simulation setup providing the simulated torque is detailed in Section 4. It allows a coarse resolution for fast simulations in the exploration phase conducted in the following of the paper as well as a fine resolution well suited for the actual solving of the optimization sub-problems resulting from this exploration.

Finally, Section 5 details the initial exploration of the design space. The emphasis is put on the strategy proposed to reduce the dimension of the optimization problem. It relies on the combination of a design parameter importance ranking together with an original clustering approach. Thanks to this exploration strategy, the initial problem that would not have been achievable with finite computing resources is decomposed into 4 sub problems of lower dimension that can be realistically addressed with conventional optimization techniques. The resolution of these smaller optimization problems is not presented to avoid impeding the paper with technical aspects of low methodological value to the scientific community.

2 Literature review

Geometrical modeling

There are very few works in literature about the 3-D modeling of Pelton turbine buckets. Generally, the models consist in a geometric parameterization based on distances and angles which correspond to hydraulic engineering parameters. However, the development of performing products relies more and more on Computer-Aided Engineering (CAE) and thus on parametric models [START_REF] Sobieszczanski-Sobieski | Multidisciplinary aerospace design optimization: survey of recent developments[END_REF]. In the case of products with complex shape as for the Pelton runners, Bézier, B-spline or Nurbs patches are usually implemented. The critical problem lies in the number of design parameters that must be large enough to provide sufficient shape diversity but small enough to allow an efficient exploration of the design space.

The approach proposed by Anagnostopoulos and Papantonis (2012) consists in defining a 2-D boundary curve of the bucket and a deepest point in order to build intermediate frames and consequently the lateral surface by interpolation. The main design parameters are the length and width dimensions of the 2-D boundary, the coordinates of the deepest point and its depth. The intermediate frames dimensions are defined with scaling factors, each intermediate frames introducing 3 additional variables. Once the lateral surface of the bucket is generated, the cutout is constructed at the intersection of the bucket surface with a toroidal beam modeling the water jet and defined by two design variables representing its axis location and radius. According to the authors, the construction of the entire inner surface of the bucket can be controlled by 19 geometric variables, which is quite small regarding its complexity. The drawback of this method is that the splitter tilt angle cannot be modified as well as the cutout geometry. This approach has been reused by [START_REF] Solemslie | A reference Pelton turbine design[END_REF] to build the NUBRS description of the Pelton bucket for their reference Pelton turbine design. [START_REF] Michálková | Imposing angle boundary conditions on B-spline/NURBS surfaces[END_REF] propose a geometrical description based on uniform bicubic B-splines surfaces composed of 7 x 4 patches to model the inner surface of the bucket. The novelty consists in trying to satisfy prescribed angle distribution along the boundary curve of the bucket which is useful to control the tangential velocities at the bucket outlet. However, as the exact solution exists in very special cases, an approximate solution is proposed. The drawback of this approach is the use of 7 x 4 bicubic patches which are determined by 84 3-D control points which potentially represents 252 parameters. Although reductions are possible, that still represents a significant number of parameters to carry out the geometrical optimization.

Numerical simulation

The deviation of a high-speed water jet by the rotating Pelton buckets is a challenging fluid mechanics problem, which involves complex geometries, moving boundaries, free surface flows and high-pressure variations. The design approach combining numerical simulation and experimental model testing, as described by [START_REF] Mack | Modernization of vertical pelton turbines with the help of cfd and model testing[END_REF], can yet not be applied in the early stages of the design process addressed in this paper. The simulation of a Pelton runner can be investigated using either Eulerian grid-based or Lagrangian particle based numerical methods. [START_REF] Mack | Numerical investigation of the flow in a Pelton turbine[END_REF] and [START_REF] Jošt | Numerical prediction of a Pelton turbine efficiency[END_REF] highlighted the grid influence on the efficiency prediction as well as the needs of significant computing resources to compute accurately a Pelton runner using the twophase homogenous model. [START_REF] Xiao | Numerical simulation of unsteady free surface flow and dynamic performance for a Pelton turbine[END_REF] performed a Volume Of Fluid (VOF) simulation of a rotating Pelton runner. The computed efficiency is a bit lower than the experiments and the dependance on the grid resolution is also highlighted. [START_REF] Marongiu | Free surface flows simulations in Pelton turbines using an hybrid SPH-ALE method[END_REF] demonstrated that particlebased methods are well suited to compute the flow in a Pelton runner. However, particle-based methods require significant computational ressources. [START_REF] Židonis | State of the art in numerical modelling of Pelton turbines[END_REF] present a comparison of different Eulerian and Lagrangian solvers applied to the simulations of rotating Pelton buckets. [START_REF] Anagnostopoulos | A fast lagrangian simulation method for flow analysis and runner design in pelton turbines[END_REF] and [START_REF] Xiao | Numerical prediction of dynamic performance of Pelton turbine[END_REF] proposed a fast Lagrangian computation to design Pelton runners. However, this method is only based on the inlet and outlet velocity vectors of the particles, which provides an estimation of the integrated pressure. Neither the whole pressure field nor the exact water sheet location can be accurately computed.

In 2015, [START_REF] Vessaz | Flow simulation of jet deviation by rotating pelton buckets using finite volume particle method[END_REF] investigated the use of Finite Volume Particle Method (FVPM) to simulate the flow in five rotating Pelton buckets. FVPM is a particle-based solver introduced by [START_REF] Hietel | A finite-volume particle method for compressible flows[END_REF][START_REF] Hietel | A finite-volume particle method for compressible flows[END_REF]. In 2009[START_REF] Nestor | Extension of the finite volume particle method to viscous flow[END_REF] extended the method to incompressible flows. This method features an Arbitrary Lagrangian-Eulerian (ALE) formulation, which means that the computing nodes can either move with the material velocity or a user-prescribed velocity. This method is able to satisfy free surface and no-slip wall boundary conditions precisely. FVPM combines attractive features of Smoothed Particle Hydrodynamics (SPH) and conventional grid-based Finite Volume Method (FVM). It also features the ability to include additional physics, such as silt laden flow erosion [START_REF] Jahanbakhsh | Development of a finite volume particle method for 3-D fluid flow simulations[END_REF]. Therefore, it appears as an appropriate solution for numerical estimations of the performance of Pelton runners in the presented study.

Optimization techniques

According to the statements from subsection 2.1, the design of a Pelton runner is likely to be high-dimensional. The numerical simulation techniques only allows to have an implicit -or black-box -evaluation of the performance of a runner at a noticeable computing cost [START_REF] Vessaz | Flow simulation of jet deviation by rotating pelton buckets using finite volume particle method[END_REF].

Straightforward implementations of usual optimization techniques are generally unsuccessful to solve such High-dimension Expensive Black-box (HEB) problems [START_REF] Shan | Survey of modeling and optimization strategies to solve high-dimensional design problems with computationally-expensive black-box functions[END_REF]. The expensiveness of the problem can -until a certain extent -be addressed by committing computing power accordingly. In their survey, [START_REF] Shan | Survey of modeling and optimization strategies to solve high-dimensional design problems with computationally-expensive black-box functions[END_REF] stated that a key in the HEB context is to implement a strategy to tackle the high-dimensionality.

No specific physical reason came out to help decomposing the Pelton design problem into several independent sub-problems of lower dimensions. Then, the most common approach to handle the high-dimensionality is to work on the design space reduction by two means:

reducing the dimension by removing some design parameters or by transforming a set of correlated variables into a new smaller set of uncorrelated variables with an acceptable loss on the performance prediction; reducing the range of some design variables to a relevant portion of the initial design space.

One popular approach to tackle high-dimensionality is to perform a Principal Components Analysis (PCA) to identify the most important combination of design parameters [START_REF] Shan | Survey of modeling and optimization strategies to solve high-dimensional design problems with computationally-expensive black-box functions[END_REF].

When such properties are not encountered, the importance of the design parameters can still be evaluated through a sensitivity analysis. Usual techniques are also limited because of the high-dimension [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF] and because of the expensiveness of the performance function Caniou (2012). In the case of HEB, one straightforward way to rank the importance of the design parameters consists in building a surrogate model f * of the performance function f Caniou (2012) with a reduced computing cost. The surrogate model is not used to predict the performance at unexplored points but it serves to run the sensitivity analysis without increasing the overall computing budget. This approach proved efficient for the problem addressed in this work. It serves as a basis for the method proposed in Section 5 to define sub design spaces in which usual optimization techniques are more likely to be applied successfully.

Bucket geometrical modeling

General approach

The proposed parametric model of a Pelton bucket is built considering 7 parameters giving the general dimensions such as length, width and depth as represented in Fig. 2. The intersection of the water jet axis with the outer edge plane of the bucket is named O and set as the origin of the bucket coordinate system. As the bucket is assumed to be symmetric with respect to the (X, Y ) plane, only the half bucket with positive coordinates along the Z axis will be described.

Then, a set of physical points with specific properties detailed in subsection 3.2 are used to split the inner surface into four bicubic Bézier patches. These patches are defined in subsection 3.3.

The Pelton bucket performance is known to be enhanced when the outlet angles β 1 decrease. This assertion is limited by the risk of heeling illustrated in Fig. 3: up to a certain angle, the water jet impinges on the outer surface of the next bucket with an associated loss of energy. To account for this heeling phenomenon, the outer surface of the bucket is defined thanks to a thickness map, as explained in subsection 3.4.

Definition of physical points

A set of physical points with specific properties of location or tangent depicted in Fig. 4 are defined on the inner surface:

-T is the tip of the bucket, located at (C 2 , -y t , 0); -Ce is the other extremity of the cutout, with the coordinate (A 2 , 0, E/2); -Cb is the bottom point of the cutout edge, where the edge tangent is normal to the Y axis, with the coordinate (x Cb , -y Cb , z Cb ); -Bt is the bottom of the inner surface where the surface normal is oriented along Y , with the coordinate (x Bt , -F, B/4); -Om is the extreme point of the outlet edge along the Z direction, with the coordinate (x Bt , 0, B/2); -Sm is the intermediate point of the inlet edge, with the coordinate (x Bt , y Sm , 0); -Se is the extreme point of the inlet edge, with the coordinate (-C 3 , y Se , 0); -Oe is the extreme point of the outlet edge along the X direction, with the coordinate (-A 1 , 0, B/4).

The inlet orientation angle α Sm is defined between the inlet edge and the X direction, as depicted in Fig. 5. Then the coordinates of the physical points Sm and Se along the Y direction can be expressed by eqs. ( 1) and (2).

y Sm = C 2 • tan α Sm -y T
(1)

y Se = C • tan α Sm -y T (2)
The definition of the physical point finally requires the four additional parameters on top of the general dimensions previously defined: α Sm , y T , x Cb , y Cb , z Cb . The fifth parameter x Bt is fixed to 0.

Definition of the inner surface

Four bicubic Bézier patches are defined on the physical points. A C 1 continuity between the four surfaces is ensured by imposing a symmetry of the four control points at each vertex shared by several surfaces.

The control points of the first patch are built according to the scheme described in Fig. 6 outlet surface in Oe is oriented with an angle β 1 around Z. The angle α Se is arbitrarily set to 45 • in order to create a smooth edge between Se and Oe. This edge is unlikely to have an influence on the hydrodynamic behavior of the bucket as it usually do not receive water [START_REF] Vessaz | Finite particle flow simulation of free jet deviation by rotating Pelton buckets[END_REF].

For the third patch, 8 out of the 16 control points are derived from the second patch to ensure continuity, as depicted in Fig. 6(c). Only three of its vertex are physical points (Oe, Om and Bt). The fourth vertex named Or and the associated edge orientation are built according to the definition given in Fig. 7. It requires three scalar parameters χ Or , φ Or and ψ Or fixed to 0.75, 0.5 and 0.075 respectively. The outlet surface is oriented with the angle β 1 .

Only the three control points around Ce remain free for the fourth patch. The other inherited control points are pictured in Fig. 6(d). The outlet edge is oriented with an angle α Ce around Y in Ce and the cutout edge in Ce is oriented directly along Y . The position of the other control points are defined by the two distances L y,Ce and L xz,Ce .

Finally, the inner surface requires 8 variables -or free parameters -and 13 fixed ratios of the general dimensions to be completely defined.

Definition of the outer surface

The outer surface is built by offsetting the patches of the inner surface along their local normal vectors with a given thickness. P uv represents a point of one patch of the inner surface with the parameters (u, v). The normal at this point is written n uv . Given a thickness function T defined on [0, 1] 2 that returns a thickness for each (u, v) parameters, the point P uv of the associated patch of the outer surface is given by:

P uv = P uv + T (u, v) • n uv (3)
For each patch of the bucket inner wall, a thickness map is defined by a set of m × n control thickness t ij . The thickness function is of the form defined in (4) with B ij , the Bernstein polynomials, basis functions of the Bézier patches.

T : [0, 1] 2 -→ R u, v -→ m i=0 n j=0 B im (u) • B jn (v) • t ij (4)
Three values t 1 , t 2 and t 3 are used to define the thickness maps associated to the four inner surfaces. In this study, the thickness are kept constant equals to 3 mm, 8 mm and 6 mm respectively. The sets of control thickness are graphically represented in Fig. 8. When generating the entire bucket by symmetry, the sampled points with negative coordinates along Z will be removed.

Closing the bucket volume

The inner surfaces and outer surfaces edges are collections of isoparametric curves. A pair of curve C i and C o with parameter written p defined on adjacent edges are considered to define the joining surface S between them as given in eq. ( 5) to yield half of a bucket.

S : [0, 1] 2 -→ R p, α -→ α • C i (p) + (1 -α) • C o (p) (5)
A symmetry with respect to the (X, Y ) plane is applied to provide the entire bucket depicted in Fig. 4.

Implementation of the model for numerical simulation

The proposed model defines a Pelton bucket as a collection of parametric surfaces S(u, v). Despite the effort made to reduce the number of free parameters, geometrical model of a bucket depend on the 19 design parameters listed in Table 1. Locating the bucket in the runner coordinate system also requires 2 additional parameters Y 0 and X 0 . Given a sampling strategy for the parameters u and v for each surface adapted to the simulation requirement, a discrete representation of the bucket can be generated for numerical simulations. Table 1 Inventory of the parameters describing the Pelton bucket geometry.

Category Parameters

General dimensions A 1 , A 2 , B, C 2 , C 3 , E, F Physical points α Sm , y T , x Cb , y Cb , z Cb Inner surface α T , α Ce , β 1 , β 1,Cb , β 1 , L y,Ce , L xz,Ce Bucket location Y 0 , X 0 4 Numerical simulations

Finite volume particle method

In FVPM, the Sheppard interpolating or shape function ψ is used to discretize the governing equations. The Sheppard function is zero-order consistent and is defined as:

ψ i (x) = W i (x) σ (x) (6)
where

W i (x) = W i (x -x i , h
) is the kernel function and σ (x) = j W j (x) is the kernel summation. The spatial resolution of the interpolation is given by the smoothing length h. In the present study, a rectangular top-hat kernel is used to compute the interaction vectors, which reads:

W i (x) = 1 x -x i ∞ ≤ h, 0 x -x i ∞ > h. (7)
The control volumes are replaced by particles and the exchange occurs through the interfaces defined by overlapping regions. For each pair of overlapping particles, two interaction vectors are defined. Their difference ∆ ij is analogous to the area vector in FVM and is defined as:

∆ ij = Γ ij -Γ ji (8)
which depends on the interaction vector between particles i and j:

Γ ij = Ω ψ i ∇W j σ dV = Ω W i ∇W j σ 2 dV. ( 9
)
Due to the complexity of shape functions, their integrations are usually approximated using quadrature rules over a large number of integration points. In 2011, [START_REF] Quinlan | Fast exact evaluation of particle interaction vectors in the finite volume particle method[END_REF] developed a new FVPM in which the integrals are computed quickly and exactly for 2-D simulations. They simplified the shape functions to circular top-hat kernels and achieved a reasonable compromise between computational cost and accuracy. Recently, Jahanbakhsh et al ( 2014) introduced rectangular top-hat kernels to compute quickly and exactly the integrals in 3-D. A 2-D example of particles interactions with rectangular support is given in Fig. 9(a). The top-hat kernel is less smooth than a bell-shaped kernel as shown by the contours of the Sheppard shape function given in Fig. 9(b). However, [START_REF] Quinlan | Fast exact evaluation of particle interaction vectors in the finite volume particle method[END_REF] demonstrated that top-hat kernel allows a fast and exact computation of the interaction vector in 2-D with a circular support. In 3-D, [START_REF] Jahanbakhsh | Finite volume particle method for 3-D elasto-plastic solid simulation[END_REF] showed that the use of top-hat kernel with a rectangular support reduces significantly the cost of the integral computations in eq. ( 9). Therefore, the latter is simplified as:

Γ ij = - m l ∆S l σ + l σ - l ( 10
)
where m is the number of partitioned rectangles, ∆S represents the surface vector of the partitions, σ -and σ + are the summation kernel inside and outside the surfaces respectively. An outline of the 2-D computation of eq. ( 10) is given in Fig. 9(c), where the rectangular partitions are simplified as lines segments. In this example, 4 segments are required to compute the summation of eq. ( 10) for particles i and j respectively. The water flow is assumed inviscid and weakly compressible. The flow motion is governed by the mass and linear momentum conservation equations:

dρ dt = -ρ∇ • C and d (ρC) dt = ∇ • σ + ρg (11)
where ρ is the density, C is the velocity vector, g is the gravity vector and σ = -pI + s is the stress tensor, which includes p the static pressure and s the deviatoric stress contribution. In the case of an inviscid flow, the deviatoric stress contribution is equal to zero. However, in the present study, an artificial viscosity is introduced to damp the numerical oscillations [START_REF] Vessaz | Finite particle flow simulation of free jet deviation by rotating Pelton buckets[END_REF]. The static pressure is computed from the barotropic equation of state:

p = ρ • a 2 7 ρ ρ • 7 -1 (12)
where ρ • is the reference density and a is the sound speed. According to the weakly compressible assumption, the sound speed is set to 10 • C max , C max being the discharge velocity of the water jet. This assumption is based on the weakly compressible approach of [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF] to ensure that density variations remain below one percent and the Mach number is limited to 0.1 all along the numerical simulation, which is usual for particle-based models derived from SPH formulation. Therefore, the weakly compressible approach allows to increase the time step value compared to a pure compressible simulation by decreasing the sound speed in the CFL condition. The time integration is performed using a second-order explicit Runge-Kuta scheme and the time step is computed by:

∆t = 0.6 • min h a + C i . (13)
In order to stabilize the numerical simulations, a correction term is applied to the mass flux following Jahanbakhsh (2014) and the AUSM + scheme of [START_REF] Liou | A sequel to AUSM: AUSM+[END_REF] is used. Moreover, a particle velocity correction is computed at each time step to ensure a uniform distribution of particles in the flow and avoid particles clustering [START_REF] Vessaz | Finite particle flow simulation of free jet deviation by rotating Pelton buckets[END_REF].

To impose the solid boundary condition, one layer of wall boundary particles is located on the interface. The wall boundary particles have the property of the fluid particles, i.e. their pressure and stress are computed from governing equations of the fluid. However, their velocities are imposed equal to the wall velocity to ensure that the wall boundary particles remain attached to the solid interface. Consequently, the force applied on the boundary is given by:

f B,i = j∈fluid (-p ij I + s ij ) • ∆ ij . (14)

Input and initial setup

In the present study, the values of the following parameters are arbitrarily set in order to have a well defined operating point for the exploration process. First, the water jet parameters are imposed as follow: the orientation of the jet is in the -X direction and its inlet is located at X = 0.185 m, Y = -0.15 m and Z = 0.0 m. The discharge velocity of the water jet is C max = 30.0 m s -1 with a diameter D 2 = 0.03 m. Second, the number of buckets of the Pelton runner is set to 20 and the rotational speed is imposed at 955 rpm.

In order to decrease the computing time, only 2 buckets are used to represent the Pelton runner. The total torque is deduced from the torque evolution in the first bucket. However, this assumption does not capture the heeling phenomenon. In this preliminary study, this phenomenon is avoided by setting a minimal β 1 value, which is large enough, and checking that the healing phenomenon does not occur for the optimized geometries. In further studies, three buckets will be used and the total torque will be deduced from the second bucket. Consequently, only 0.013 seconds are simulated, which corresponds to a rotation angle of 75 • , and is sufficient to compute the torque evolution in the first bucket. The initial setup of the simulation is presented in Fig. 10, which includes:

the geometry of the bucket; two additional parameters X 0 and Y 0 which set the bucket location according to the X and Y Cartesian coordinate respectively; the initial rotation θ 0 which is deduced from the bucket tip in order to obtain an angle of 0.5 • between its location and its first impact through the water jet; and an initially developed portion of the water jet.

The numerical simulations are performed with the FVPM solver SPHEROS developed by [START_REF] Jahanbakhsh | Implementation of a parallel SPH-FPM solver for fluid flows[END_REF]. An example of the SPHEROS results is given in Fig. 11. The particle-based representation uses the instantaneous wall pressure field to render the buckets particles. The visualization of the results is performed using the rendering software ParaView [START_REF] Ayachit | The ParaView Guide: A Parallel Visualization Application[END_REF].

Torque computation

During the simulation, the torque is computed for each bucket and at each time step according to:

T = i∈wall R i × f B,i (15)
where R i is the radius between the runner axis and the particle position X i . The evolution of the torque in each bucket as well as the total torque applying on the two buckets are given in Fig. 12 for the finest particles resolution investigated, i.e. D 2 /X ref = 50. In order to set the objective of the optimization process, the mean torque applied on the runner has to be evaluated from the torque evolution for bucket 1. Therefore, the torque evolution for bucket 1 is resampled according to a given ∆θ = 0.025 • increment. Then, the torque applied on buckets 2 to 20 are deduced by shifting the torque evolution of bucket 1. Finally, the total torque is evaluated by summing the torque evolution of the 20 buckets. An example of the runner torque evolution is given in Fig. 13 for the particles resolution D 2 /X ref = 50. The mean value, as well as the standard deviation, are computed to obtain global variables for the optimization process.

The convergence of the results according to the spatial discretization is shown in Fig. 14. Indeed, the FVPM ensures the convergence of the results thanks to its conservative and consistent formulation [START_REF] Vessaz | Flow simulation of jet deviation by rotating pelton buckets using finite volume particle method[END_REF]. This convergence is also highlighted by the mean runner torque in Table 2.

However, increasing the spatial resolution also increases drastically the computing time required for the simulations, which is highlighted in Table 2 for the five particles resolutions investigated. For the following optimization process, a coarse resolution D 2 /X ref = 10 is selected in order to evaluate many different bucket geometries in a reasonable computing time. Choosing a consistent design space, i.e the explored design space generates buckets shaped geometries, the simulations uncertainty is assumed to be identical for all the simulations because the physics of the flow remain the same.

The drawback of using a coarse resolution is the increase of numerical noise as highlighted by the standard deviation in Table 2. In the present paper we focus on the exploration of the design space in order to reduce the dimension of the problem. However, a finer resolution with reduced noise and uncertainty has to be used in further studies to actually solve the reduced optimization problem.

Design space exploration

Design space definition

The design parameters listed in Table 1 are gathered in a design parameter vector x of dimension z P = 21. The Table 3 shows the lower bound x i,min and the upper bound x i,max of the explored range for the associated

T [N.m] D 2 / X ref = 10 D 2 / X ref = 30 D 2 / X ref = 50
Fig. 14 Influence of the spatial discretization on the phase history of torque for bucket 1. geometrical parameter. The design space Ω is defined by eq. ( 16).

Ω = z P i=1 [x i,min , x i,max ] (16)
The explored range has been arbitrarily defined. This process is made easy thanks to the physical meaning associated with the parameters. Nevertheless, the consistency of this design space has been checked by sampling N test = 500 random input vectors and visualizing the associated runners without noticing any unrealistic one.

The performance T (x) representing the mean torque applied on the runner defined by the parameter vector A 1

[mm] 30.0 60.0 x over one rotation is evaluated through FVPM simulations with the setup described in section 4. The initial exploration of the design space was conducted with a Halton sequence containing N sp = 2000 design points. The main advantage of using a Halton sequence for the exploration lies in the ability to dynamically increase the size of the exploration sample if needed. Every new explored point of the design space will improve the uniformity of the sample [START_REF] Halton | Algorithm 247: Radical-inverse quasirandom point sequence[END_REF]. [START_REF] Iooss | Numerical study of the metamodel validation process[END_REF] reported the efficiency of such low discrepancy sequences for the construction and the validation of surrogate models in high-dimension while pursuing similar objectives of global and uniform covering of the entire design space for exploration.

It required 50 hours of computation distributed on 20 nodes with 2 Ivy Bridge Intel Xeon E5-2650 v2 processors with 8 cores each. The average performance of the investigated population reaches 47.5 N.m. The best point provided 60.3 N.m while the worst leaded 19.6 N.m.

Dimension reduction

The exploration sample showed the high-non linearity of the performance function. No correlated effect of the design parameters was found either. The dimension reduction performed therefore only relies on sensitivity analysis.

A cubic Multivariate Adaptive Regression Spline model (MARS) [START_REF] Friedman | Multivariate adaptive regression splines[END_REF] with 108 non-constant basis functions has been built with the 2000 explored points. The MARS model has been selected for its ability to outperform other surrogate modeling techniques, such as neural networks, polynomial chaos expansion, support vector regression or Kriging, when the input dimension is higher than a dozen [START_REF] Andolfatto | Assistance à l'élaboration de gammes d'assemblage innovantes de structures composites[END_REF]. In the present case, the input dimension is 21.

As for polynomial chaos expansion, the MARS surrogate modeling framework is also well suited to evaluate the importance of each design parameter on the performance function using only the initial exploration sample. On one hand, the model consists in a sum of piece-wise cubic spline basis functions. It is therefore possible to estimate the loss of quality of a model containing all the basis functions but the ones involving one design parameters. Friedman proposed to measure the importance of a design parameter through the difference of the Generalized Cross Validation error (GCV) between the full identified model and a model without basis functions involving this design parameter [START_REF] Friedman | Multivariate adaptive regression splines[END_REF]. The GCV error between the mean torque T and its surrogate model T * is computed according to eq. ( 17), where z M is the number of parameters of the MARS model.

GCV = 1 N sp Nsp j=1 T ( j x) -T * ( j x) 2 1 - z M N sp 2 (17)
Removing basis functions involving one design parameter decreases the number of parameters z M of the model but increases the difference between the actual values T and the modeled values T * . Therefore, a high increase of the GCV when a design parameter is removed implies a high importance of this design parameter.

On the other hand, the analytic expression of the model allows to perform directly an ANalysis Of VAriance (ANOVA). The impact of each design parameter on the performance over the design space is estimated according to the cumulative standard deviation σ of the basis functions involving this parameter. It can be interpreted in a manner similar to a standardized regression coefficient in a linear model [START_REF] Friedman | Multivariate adaptive regression splines[END_REF].

With this strategy, the surrogate model is never used to predict the performance at unexplored points. The prediction error is therefore not a direct matter of concern. Yet, the quality of the surrogate model must still be assessed to ensure that all the trends of the performance function are well captured. It can be measured according to the coefficient of determination R 2 defined in eq. ( 18). The strategy is only applicable if R 2 is close enough to one.

R 2 = 1 - Nsp j=1 T ( j x) -T * ( j x) 2 Nsp j=1 T ( j x) -T ( j x) 2 (18)
with:

T ( j x) = 1 N sp Nsp j=1 T ( j x) (19)
It the present case, the coefficient of determination R 2 reached a value of 0.9773, which has been considered sufficient.

Table 4 presents the increases of generalized cross validation error ∆GCV and the cumulative standard deviation σ related to each design parameters. The parameters are sorted in decreasing importance and the ranks R GCV and R σ of each parameter are also presented in Table 4. The two importance indicators provide almost the same ranking between parameters. The five parameters which have the most influence are related to the runner diameter, cutout shape and bucket depth. The seven last parameters have a negligible influence on the performance function. The associated dimensions of the design space can therefore be left unexplored in further studies. These parameters are related to the physical point Ce, which is not an active surface for the torque generation. The inlet and outlet angles, β 1 and β 1 , feature also a negligible influence due to the well chosen range of exploration, i.e. as small as manufacturable and to avoid the heeling phenomenon.

Range reduction thanks to clustering

The range reduction aims at identifying shrunk area of the design space in which solving an optimization problem will require less efforts. Many methods have been reported in the literature [START_REF] Shan | Survey of modeling and optimization strategies to solve high-dimensional design problems with computationally-expensive black-box functions[END_REF], but most of them require a dynamic sampling of the original design space. In the present case, a method based on the already explored design points is proposed to avoid supplementary computation. It consists in the identification of areas of high performance within the design space. To do so, the proposed algorithm yields clusters of runners with similar geometries and high performance. The design points in each cluster are further analyzed to define a sub design spaces of reduced dimension and reduced range in which an optimization problem can be solved. For the proposed approach, only the N best runners and their associated design points j x j∈1,...,N are considered. A distance matrix D is computed. Each term of the matrix is defined according to eq. ( 20).

d jk = j x r -k x r 2 (20)
where the standardized design point j x r is the image of j

x by the linear application that maps the design space Ω onto [0, 1] 21 . The matrix is symmetric with zeros on the diagonal. An initial graph with N nodes representing the N runners is drawn without initial arc. The set of arcs in the graph is denoted P. Then, the closest runners j and k satisfying eq. ( 21) are sequentially searched and connected in the graph until each runner has at least one arc. To define a new reduced design space associated to each cluster, the range ∆ r,i of the reduced coordinates x r,i associated with each design parameter i is computed. If the range falls below an arbitrary threshold ∆ min , meaning that the associated design parameter is almost constant within the cluster, the associated dimension no longer requires to be explored. Otherwise, the associated dimension must still be explored in the vicinity of the domain occupied by the design points of the cluster.

The design parameters with a negligible impact on the performance function identified in the previous subsection can be fixed at the mean value within the cluster.

The Table 5 shows the four sub design spaces built with ∆ min fixed at 5%. Four design problems of lower dimensions -from 7 to 11 instead of the initial 21 -and with smaller range can be formulated based upon these results.

Conclusion

The novel contribution of this paper lies in proposing a framework addressing the major difficulties toward the design optimization of a Pelton runner. It is chained along three main links.

First, the proposed parametric model of a Pelton bucket consists of four bicubic Bézier patches with C 1 continuity. It requires only 21 parameters while keeping enough degrees of freedom to conduct the geometrical optimization. It also features straightforward discretization capabilities to be automatically linked to the numerical simulation solver.

Then, the FVPM simulation is a state-of-the-art convenient and accurate tool to capture the deviation a water jets by rotating Pelton buckets. The total runner torque is evaluated from the torque evolution in a single bucket to reduce the computing time. However, a coarse resolution has to be selected for the numerical simulations in order to evaluate many different bucket geometries in a reasonable computing time for the exploration and dimension reduction purpose pursued in the last part of the paper. But the same simulation setup with a fine resolution can be applied for solving the reduced optimization sub-problems as in [START_REF] Vessaz | Design optimization of a 2D blade by means of milling tool path[END_REF].

Addressing the geometrical modeling and the evaluation of the runner's performance through numerical simulation leads to a High-dimension Expensive Blackbox problem. Such problems can not be addressed directly with usual optimization technique. Therefore, the last aspect of the presented work focuses on reducing the dimension of the problem and the range of the ex- plored parameters. The initial design space is explored in order to identify the design parameters with the highest impact on the performance and the areas of the design space with the highest performance. The original clustering approach presented yields a set of optimization problems of lower dimensions with design spaces covering lower ranges that can therefore be addressed with conventional optimization techniques. This clus- tering approach provides the last link of the proposed chain toward the design optimization of a Pelton runner. For the presented test case, four sub design spaces of dimension 7 to 11 are obtained, each of them leading to an optimization problem noticeably easier to solve than the initial one.

The next steps consists in actually solving these resulting optimization problems. As the performance of the runners in the sub design spaces are likely to vary less than in the initial space, it may be necessary to refine the particles resolution in order to increase the reliability of the simulated torque, even if it will induce an increase of the computing time. The proposed geometrical model based on an analytic description of the runner also allows the easy introduction of constraints related to manufacturing. Furthermore, the FVPM simulation framework is suitable for integrating other physical phenomenon such as mechanical fatigue due to cyclic loads on the buckets or erosion due to silt laden flows.

Fig. 2

 2 Fig. 2 General dimensions of the bucket.

Fig. 4

 4 Fig. 3 Illustration of the heeling phenomenon.

Fig. 5

 5 Fig. 5 Definition of the inlet orientation angle α Sm .

Fig. 6 Fig. 7

 67 Fig. 6 Control points defining the four bicubic Bézier patches.

Fig. 8

 8 Fig.8Variable thickness map defining the outer surface as an offset of the inner surface.

Fig. 9

 9 Fig.9Rectangular support kernels and overlapping regions (a) plotted with the contour of Sheppard shape function for the particle i and top-hat kernels (b), outline of the intersection volume between particles i and j (c)

Fig. 10

 10 Fig.10Outline of the initial setup for the numerical simulations.

Fig. 11

 11 Fig. 11 FVPM simulation of two rotating buckets: particlebased representation (up) and free surface reconstruction of the water sheet (down).

Fig. 12

 12 Fig.12Phase history of the torque for each bucket (colors) and total torque (black) for a spatial discretization of D 2 /X ref = 50.

  Fig.13Phase history of the runner torque evaluated from bucket 1.

  clusters of connected runners with similar geometries.With N = 10, the average distance between point is 1.612. The clustering graph with the smallest distance between runners is presented in Fig.15. It yields four clusters. The associated buckets are depicted in figures 16 to 19.

Fig. 15

 15 Fig. 15 Clustering graph representing the N = 10 best runners identified by their index in the explored sample.

Fig. 16

 16 Fig.16Buckets from the runners of the cluster 1.

  Fig.17Buckets from the runners of the cluster 2.

Fig. 18

 18 Fig.18Buckets from the runners of the cluster 3.

Fig. 19

 19 Fig.19Buckets from the runners of the cluster 4.

Table 2

 2 Influence of the spatial discretization on the mean runner torque, standard deviation and computing time.

	D 2 /X ref	mean(T ) std(T ) Computing time
	[-]	[N]	[N]	[hour]
	10	60.31	8.73	0.5
	20	64.38	3.44	4.5
	30	66.64	2.32	21.2
	40	68.32	2.37	61.2
	50	69.62	2.29	164.7

Table 3

 3 Definition of the design space Ω.

	Parameters	Unit	x i,min	x i,max

Table 4

 4 Relative importance of the design parameters estimated on the MARS model of the performance function, measured thanks to the Generalized Cross Validation error GCV and the cumulative standard deviation σ.

	Parameter ∆GCV	R GCV	σ	Rσ
	X 0	332.06	1	21.87	1
	C 2	24.26	2	5.66	3
	x Cb	22.77	3	6.14	2
	y Cb	20.41	4	3.48	8
	F	18.20	5	4.23	4
	z Cb	16.86	6	4.14	5
	y T	16.33	7	3.97	7
	E	14.63	8	4.01	6
	A 2	10.76	9	2.86	9
	B	10.44	10	1.87	11
	α Sm	8.74	11	1.91	10
	Y 0	2.80	12	0.55	12
	C 3	1.51	13	0.29	13
	β 1,Cb	1.37	14	0.19	14
	A 1	0.00	15	0.00	15
	α Ce	0.00	15	0.00	15
	α T	0.00	15	0.00	15
	β 1	0.00	15	0.00	15
	β 1	0.00	15	0.00	15
	L y,Ce	0.00	15	0.00	15
	L xz,Ce	0.00	15	0.00	15

Table 5

 5 Sub design spaces associated to each cluster defined by fixed value or range of each design parameter.

		Cluster 1	Cluster 2	Cluster 3	Cluster 4
	Parameter	min	max	min	max	min	max	min	max
	X 0	132.1	137.1	146.1		142.6	145.3	138.4
	C 2	52.4	59.6	36.9	40.5	35.9	49.3	60.0
	x Cb	39.8	30.2		39.9	46.5	31.3	33.3
	y Cb	-16.8	-14.3	-15.0		-11.5	-12.4
	F	27.1	31.0	29.0	30.2	16.1	18.5	29.0	32.5
	z Cb	13.2	20.6	13.0		13.4	18.2	22.6
	y T	39.8	-9.4	-8.6	-9.1	-6.9	-8.6	-6.3
	E	37.0	41.5	40.1	41.6	47.0	42.5	47.1
	A 2	43.9	52.8	31.9	54.1	32.8	42.6	34.9	48.2
	B	77.4	93.4	82.9	90.9	83.1	94.7	82.4	106.4
	α Sm	39.8	3.0		4.0	5.4	0.0	1.6
	Y 0	11.8	19.4	20.1		6.9	7.4
	C 3	21.3	41.5	21.0	23.7	38.0	41.5	35.3
	β 1,Cb	50.4	54.0	63.5		55.8	57.9	59.0
	α T	26.6	31.1		28.7	29.6
	β 1	7.6	11.1		9.5	7.6
	β 1	4.5	7.4		3.4	5.3
	α Ce	11.2	-1.4		12.9	4.4
	L y,Ce	6.1	5.1		3.9	2.7
	L xz,Ce	15.2	14.4		15.3	12.9
	A 1	51.2	46.9		42.6	43.1
	Reduced dimension	11		7		10		8