D Djafer
email: djdjafer@gmail.com

A Irbah

M Zaiani
email: zaianimo@gmail.com

Identification of clear days from solar irradiance observations using a new method based on the wavelet transform

Keywords: solar radiation, turbidity parameters, clearness index, wavelet transform

published or not. The documents may come

Identification of clear days from solar irradiance observations using a new method based on the wavelet transform

Djelloul Djafer, Abdanour Irbah, Mohamed

Introduction

Solar radiation passing through the Earth atmosphere keeps useful information on the medium. The solar flux is attenuated by scattering, air mass, clouds, aerosols and various reflections. Physical models that explain the solar flux and its temporal variations measured on the ground are of great importance to know the atmosphere state when observing, especially to know the components that are responsible of flux fluctuations. The analysis of these fluctuations is a passive way to probe the atmosphere in contrast to the use of specific instruments such as lidars and radars. Observations during clear days allow detecting certain atmospheric constituents such as aerosols and track their evolution, which may impact Earth's radiative budget and the climate. Clear day observations are also required for models to derive and predict solar irradiance especially where it cannot be measured. These models are needed to obtain the correct design and output of solar power plants in case of clear sky conditions. This accrued interest comes from the huge investments in solar applications in most of countries where the radiative solar resources need to be accurately characterized either measured or modeled. Validating candidate models for this kind of task is an essential step in the process to provide analysts with sufficient background information on their detailed performance. Several works studying and comparing broad band irradiance models are reported in the literature [START_REF] Gueymard | Clear-sky irradiance predictions for solar resource mapping and large-scale applications: Improved validation methodology and detailed performance analysis of 18 broadband radiative models[END_REF]. These models with high-performance and good accuracy are, as previously mentioned, of great interest for clear sky conditions. To assess the performance and improve the validation of such models it is necessary to have clear sky situations i.e. data free of noise and cloudy periods, to obtain valid performance results [START_REF] Gueymard | Clear-sky irradiance predictions for solar resource mapping and large-scale applications: Improved validation methodology and detailed performance analysis of 18 broadband radiative models[END_REF]. The commonly used tool to characterize sky conditions (classification of day types) over a particular site is the clearness index k t [START_REF] Alves | Effects of Sky Conditions Measured by the Clearness Index on the Estimation of Solar Radiation Using a Digital Elevation Model[END_REF][START_REF] Babatunde | Characteristics Variation of Total Solar Radiation at llorin, Nigeria[END_REF][START_REF] Ideriah | Sky conditions at Ibadan during 1975-1980[END_REF][START_REF] Kumar | Estimation of global radiation using clearness index model for sizing photovoltaic system[END_REF][START_REF] Khem | Estimation of Global Solar Radiation Using Clearness Index and Cloud Transmittance Factor at Trans-Himalayan Region in Nepal[END_REF][START_REF] Kuye | Analysis of solar radiation data for Port Harcourt, Nigeria[END_REF][START_REF] Liu | The interrelationship and characteristic distribution of direct, diffuse and total solar radiation[END_REF][START_REF] Mellit | Methodology for predicting sequences of mean monthly clearness index and daily solar radiation data in remote areas: Application for sizing a stand-alone PV system[END_REF][START_REF] Mellit | Methodology for predicting sequences of mean monthly clearness index and daily solar radiation data in remote areas: Application for sizing a stand-alone PV system[END_REF][START_REF] Okogbue | Review Hourly and daily clearness index and diffuse fraction at a tropical station, Ile-Ife, Nigeria[END_REF][START_REF] Santos | Methodology for generating daily clearness index index values K t starting from the monthly average daily value . Determining the daily sequence using stochastic models[END_REF]. On a clear day, the atmosphere causes a reduction of the extraterrestrial solar input by about 30 per cent and to nearly 90 percent in a very turbid (cloudy) day. k t ranging values depend from one author to the other. Alves et al. [START_REF] Alves | Effects of Sky Conditions Measured by the Clearness Index on the Estimation of Solar Radiation Using a Digital Elevation Model[END_REF] reported that a cloudy sky corresponds to 0 ≺ k t ≺ 0.3, a partially cloudy sky to 0.3 ≤ k t ≤ 0.65 and a clear sky to 0.65 ≺ k t ≺ 1. The last definition was also adopted by Gueymard [START_REF] Kuye | Analysis of solar radiation data for Port Harcourt, Nigeria[END_REF]. Following Bendt et al. [START_REF] Bendt | The frequency distribution of daily insolation values[END_REF], clear sky is when 0.5 ≤ k t ≤ 0.85, greater than 0.5 for Ahmed et al. [START_REF] Ahmed | Estimation of Global and Diffuse Solar Radiation for Hyderabad, Sindh, Pakistan[END_REF] and greater or equal than 0.7 for Molineaux et al. [START_REF] Molineaux | Direct luminous efficacy and atmospheric turbidity -improving model performance[END_REF], Li and Lam [START_REF] Li | An analysis of climate parameters and sky condition classifications[END_REF], Li et al. [START_REF] Li | Overcat sky conditions and luminance distribution in Hong Kong[END_REF] and Eftimie [START_REF] Eftimie | Beam Horizontal Irradiance Simulation for Bra?ov Urban Area -Clear Sky Model[END_REF]. For Iqbal [START_REF] Iqbal | An Introduction to Solar Radiation[END_REF], a clear sky is defined when k t ranges between 0.7 and 0.9 while Reindl et al. [START_REF] Reindl | Diffuse fraction correlation[END_REF] proposed k t > 0.6. The clear sky of most tropical regions for Ndilemeni et al. [START_REF] Ndilemeni | Evaluation of clearness index of Sokoto Using Estimated Global Solar Radiation[END_REF] corresponds for 0.68 ≤ k t ≤ 0.75. The choice of k t interval values is also different from one site to another according to Mellit et al. [START_REF] Mellit | Methodology for predicting sequences of mean monthly clearness index and daily solar radiation data in remote areas: Application for sizing a stand-alone PV system[END_REF]. In addition it varies also from day to day and from month to month as stated by Serban [START_REF] Serban | Estimating Clear Sky Solar Global Radiation Using Clearness Index, for Brasov Urban Area[END_REF] and Ahmed et al. [START_REF] Ahmed | Estimation of Global and Diffuse Solar Radiation for Hyderabad, Sindh, Pakistan[END_REF]. We see clearly that there is no unique limit value of k t neither a direct method for its determination. Therefore, its choice may be controversial to discriminate between clear and turbid days. A bad value will affect the number of clear and turbid days by not selecting the good days. Consequently, the performance of modeling and analysis of solar irradiance data may highly be k t dependent. For this raison, the problem of k t choice is first considered in the present work proposing a method that easily separate clear and cloudy periods. This uses the wavelet transform properties that automatically computes a threshold value to discriminate between clear and turbid days. It is chosen according to the noise level and cloud signature present in the recorded data of solar irradiance. First, we present the algorithm based on the k t mean value of the clearness index to distinguish between a clear and a non clear day. Then, we introduce the wavelet transform method that we compare with the k t mean one considering the number of clear days obtained from five years of data recorded at Ghardaïa city in Algeria [START_REF] Djafer | Estimation of atmospheric turbidity over Ghardaa city[END_REF]. Finally, we base our comparison on the errors induced onto the GSR (Global Solar Radiation) model and in the estimation of turbidity parameters.

The clearness index method to select clear and turbid days

The clearness index is usually used to determine the clear days to perform model comparison and other studies. It was introduced by Liu and Jordan [START_REF] Liu | The interrelationship and characteristic distribution of direct, diffuse and total solar radiation[END_REF] as a parameter that accounts for stochastic property conditions for a given site [START_REF] Peled | Evaluation of solar radiation properties by statistical tools and wavelet analysis[END_REF]. Interval values for the clearness index are used to separate clear and turbid days (see Section 1). They change however from one site to another leading to misinterpretation of the results especially when authors are comparing and studying the clear sky empirical models. We develop then an algorithm based on the instantaneous clearness index to determine automatically complete clear days from a huge data set. The different steps of the algorithm are as follows: 2

1. Selection of global solar radiation records of a given day where the Sun elevation is greater than 10 • . This restriction is due to the thin haze in the early morning or late afternoon that may be present. This may lead to consider clear days as not, if any. 2. Calculation of the instantaneous clearness index k t . It is defined over time t between sunrise and sunset as the ratio between the terrestrial global solar radiation G on a horizontal surface and the extraterrestrial one G 0 . Its mathematical expression is:

k t = G G 0 ( 1
)
G 0 in W/m 2 is given by [START_REF] Duffie | Solar engineering of thermel processes[END_REF]:

G 0 = I sc * 1 + 0.033 cos 360N 365 × [cos(φ) cos(δ) cos(ω) + sin(φ) sin(δ)] (2)
where I sc = 1367 W/m 2 is the solar constant and N the day number in the year (N=1 is the first and N=365 the last day in the year). φ, δ, ω are respectively the latitude of the location, the solar declination angle and the hour angle at Sun rise in degrees.

The data used in this work were collected at the Unit of Applied Research in Renewable Energy (Ghardaïa, Algeria) between (2005)(2006)(2007)(2008). The instruments and methods of the data collection are described in detail in Djafer and Irbah [START_REF] Djafer | Estimation of atmospheric turbidity over Ghardaa city[END_REF]. The key to know is that the three components of solar radiation are recorded every 5 minutes since 2004 together with the temperature and humidity. The instruments that measure the direct, the global and the diffuse solar radiation components are from EKO instruments (http://eko-eu.com/). They are usually cleaned two or three times per week according to the weather conditions. They are also calibrated each three years at the meteorological station of Tamanrasset (latitude = +22.78 o , longitude = +5.51 o , altitude = 1270 m) where the solar irradiance is measured since 1970. The histograms of the instantaneous clearness values obtained with these data are shown in Figure 1. Each histogram in the figure is computed over one year for the period 2005 to 2008. The values of the statistical analysis are given in Table 1. We note that the annual mean and median values of k t are 0.67 and 0.73. We tried several values of k t between 0.68 and 0.74 to select the clear days from the data set and we finally opted for k t = 0.70 as the optimal limit to define them. We should point out that we may have few days classified as clear with this choose when the clearness index are greater than 0.70 (see upper curves in Figure 2) and other omitted if it is lower (see bottom curves in Figure 2). This is a source of error for the models that use the results coming from this kind of data analysis because of the inefficiency to select the boundary parameter of the clearness index.

Discrimination between clear and turbid days using the wavelet method

The wavelet transform decomposes signals simultaneously in the time and frequency domains. This allows to better control their main properties and to process the noise. The discrete version of the wavelet transform is implemented with several algorithms among them the à trous algorithm [START_REF] Holschneider | On the wavelet transformation of fractal objects[END_REF][START_REF] Shensa | The discrete wavelet transform: wedding the à trous and Mallat algorithms[END_REF][START_REF] Starck | Image Processing and Data Analysis: The Multiscale Approach[END_REF]. We will use the wavelet transform implemented with the à trous algorithm to automatically detect the signatures of clouds and the noise magnitude in the daily solar irradiance. We chose this algorithm since the transform is known at each decomposition level, its evolution is followed from one scale to the other and it is easy to implement [START_REF] Djafer | Processing Method effect on Sun Diameter Measurement with Solar Astrolabe[END_REF][START_REF] Djafer | Processing Method Effects on Solar Diameter Measurements: Use of Data Gathered by the Solar Disk Sextant[END_REF]. The decision if we have or not a clear day is taken from the decomposition analysis considering some We first start by recalling the principle of the wavelet transform decomposition using the Mexican hat wavelet and the à trous algorithm (see Fig. 3). The first level decomposition of the measured global solar irradiance component (Fig. 3a) generates two signals, which are the approximation and the detail. The first one (Fig. 3b) is a smoothed version of the main signal (low frequencies) while the second one (Fig. 3c) contains the high frequency components. The noise and signature of clouds of the measured global solar irradiance are present in the high frequency components, if any. The same decomposition procedure is repeated on the first level approximation of the signal giving rise to the second approximation and detail (see Fig. 3d and Fig. 3e). It is at this second level decomposition that the decision is taken if the solar irradiance data of interest correspond to a clear or turbid day. We look if the signal detail values are smaller or larger than a determined threshold (Fig. 3e). This threshold is estimated from the multiresolution decomposition support according to Starck et al. [START_REF] Starck | Image Processing and Data Analysis: The Multiscale Approach[END_REF] principle. The multiresolution support M of a signal describes in a logical way whether or not a signal I contains information at the time x for a given scale j of the wavelet decomposition. If M I ( j, x) = 1 or is true then I contains information at the time x for the scale j. The wavelet transform of a signal with the à trous algorithm produces then a set w j for each scale j. The I 0 original signal can be expressed as the sum of its I p smoothed form and all wavelet coefficients of all p decomposition planes:

I 0 = I p + p j=1 w j (3)
Any value of the signal at x time can be also expressed in the same way as follow:

I 0 (x) = I p (x) + p j=1 w j (x) (4)
The multiresolution support M will be obtained by detecting the significant coefficients at each scale. It is defined by: The algorithm to create the multiresolution support is as follows:

M( j, x) = 1 if w j (x) is significant M( j, x) = 0 if w j (x) otherwise 5
1. We compute the wavelet transform of the signal I. 2. We estimate the noise standard deviation at each scale. 3. The binarization at each scale leads to the multiresolution support.

The estimation of the threshold value which is the noise standard deviation σ I is then computed by the following iterative steps:

1. We estimate the standard deviation σ 0 I of the noise of the signal I 0 2. We compute p scales of the wavelet transform of I 0 with the à trous algorithm and then the equation:

I 0 (x) = c p (x) + p j=1 w j (x) (5)
where w j are the wavelet coefficients in each scale and c p is the low frequency part of I 0 . The noise in c p is negligible. 3. We set the iterative step n to 0. 4. We compute the multiresolution support M which is derived from the wavelet coefficients and from σ n I . 5. We select the x time values which belong to the set S : x ∈ S if M( j, x) = 0 for all j in 1...p. S is the set of x time values, which are only attributed to the noise.

6. For all the selected x time values, we compute the values I(x) -c p (x) and the standard deviation σ n+1 I of these values. 7. We increase the iterative step: n = n + 1. This method converges in few iterations giving an estimation of the noise standard deviation, our threshold value. This value obtained for a collected GSR is plotted in Figure 3e together with the second decomposition detail. We have then a clear day if the detail coefficients are less than the constant threshold and it is considered as a turbid day otherwise. The upper curves in Figure 2 show the case of turbid days while clear ones are plotted in the bottom curves in Figure 2.

if

Comparison of the results obtained from the 2 methods

We apply the two methods presented in Section 2 and 3 to process the four years (2005-2008) of global solar irradiance data recorded at Ghardaïa. The monthly average number of clear days we find, is plotted in Figure 5 (left curves). We observe that the number of clear days obtained with the wavelet method is in average lower than what we have using the clearness index criteria. 490 clear days are obtained with the wavelet method while we find 623 with the clearness criteria method. A difference of 143 days exists when we process the same whole data set with the 2 methods. This difference is significant in statistical studies and of course depends on the limit value of the modified clearness index. The widely used value 0.70 as seen in see Section 2 was taken in our analysis. Since the clearness index method is based on the average of hourly values of the solar irradiance, this makes it very sensitive to the signature of clouds and aerosols. A turbid day may easily classified as clear due to this mean calculation if the signature of clouds and aerosols present in the solar irradiance component is smaller than few tenths of W/m 2 . This is confirmed by Lopez and Battles [START_REF] Lópze | Estimate of the atmopsheric turbidity from three broad-band solar radiation algorithms. A comparative study[END_REF] who wrote that some data may erroneously be included as a totaly cloud-free atmosphere if a simple radiometric criterion such as k t > 0.7 is used. It is then expected that the clearness criteria method is not sensitive for instance to thin cloud covers and cirrus clouds, which have small effects on the measured global irradiance. Our method has the advantage to detect these fluctuations in the global solar irradiance component. The signatures of either cirrus or thin clouds such as shown in Figure 2 will be transformed in distinguishable frequencies of different magnitudes in the detail coefficients as illustrated in the wavelet transform decomposition of Figure 3. Indeed, the day is classified as clear or turbid according to a constant value calculated automatically following the importance of the perturbations. Nevertheless, we observe in Figure 5 that the curves obtained with the two methods have very similar shapes. We note that the number of clear days is lower during months between July and October. The reason is that the sky in the Ghardaïa region is affected by clouds and sandy winds during these months [START_REF] Djafer | Estimation of atmospheric turbidity over Ghardaa city[END_REF]. The number of clear days per year obtained with the two methods is given in Table 2 and plotted in Figure 5 (right curves). The corresponding percentage of clear days relatively to the whole data set is also in Table 2. The first row gives the number and percentage of clear days obtained with the clearness index criteria and the second one using the wavelet method. As before, we observe from the right curves of Figure 5 that the two curves have similar shapes. We note a maximum of clear days in 2007 for the two methods. This confirms previous results of the Linke turbidity factor obtained from the same data set [START_REF] Djafer | Estimation of atmospheric turbidity over Ghardaa city[END_REF]. We found that it had a lower value during this year. The annual number of clear days is however relatively constant for the other years when the wavelet method is used while it is not the case for the clearness criteria method. We finally compute the GSR using the Iqbal C model, that will be described in the following section, to see how it is impacted by the clear day number estimated by the two methods. We fit with this model each recorded GSR data that corresponds to a clear day. The analysis consists then in the calculation of several fit quality parameters such as the mean absolute error (maee), the root mean squared error (rmsee), the mean absolute relative error (maree), the root mean squared relative error (rmsree), the mean squared relative error (msree), the mean absolute percentage error (mapee), the coefficient of determination (r2) and the correlations coefficient (cf). The values of these parameters are given in Table 3. The wavelet method clearly provides best values for each quality parameter than the clearness index method. We conclude at this stage that our method is more efficient and is not dependent of the site and observations. We will now analyze in the next Section how the estimation of the turbidity parameters is also affected by the number of clear days used to build the GSR model.

Estimation of the turbidity parameters

We present first the innovative approach we use to estimate the turbidity parameters with a GSR model. These parameters are the Angstrom coefficient β, the average surface albedo ρ g , the wavelength Angstrom exponent α, the forward scatterance F c and the aerosol single scattering albedo w 0 . According to Gueymard [START_REF] Gueymard | Clear-sky irradiance predictions for solar resource mapping and large-scale applications: Improved validation methodology and detailed performance analysis of 18 broadband radiative models[END_REF], the Iqbal C model offers a better accuracy than the other more conventional models and we will consider it in the following. We first recall hereafter the main equations of this model described in detail in Iqbal [START_REF] Iqbal | An Introduction to Solar Radiation[END_REF]. The direct normal irradiance I n (W/m 2 ) is given by:

I n = 0.9751I sc E 0 τ 0 τ g τ w τ r τ a (α, β) (6)
where τ 0 , τ g , τ w , τ r and τ a (α, β) are respectively the ozone, gas, water, Rayleigh and aerosol scattering transmittances. I sc and E 0 are respectively the solar constant and the eccentricity correction factor. The aerosol scattering transmittance, which depends of the relative optical air mass m a , the Angstrom coefficient β and wavelength Angstrom exponent α, is given by [START_REF] Canada | Determination of Angström's turbidity coefficient at Valencia[END_REF]:

τ a (α, β) = (0.12445α -0.0162) + (1.003 -0.125α) -βm a (1.089α+0.5123) (7)
We are interesting by expressing the global irradiance I t on a horizontal surface, which is the parameter we measure with our instruments. It is the contribution of of 2 irradiance terms:

I t = I nh + I d ( 8
)
where I nh is the normal irradiance on an horizontal plane and I d is the horizontal diffuse irradiance. The normal irradiance I nh (W/m 2 ) is given by:

I nh = I n sin(h) (9)
where h is the elevation angle of the Sun in degrees.

The horizontal diffuse irradiance I d (W/m 2 ) is a combination of the three individual components, which are the Rayleigh I dr (W/m 2 ) and aerosols I da (W/m 2 ) scattering after the first pass through the atmosphere, and the multiple reflection processes I dm (W/m 2 ) between the ground and sky:

I d = I dr + I da + I dm ( 10
)
I dr which depends on aerosol single scattering albedo w 0 , is given by:

I dr (w 0 ) = 0.79 2 τ 0 τ g τ w (1 -τ r ) I sc sin (h) (1 -m a + m 1.02 a )
τ aa [START_REF] Eftimie | Beam Horizontal Irradiance Simulation for Bra?ov Urban Area -Clear Sky Model[END_REF] where

τ aa = (1 -(1 -w 0 )(1 -m a + m 1.06 a )(1 -τ a )
is the direct radiation transmittance due to aerosol absorptance. I da is related to the forward scatterance F c :

I da (F c ) = 0.79τ 0 τ g τ w τ a F c (1 -τ aa ) τ aa I sc sin (h) (1 -m a + m 1.02 a ) (12)
I dm related to the ground albedo ρ g , is given by:

I dm (ρ g ) = (I nh + I dr + I da ) ρ g ρ a (1 -ρ g ρ a ) (13)
where ρ a is the albedo of the cloudless sky, which can be computed with:

ρ a = 0.0685 + (1 -F c )(1 - τ a τ aa ) (14)
The (1 -F c ) term corresponds to the back-scatterance. The second term on the right hand side of Equation 14 represents the albedo of cloudless skies due to the presence of aerosols, whereas the first term is the albedo of clean air.

The global irradiance I t on a horizontal surface is given by:

I t = I nh + I d = (I n sin (h) + I dr + I da ) 1 1 -ρ g ρ a ( 15
)
We fit the global irradiance component recorded during clear days we extract from the global data set with the Iqbal C model given by Equation 15. The method consists to solve a nonlinear fitting problem in the least-squares sense i.e. we look for the x-vector coefficients (β, ρ g , α, w 0 , F c ) that minimize the following residual function:

I qbal (x) -I tr 2 = i I qbal (x i ) -I tr 2 ( 16
)
where I tr is the recorded global solar irradiance component and I qbal (x) = I t (β, ρ g , α, w 0 , F c ) the Iqbal C model. The result will be a vector x 0 = (β 0 , ρ g 0 , α 0 , w 00 , F c0 ) that best fits the considered solar irradiance component. The left curve of Figure 6 plots for instance a recorded global solar irradiance component of a clear day superposed to its fit by the Iqbal C model. We will now apply this process to our whole solar data set to estimate the temporal variations of the turbidity parameters given by the x-vector components.

The turbidity parameter we consider first is the x-vector component β, which is the Angstrom coefficient. Its monthly average values during the four years is shown by the right curve in Figure 6. We note that this parameter has maximum values during the summer months and minimum ones in winter. This result confirms what was found and stated in our previous paper [START_REF] Djafer | Estimation of atmospheric turbidity over Ghardaa city[END_REF]. It is explained by a hot summer climate and winds of the south sectors (Sirocco) that characterize the region of Ghardaïa. This kind of winds brings particles of dust and sand with them, which directly impact the Angstrom coefficient. We also note that the wavelet method gives higher values compared to the clearness criteria one, especially between May and August. The daily difference between the values obtained by the two methods reaches 0.291 whereas it is 0.010 for the monthly difference or about 15%. The second x-vector component we consider is the surface albedo ρ g . Its monthly average values computed over the four years are plotted in Figure 7 (left curves). The curves obtained from both methods have similar shapes with locally some differences reaching 0.02 or 5%. The daily difference between the values obtained by the two methods can reaches 0.5. We observe that the surface albedo has higher values in spring.

The third component of the x-vector of the model is the Angstrom exponent α. The right curves in Figure 7 show the monthly average values of the Angstrom exponent computed over the four years of solar data. The two methods provide similar results. The mean value of alpha is close to that obtained in 2005 with MODIS data [START_REF] Zaiani | Atmospheric Turbidity Study Using Ground and Orbit Data[END_REF]. A small difference of about 5% is however observed in August and September. The α parameter has high values around July indicating that a relatively high ratio exists between small and large particles of the aerosols [START_REF] Lópze | Estimate of the atmopsheric turbidity from three broad-band solar radiation algorithms. A comparative study[END_REF] in opposite to August. The daily difference between the values obtained with the two methods reaches 2.279 whereas the monthly one is 0.085.

The monthly average values of the fourth x-vector component i.e. the forward scattering parameter F c , is plotted in Figure 8 (left curves). The curves representing F c versus time obtained from both methods are very close. We note that the shape of this parameter behaves in an opposite way compared to the Angstrom coefficient (β) one. This means that aerosol particles are present when β is high and consequently there is weak forward scattering. The daily difference between the values obtained with both methods reaches 0.461 whereas the monthly one is 0.045. Finally, The right curves in Figure 8 plot the monthly average values of the aerosol single scattering albedo parameter w 0 , which is the fifth x-vector component. We observe also that both curves are very close. w 0 has the same shape as the Angstrom coefficient β in opposite with the forward scattering parameter F c . We found that the daily difference between the values obtained by the two methods reaches 0.6 whereas the monthly difference is 0.015.

Conclusion

A new method is presented in the present work to extract clear observation days from solar irradiance measurements and consequently estimate the turbidity parameters. It is based on the wavelet transform properties and allows the determination of all clear days from a huge data set which are useful for modeling studies. This wavelet method is compared to the well-known one called the clearness index method. The number of clear days found by the two methods differs with 143 days when we process four years of data (2005)(2006)(2007)(2008) of global solar irradiance recorded at the Unit Research of Renewable Energies located in Ghardaïa in the south of Algeria. This difference is highly dependent on the limit value used with the clearness index method. We can point out that we may classify few days as clear with this limit value and other ones as not due to the fixed threshold of k t that depends on the location where measurements are performed. Our method avoid this problem and allows to localize the signatures of noise, aerosols and clouds in time and frequency and then distinguish between clear and turbid days independently of the observation location. We used in a second step the Iqbal C model to quantify the effects of the difference of clear days obtained with the two methods. We studied first how the model fits our measurements and secondly we presented a new approach to estimate the turbidity parameters, which are the Angstrom coefficients (β), the wavelength exponent (α), the aerosol single scattering albedo (w 0 ), the forward scatterance (F c ) and the average surface albedo (ρ g ). The statistical analysis revealed that there are some differences between the results to model the GSR and in the estimation of the turbidity parameters. It should be noted that the difference between the two methods could be more significant if another threshold rather than 0.70 was used for the clearness index method. We can conclude that our new method is more efficient and is site independent contrarily to the clearness index method. The other new approach we proposed in this work to estimate the turbidity parameters is also promising and will be used in a near future to study the temporal variations of the turbidity parameters at Ghardaïa. Finally, we notify the persons interested by our method that the algorithms of our method are available.

FrequencyFigure 1 : 2 Figure 2 :

 122 Figure 1: Histogram of the instantaneous clearness index k t for each year.4

Figure 3 :

 3 Figure 3: (a) Principle of the wavelet transform decomposition till level 2 of a global solar irradiance component using the à trous algorithm and the Mexican hat wavelet. (b) First approximation and (c) detail of the decomposition. (d) Second approximation and (e) detail of the decomposition with the threshold (red lines) computed with Starck's algorithm.

  error accuracy) then we go to step 4

Figure 4 :

 4 Figure 4: Flowshart for the threshold value estimation to estimate clear days

Figure 5 :

 5 Figure 5: Left curves: Average number of clear days per month for 2005, 2006, 2007 and 2008 obtained with the two methods. Right curve: Number of clear days per year obtained by the two methods.

Figure 6 :Figure 7 :

 67 Figure 6: Left curve: A recorded solar irradiance component of a clear day (full line) superposed to its fit obtained with Iqbal C model (dashed line). Right curves: Monthly average of the Angstrom coefficient obtained from Iqbal C model built with clear days obtained from clearness index and wavelet methods. The average is over the four years of solar data.

Figure 8 :

 8 Figure 8: Left curves: Monthly average values over the four years of data of the forward scatterance (F c ) considering the IqbalC model and the clear days obtained with by the clearness index criteria and the wavelet method. Right curve: Monthly average of the aerosol single scattering albedo (w 0 ) using the same analysis and four years of data.

Table 1 :

 1 Statistics of the instantaneous clearness index over Ghardaïa city from to 2008.

	Year	2005	2006	2007	2008
	Mean	0.67	0.67	0.67	0.68
	Standard deviation	0.17	0.16	0.17	0.16
	Median	0.73	0.73	0.74	0.74
	25 th Percentile	0.61	0.60	0.61	0.62
	50 th Percentile	0.73	0.73	0.74	0.74
	70 th Percentile	0.78	0.78	0.79	0.78
	Semi-Interquartile	0.09	0.09	0.09	0.08
	deviation
	criterions.

Table 2 :

 2 NB number of clear days per year obtained with the two methods and the corresponding percentage relatively to the whole data set.

	Year	2005	2006	2007	2008
	Clearness index
	NB	134	160	174	165
	Percentage (%)	36.74	43.84	47.67	45.21
	Wavelet method
	NB	113	120	147	110
	Percentage (%)	30.96	32.88	40.27	30.14

Table 3 :

 3 The annual average errors (see text for the notation).

	year	maee	rmsee	maree	rmsree msree	mapee r2	c f	method
	2005	8.786	11.479	0.021	0.035	0.002	2.081	0.995	0.998	Wavelet
		12.196	17.004	0.031	0.060	0.011	3.075	0.991	0.996	index
	2006	10.929	13.696	0.026	0.045	0.004	2.590	0.994	0.998	Wavelet
		15.316	21.532	0.044	0.093	0.026	4.353	0.983	0.992	index
	2007	8.525	11.030	0.020	0.035	0.002	2.010	0.997	0.999	Wavelet
		11.855	16.920	0.031	0.061	0.011	3.080	0.987	0.994	index
	2008	15.519	18.504	0.039	0.060	0.008	3.931	0.957	0.997	Wavelet
		19.247	25.240	0.050	0.095	0.022	5.001	0.946	0.992	index

acknowledgements

We would like to thank and acknowledge the support of the team at the Research Unit on Renewable Energies of Ghardaïa whose collected the data used in the present study.


Philippe Helluy
email: helluy@math.unistra.fr

Erwan Olivier Hurisse
email: olivier.hurisse@edf.fr

Erwan Le Coupanec
email: erwan.lecoupanec@edf.fr

Verification of a two-phase flow code based on an homogeneous model

Keywords: Homogeneous model, two-phase flows, verification

come

Verification of a two-phase flow code based on an homogeneous model

Introduction

The production of electric power through pressurized water nuclear reactor, as many other industrial processes, is associated with the control of the mass transfer phenomena between liquid water and its vapor. The numerical simulation of two-phase flows including phase change is then an important challenge though it is a difficult task. For this purpose, several kinds of models have been proposed for approximately 50 years. The model retained in the sequel has been proposed quite recently [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF]. This two-phase flow model accounts for the complete thermodynamical disequilibrium -namely the pressure, temperature and chemical potential disequilibrium -but makes the assumption that both phases have the same velocity. It is based on the Euler system of equations supplemented by a complex pressure law describing the behaviour of the mixture of the two phases. This pressure law depends on three mixture fractions which define the thermodynamical disequilibrium between the phases, and the return to equilibrium is achieved by three source terms. In the sequel, the verification test cases do not account for these source terms and the focus is put on the convective terms of the model.

In [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF], the behaviour of some explicit first-order schemes has been compared on two Riemann problems, which are unsteady test cases. The aim of the present paper is to extend this verification process to a steady-state analytical test-case and to investigate the behaviour of a first-order semi-implicit scheme, based on a fractional step algorithm [START_REF] Yanenko | Méthode à pas fractionnaires. Résolution de problèmes polydimensionnels de physique mathématique[END_REF]. The latter class of schemes is known to handle more accurately steady-state situations than the first-order explicit schemes. Besides this semi-implicit scheme, the VFRoe-ncv scheme and the Rusanov scheme already used in [START_REF] Buffard | A sequel to a rough Godunov scheme: application to real gases[END_REF] are retained here. In order to improve the accuracy of the VFRoe-ncv scheme on the contact waves, a centered-pressure correction is proposed on the basis of the work [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF].

Two analytical test cases have been retained here. Both are representative of industrial applications at high pressure and temperature. First, as in [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF] the different schemes have been used to compute approximate solutions of a Riemann problem. This unsteady test case allows to test the scheme for breach simulations which are part of the validation test cases for industrial codes in the nuclear domain. The second test case is a rough representation of the way a core or a steam generator operates. It consists in computing the steady state reached by a mixture of liquid and water which is heated by an external source. The solution proposed in [START_REF] Hurisse | A new verification test case for the compressible module of Code Saturne based on a non trivial steady-state solution for the Euler system with energy[END_REF] for compressible single-phase flow is thus extended to two-phase flows without return to equilibrium.

The homogeneous model

In the sequel the two-phase-flow homogeneous model [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF][START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF][START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF] is based on the Euler set of equations. The thermodynamic behaviour of the mixture is ruled by two ingredients: three mixture fractions Y = (α, y, z) , which are respectively the volume fraction, the mass fraction and the energy fraction; and a mixture pressure law P = P(Y, τ, e) defined through the Gibbs relation on the entropy of the mixture. The thermodynamic disequilibrium and the return to equilibrium are taken into account by source terms Γ Y on the fraction equations. The whole model then reads:

                       ∂ ∂t (ρY ) + ∂ ∂x (ρU Y ) = ρΓ Y ∂ ∂t (ρ) + ∂ ∂x (ρU ) = 0 ∂ ∂t (ρU ) + ∂ ∂x ρU 2 + P = 0 ∂ ∂t (ρE) + ∂ ∂x (ρU E + U P ) = 0. ( 1
)
The velocity of the mixture is denoted by U . No kinematic disequilibrium is taken into account, then U also represents the velocity of each phase. We note: ρ = 1/τ the density of the mixture, where τ is the specific volume, E = e + U 2 /2 the specific total energy of the mixture and e the specific internal energy of the mixture. The pressure law P is:

P(Y, τ, e) = α T l P l + (1-α) Tv P v z T l + 1-z Tv , (2)
where the phasic pressures P k = P k (τ k , e k ) and temperatures T k = T k (τ k , e k ) must be given by the user, with e l = ze/y, e

v = (1 -z)e/(1 -y), τ l = ατ /y, τ v = (1 -α)τ /(1 -y) if Y represents the fraction of phase l.
The source terms Γ Y are defined in agreement with the second principle of thermodynamics [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF][START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF][START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF], but they are not considered here for the proposed verification test cases and we set here: Γ Y = 0. In the sequel, a Stiffened Gas equation of states will be used for each phase:

P k (τ k , e k ) = e k τ k (γ k -1) -Π k γ k and T k (τ k , e k ) = e k -Π k τ k C v,k ,
where γ k , Π k and C v,k are parameters. For these Stiffened Gas equations of states, the model ( 1) is hyperbolic if the phasic temperatures are non-negative. Values of the thermodynamical coefficients will be given in section 4. More details on the whole model can be found in [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF][START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF].

Numerical schemes

The overall numerical scheme follows a fractional step algorithm in which convective terms and source terms are treated successively [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF]. Since source terms have been omitted in the sequel, we only describe below the numerical schemes that have been implemented to account for the convective part of the model [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF]. The numerical approximations are obtained by using two different types of first-order finite-volume schemes [START_REF] Eymard | The finite volume method[END_REF]:

• approximate godunov schemes which are explicit,

• and a semi-implicit fractional step scheme.

Considering the first class of schemes, three different schemes have been investigated. As proposed in [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF] a VFRoe-ncv scheme [START_REF] Buffard | A sequel to a rough Godunov scheme: application to real gases[END_REF] is used to compute the numerical fluxes at the interface between two cells. The linearised problem at the interface is solved by considering the variables (Y, τ, U, P ). This algorithm is supplemented with the entropic correction proposed in [START_REF] Helluy | A simple parameterfree entropy correction for approximate Riemann solvers[END_REF]. Moreover, a partial WFRoe scheme [START_REF] Gallouët | A new approximate Godunov scheme with application to dense gas-solid flows[END_REF] is performed to improve the prediction of the speed of the contact wave in the linearised problem at the interface between two cells. Otherwise, non-symmetric behaviours have been observed on some test cases, in particular for bubble collapse. This scheme will be nicknamed here VFROENCV. In order to improve the accuracy of this scheme, another version has been implemented on the basis of the idea of [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF]. The interfacial pressure issued from the linearised Riemann problem is mixed with a centered pressure chosen as the mean of the pressure of the left (L) and right (R) cells. The blending function f b (M ) for the pressure is based on the local Mach number M = |u LR |/max(C L , C R ), where u LR is the interfacial WFRoe-ncv velocity. C L and C R are respectively the left and right sound speeds. We have chosen the function:

f b (M ) = 1-cos(ΠM/M lim ) 2 , if M ≤ M lim , 1, otherwise. (3)
The modified interfacial pressure P * LR used to compute the numerical fluxes is then:

P * LR = f b (M )P LR + (1 -f b (M )) P L + P R 2 ;
where P LR is the interfacial pressure computed by the VFROENCV scheme, P L and P R the left and right pressures. The parameter M lim is the limit mach number above which the pressure correction is not activated. The influence of this parameter on the results will not be discussed here even if it is an important parameter. The results of section 4 have been obtained with M lim = 0.25. This scheme will be nicknamed here VFROENCV-C. More details about VFROENCV and VFROENCV-C schemes can be found in appendix B. The Rusanov scheme [START_REF] Rusanov | Calculations of interaction of non-steady shock waves with obstacles[END_REF] is also considered as it remains a reference scheme, it will be denoted by RUSANOV.

The VFRoe-ncv scheme is an explicit scheme. In order to maintain the stability of the numerical approximations a CFL condition must be fulfilled by the time-step. Moreover, it is known that for low mach-number situations discrepancies may be observed. These limitations can be strong drawbacks for industrial applications.

The semi-implicit scheme [START_REF] Martin | Modeling fluid flows in obstructed media[END_REF] relies on a fractional step method [START_REF] Yanenko | Méthode à pas fractionnaires. Résolution de problèmes polydimensionnels de physique mathématique[END_REF] which consists in 4 main steps: a density step, a momentum step, a total energy step and a mixture fraction convection step. The density step consists in solving the continuity equation turned into a pressure equation (thanks to a linearization of the density). This step yields a predicted pressure and allows to update the mass flux and the density. The momentum step then consists in solving the momentum balance equation considering the explicit gradient of the predicted pressure and the updated mass flux. The total energy step finally consists in solving the total energy balance equation considering the updated mass flux. The internal energy can at last be deduced and is taken into account together with the updated density (from the density step) in order to update the pressure. This scheme is detailed in appendix C. It will be nicknamed here SEMI-IMP.

Numerical results

The behavior of the schemes is compared on two different test cases: an unsteady test case and a steady-state test case. For both tests we use the same parameters for the phasic equation of states: The approximations obtained by the schemes on uniform meshes have been compared with the analytical solution in terms of the L 1 error:

Cv, v = 4.477815802223535 10 3 J/kg/K Cv, l = 1.395286166711847 10 3 J/kg/K γ v = 1.084875362318841 γ l = 1.665128030303030 Π v = 0.0 Π l = 3.
err(φ, T ) = N i=1 |φ approx (T, x i ) -φ exact (T, x i )| N i=1 |φ exact (T, x i )| , (5)
for the variables φ = {α, y, z, ρ, U, P }. This error is computed at the final time t for the unsteady test case and when convergence in time is assumed to be reached for the steady case (see section 4.2 for details on the convergence criterium).

A Riemann problem test case

As in [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF], we consider here a Riemann problem consisting in a U -contact wave and a U + C-shock wave. The U -C-wave is a ghost wave. The domain is [0, 1 m] and the initial discontinuity is located at schemes: RUSANOV, VFROENCV, VFROENCV-C and SEMI-IMP. The meshes have from 50 to 200000 regular cells, except for the SEMI-IMP scheme (up to 50000 cells). The error defined in ( 5) is plotted on figures 3 and 4 and the different approximate solutions are compared on figure 5 for a mesh with 800 cells and a time step chosen such that CF L U +C = 0.5.

x d = 1/2 m.
Considering figure [START_REF] Chorin | Numerical solutions of the Navier-Stokes equations[END_REF], it can be noted that all the schemes have almost the same accuracy on coarse meshes, except the RUSANOV scheme that is less accurate for α and τ . Nonetheless, when the mesh is refined, the VFROENCV-C scheme widens the gap with the other schemes for the variable P , U and τ . Up to 200000 cells the asymptotic order of convergence 1/2 is not yet reached with VFROENCV-C for P and U . On the last three meshes the convergence rate is close to 1 for P and U , which indicates that the scheme is very accurate on the contact wave. The SEMI-IMP scheme with a time step chosen such that CF L U +C = 0.5 is close to the VFROENCV scheme. Moreover, figure [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF] shows that increasing the time-step for SEMI-IMP leads to a loss of accuracy.

Besides the error analysis, one can observe on figure (5) that the schemes have very different behaviours when a discontinuity is approximated, see figure [START_REF] Ferrand | Une formulation intégrale pour la modélisation d'écoulements fluides en milieu encombré[END_REF].

• The U -C ghost wave (on the left on figure [START_REF] Ferrand | Une formulation intégrale pour la modélisation d'écoulements fluides en milieu encombré[END_REF]).

From an analytical point of view all the variables remain constant when passing through the wave U -C, but for all the tested schemes a "numerical wave" U -C can be observed. These discrepancies disappear when the mesh is refined, but they can be noticeable on coarse meshes. The VFROENCV-C scheme shows very strong oscillations with respect to the other schemes, and the SEMI-IMP scheme at CF L U +C = 0.5 is very close to the VFROENCV scheme which remains quite accurate.

• The contact wave (in the middle on figure [START_REF] Ferrand | Une formulation intégrale pour la modélisation d'écoulements fluides en milieu encombré[END_REF]). Through the contact wave, pressure and velocity are constant in the analytical solution. The VFROENCV scheme presents small discrepancies on this wave. The accuracy is then very good despite the fact that the pressure law does not fulfil the requirement to ensure the preservation of P and U as exposed in [START_REF] Gallouët | Some recent finite volume schemes to compute Euler equations using real gas EOS[END_REF]. The VFROENCV-C scheme generates oscillations on both sides of the contact wave, whereas the SEMI-IMP scheme presents a pressure over-shoot combined with small oscillations.

• The shock wave (on the right on figure [START_REF] Ferrand | Une formulation intégrale pour la modélisation d'écoulements fluides en milieu encombré[END_REF]).

Both SEMI-IMP and VFROENCV-C approximations have oscillations at the shock location, but their order of magnitude is 100 times less for the former than for the latter. They are thus almost negligible for the SEMI-IMP scheme in comparison with the pressure jump through the schock wave (∼ 36 10 5 P a).

Zoom for the pressure : ghost wave (left), contact wave (middle), shock wave (right).

0,62 0,64 0,66 0,68

A steady state test case

This test case consists in heating a mixture of water and vapor, assuming that the source terms for the fractions are equal to zero. The heating of the mixture is ensured by an external source term in the energy equation. The computation of the corresponding steady state for system (1) is reported in appendix A. Due to the complex mixture equation of state of the model, the solution at each x is obtained by the mean of a Quasi-Newton algorithm. The boundary conditions are Dirichlet conditions filled with the exact solution. The initial conditions are also the exact solution for all x. The iterations are stopped when the relative L 1 -norm between the variables Ψ ∈ {Y, τ, U, P } at iteration n and (n -1) is less than 10 -12 :

∀Ψ ∈ {Y, τ, U, P }, N i=1 |Ψ approx (t n , x i ) -Ψ approx (t n-1 , x i )| N i=1 |Ψ approx (t n-1 , x i )| < 10 -12 .
The domain is x ∈ [0, 8 m]. The meshes contain from 60 to 12800 cells, except for the SEMI-IMP scheme. Due to the computation time, the finer mesh contains 6400 cells. The source term Φ is applied for 2 m < x < 6 m and it is uniform within this interval:

Φ(x) = 5.0 10 9 W/m 3 , if 2 m < x < 6 m, 0, elsewhere. (6)
At the inlet x = 0 the reference state is: ρ 0 = 728.0 kg/m 3 ; U 0 = 5.0 m/s; P 0 = 7.69 10 6 P a; α 0 = 1.38006860749359 10 -1 ; y 0 = 7.84348344805274 10 -3 ; z 0 = 1.51281566759212 10 -2 ; [START_REF] Martin | Modeling fluid flows in obstructed media[END_REF] The analytical solution, which is computed as explained in appendix A, is plotted on figure 7. From a numerical point of view, the heating source term Φ is integrated using a fractional step approach [START_REF] Yanenko | Méthode à pas fractionnaires. Résolution de problèmes polydimensionnels de physique mathématique[END_REF]. The convection terms are first computed and the ODE system including Φ is then integrated by an explicit Euler scheme. For the present steady-state test case with a constant source term, a well-balanced approach could be helpful to improve the accuracy of the VFROENCV scheme. But we are also interested in non constant source terms and unsteady situations for which the fractional step approach remains efficient [START_REF] Gallouët | Well-balanced schemes versus fractional step method for hyperbolic systems with source terms[END_REF]. We thus focus here on SEMI-IMP and VFROENCV-C schemes which improve the prediction of the steady-state solution as shown on figures (8), ( 9), ( 10), [START_REF] Godlewski | Numerical analysis for hyperbolic systems of conservation laws[END_REF] and [START_REF] Helluy | A simple parameterfree entropy correction for approximate Riemann solvers[END_REF].

It should be first mentioned that for all schemes the error between the approximated fractions and the analytical solutions is equal to the round-off errors and is thus not considered in the sequel. The convergence curves are plotted on figure [START_REF] Gallouët | A new approximate Godunov scheme with application to dense gas-solid flows[END_REF]. One can note that on coarse meshes VFROENCV-C is less accurate than VFROENCV for the velocity variable, but more accurate for the pressure. As for the problem of section 4.1, the VFROENCV-C scheme still shows an order of convergence of the pressure and velocity greater than the asymptotic order. Due to the discontinuity of the function Φ, a stationary contact wave is introduced in the system of equations, which leads to an asymptotic order of convergence of 1/2. The latter is recovered on almost all the variables of figure [START_REF] Gallouët | A new approximate Godunov scheme with application to dense gas-solid flows[END_REF], except for the SEMI-IMP scheme for which computations should be pursued on finer meshes.

The results for the different schemes have been plotted on figure (8) for a mesh with 400 cells. The RUSANOV and VFROENCV schemes introduce unphysical pressure gaps at the inlet and outlet of the heating zone (where the source term is discontinuous). These gaps disappear when the mesh is refined (see figure [START_REF] Gallouët | Well-balanced schemes versus fractional step method for hyperbolic systems with source terms[END_REF]), but on industrial meshes they can lead to unexpected vaporisation/condensation when dealing with mass transfer. The SEMI-IMP (CF L U +C = 25) and the VFROENCV-C schemes do not present such important pressure gaps, but oscillations appear arround the location of the sudden change of the heating source term. The amplitude of these oscillations is smaller for the SEMI-IMP scheme, but a more important error is observed in the bulk of the heated zone when compared with the VFROENCV-C scheme. This explains why the error on the pressure variable is greater with the SEMI-IMP scheme compared with the error obtained for the VFROENCV-C (see figure [START_REF] Gallouët | A new approximate Godunov scheme with application to dense gas-solid flows[END_REF]).

Conclusion

Four schemes have been tested on two different verification test-cases, which are sketches of industrial configurations. Though the VFROENCV-C scheme presents the best results in terms of the error analysis, the oscillations that appear where discontinuities are computed could be a drawback for industrial purposes. They could for instance be confusing in the case of a configuration where physical oscillations might naturally appear. There are also oscillations in the results obtained with the SEMI-IMP scheme, but their amplitude is smaller on the proposed test cases. On the contrary, the VFROENCV scheme does not generate oscillations, but the main drawback is that it is globally less accurate than the latter two schemes. Some additional investigations should be performed. In particular a quantitative comparison of the error with respect to the CPU-time would be worthwhile. Actually, if the SEMI-IMP and VFROENCV-C schemes have lower approximation errors for a given mesh size with respect to VFROENCV and RUSANOV, the CPU-time is far greater on the steady-state case of section 4.2. This difference of CPU-time is less important for the unsteady test of section 4.1. Some additional comparisons could also be carried out:

(i) different values of the parameter M lim of the VFROENCV-C scheme ;

(ii) different values of the time-step for the SEMI-IMP scheme ;

(iii) a regular source term Φ ;

(iv) uniform initial conditions and an outlet boundary conditions at x = 8.

International Journal on Finite Volumes 6 Appendix A: A class of analitycal steady-state solutions for compressible two-phase heated flows

In the following, we extend the solution obtained for single-phase flows in [START_REF] Hurisse | A new verification test case for the compressible module of Code Saturne based on a non trivial steady-state solution for the Euler system with energy[END_REF] for two-phase flows.

We restrict ourselves to the steady-state solutions of (1) without the source terms S Y . The variables thus only depend on the space variable x and their space evolution follows the ODE system:

                     d dx (ρU Y ) = 0 d dx (ρU ) = 0 d dx ρU 2 + P = 0 d dx (U (ρE + P )) = Φ, (8)
where the thermodynamical behaviour of the mixture is described by a law linking Y , ρ, e and P . We set here e = (Y, ρ, P ). The source term Φ(x) denotes an external source of heat and it only depends on the space variable x.

The second equation of [START_REF] Gallouët | Some recent finite volume schemes to compute Euler equations using real gas EOS[END_REF] states that the mass flow rate ρU is uniform, so that we denote it by D 0 = ρU . The system (8) can thus be written:

                     D 0 d dx (Y ) = 0 U = τ D 0 D 0 d dx (U ) + d dx (P ) = 0 D 0 d dx (E) + d dx (U P ) = Φ. (9)
We assume that D 0 = 0, so that the fractions Y are uniform: Y = Y 0 . Thanks to the second equation of ( 9), we replace U by τ in the third and fourth equation of [START_REF] Gallouët | A new approximate Godunov scheme with application to dense gas-solid flows[END_REF]. We get the two equations:

     D 2 0 d dx (τ ) + d dx (P ) = 0 D 0 d dx (e) + D 3 0 2 d dx τ 2 + D 0 d dx (τ P ) = Φ. (10)
International Journal on Finite Volumes

Since the mass flow rate D 0 is uniform and since the source term Φ only depends on x, we can integrate this subsystem as follows on [x 0 , x]: D 2 0 (τ -τ 0 ) + (P -P 0 ) = 0 D 0 (e -e 0 ) +

D 3 0 2 (τ 2 -τ 2 0 ) + D 0 (τ P -τ 0 P 0 ) = x x 0 Φ(s)ds, (11)
where a subscript 0 denotes the value at x 0 which must be given by the user. The symbols without subscript denote the values at a point x. Due to our choice for the equation of states, we have: e = (Y, ρ, P ) and e 0 = (Y 0 , ρ 0 , P 0 ).

The first equation of [START_REF] Godlewski | Numerical analysis for hyperbolic systems of conservation laws[END_REF] gives explicitly τ as a function of P and since we have Y = Y 0 , the second equation of ( 11) is a non linear equation for P :

τ = τ 0 -(P -P 0 ) D 2 0 D 0 ( (Y 0 , 1 τ , P ) -e 0 ) - (P 2 -P 2 0 ) 2D 0 = x x 0 Φ(s)ds. (12)
Hence, solving this equation for P allows to compute the value for τ (with the first equation of ( 12)), and then, since D 0 = ρ 0 U 0 = ρU , we get U . Obviously, depending on the equation of state , Φ and D 0 , the problem can have one, several or no solution. Moreover it can be complex to solve and it may require the use of numerical algorithm to compute P at each x.

7 Appendix B: The VFROENCV-C and VFROENCV scheme

The VFROENCV and the VFROENCV-C schemes rely on the VFRoe-ncv scheme using the variable (Y, τ, U, P ) [START_REF] Buffard | A sequel to a rough Godunov scheme: application to real gases[END_REF]. Let us first recall the latter algorithm.

For the sake of simplicity, we restrict to regular meshes of size ∆x such that:

∆x = x i+ 1 2 -x i-1 2
, i ∈ Z. We denote ∆t the time step, where ∆t = t n+1 -t n , n ∈ N. In order to approximate solutions of the exact solution W ∈ R p of the conservative hyperbolic system:

∂ ∂t (W ) + ∂ ∂x (F (W )) = 0 W (x, 0) = W 0 (x)
with F (W ) in R p . Let W n i be the approximate value of 1 ∆x

x i+ 1 2 x i-1 2 W (x, t n )dx.
Integrating over [x i-1 2

; x i+ 1 2
] × [t n ; t n+1 ] provides:

W n+1 i = W n i - ∆t ∆x φ n i+ 1 2 -φ n i-1 2
The numerical flux φ n i+

     ∂ ∂t (Z) + B( Ẑ) ∂ ∂x (Z) = 0 Z(x, 0) = Z L if x < 0 Z R otherwise (13)
with the initial condition: Z L = Z(W i ) and Z R = Z(W i+1 ). The matrix:

B(Z) = (W ,Z (Z)) -1 A(W (Z)) W ,Z (Z) (A(W )
is the Jacobian matrix of flux F (W )). Once the exact solution Z * ( x t ; Z L , Z R ) of this approximate problem ( 13) is obtained, the numerical flux is defined as:

φ(W L , W R ) = F (W (Z * (0; Z L , Z R )))
Let us set l k , λ k and r k , k = 1, ..., p, left eigenvectors, eigenvalues and right eigenvectors of matrix B( Ẑ) respectively. The solution Z * ( x t ; Z L , Z R ) of the linear Riemann problem is defined everywhere (except along x t = λ k ):

Z * x t ; Z L , Z R = Z L + x t > λ k ( t l k .(Z R -Z L )) r k = Z R - x t < λ k ( t l k .(Z R -Z L )) r k
or in a slightly different form:

Z R -Z L = k=1,p ( t l k .(Z R -Z L )) r k = k=1,p αk r k (14)
setting:

αk = t l k .(Z R -Z L )
The only remaining unknown is the mean Ẑ which must comply with the condition:

Ẑ(Z l = Z 0 , Z r = Z 0 ) = Z 0
The standard average which is used is:

Ẑ(Z L , Z R ) = (Z L + Z R )/2
The explicit form of the approximate Godunov scheme will be under conservative form:

W n+1 i -W n i + ∆t ∆x F (W (Z * (0; Z n i , Z n i+1 ))) -F (W (Z * (0; Z n i-1 , Z n i ))) = 0 Remark 7.1 A different prediction may be obtained using instead: Ẑ(Z L , Z R ) = Z * (0; Z L , Z R )
where the approximate value at the interface Z * (0; Z L , Z R ) is obtained solving [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF] with:

Ẑ(Z L , Z R ) = (Z L + Z R )/2
This corresponds to the WFRoe-ncv scheme [START_REF] Gallouët | A new approximate Godunov scheme with application to dense gas-solid flows[END_REF].

The two shemes VFROENCV and VFROENCV-C proposed in the sequel are based on the variable change Z = (Y, τ, U, P ). Moreover, the velocity value of the average state Ẑ(Z L , Z R ) = (Z L + Z R )/2 is replaced by the value of the velocity U * (0; Z L , Z R ) which is computed through [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF] as in the WFroe-ncv scheme.

At last, for the VFROENCV-C, the pressure P * (0; Z L , Z R ) obtained by solving (13) is blended with the centred pressure (P L + P R )/2 as proposed in [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF]. The blending function f b (M ) for the pressure is based on the local Mach number

M = |U * (0; Z L , Z R )|/max(C L , C R ),
where C L and C R are respectively the left and right sound speeds. We have chosen the function:

f b (M ) = 1-cos(ΠM/M lim ) 2 , if M ≤ M lim , 1, otherwise.
The modified interfacial pressure P * LR used to compute the numerical fluxes is then:

P * LR = f b (M )P * (0; Z L , Z R ) + (1 -f b (M )) P L + P R 2 ;
where P L and P R are the left and right pressures. The parameter M lim is the limit mach number above which the pressure correction is not activated.

In fact, the momentum Q n+1 is approximated by the momentum Q * obtained through the following simplified momentum balance equation:

∂ t Q + ∇P (ρ, s) = 0, (15)
which yields:

Q * = Q n -∆t n ∇P n+1,-. (16)
The finite volume scheme becomes:

vol (Ω) ρ n+1,--ρ n + ∆t n Γ Q * • ndΓ = 0.
It is then assumed that, for a constant entropy, we have:

δP = c 2 n δρ, (17)
with c 2 n = c 2 (P n , ρ n ) (see remark 8.1). This linearization is then used to obtain the form:

vol Ω φ 1 (c 2 ) n (P n+1,--P n ) + ∆t n Γ Q * .ndΓ = 0. ( 18
)
Since Q * depends on P n+1,-through equation ( 16), the pressure P n+1,-can be computed from equation [START_REF] Temam | On an approximate solution of Navier-Stokes equations by the method of fractional step[END_REF]. This step then allows to obtain P n+1,-, from which ρ n+1,-can be deduced using the relation [START_REF] Rusanov | Calculations of interaction of non-steady shock waves with obstacles[END_REF]. This update of the density is mandatory in order to ensure the mass conservation. ∂s . In the mass balance step, it is thus assumed that the entropy remains constant.

Momentum balance step

In this second step, the velocity is computed implicitely, while the density and the pressure are known from the previous step; the total energy and the mixture fractions remain steady here. The momentum balance equation is integrated in time between t n and t n+1 and in space over Ω:

t n+1 t n Ω ∂Q ∂t + ∇ • u ⊗ Q + ∇P dΩdt = 0.
As for the mass balance step, the different terms of the equations are treated separately. By applying Green formula and using an upwind scheme for the convective term, the following time scheme is obtained:

vol (Ω) Q n+1,--Q n + ∆t n Γ u Q.n n+1,-dΓ + ∆t n Γ P n+1,-ndΓ = 0.
This step allows to obtain u n+1,-, and thus Q n+1,-= ρ n+1,-u n+1,-.

Remark 8.2 It is possible to rewrite

Q n+1,--Q n = ρ n (u n+1,--u n ) + u n+1,-(ρ n+1,--ρ n )
considering that the increment (ρ n+1,--ρ n ) is known from the previous step.

Energy balance step

In this third step, the total energy is computed with a semi-implicit scheme, while the velocity, the pressure and the density are known from the previous steps. The mixture fractions Y still remain constant in this step. The energy balance equation (written in terms of the total energy) is integrated in time between t n and t n+1 and in space over Ω:

t n+1 t n Ω ∂E ∂t + ∇ • (u(E + P )) -Φ dΩdt = 0,
where Φ stands for the external energy source term.

Once again, the terms of the equation are treated separately and Green formula is applied. Using an upwind scheme for the convective term, the following time scheme is deduced:

vol (Ω) (E n+1,--E n ) + ∆t n Γ (Q * • n) E + P ρ n+1,- dΓ (19)
-∆t n vol (Ω) Φ n+1,-= 0.

We use the mass flux Q * computed during the mass balance step in order to ensure the consitency with the mass balance equation.

This step allows to obtain the variable E n+1,-.

Convection of the mixture fractions Y

In this fourth step, the mixture fractions Y are simply convected using an upwind scheme with the mass flux Q * computed during the mass balance step:

vol (Ω) ρ n+1,-(Y n+1,--Y n ) + ∆t n Γ (Q * • n)Y n+1,-dΓ = 0.
We use the mass flux Q * to ensure the consitency with the mass balance equation.

This step allows to obtain the variable Y n+1,-.

Update of the variables at t = t n+1

At the end of each time step, the variables are updated from t n+1,-to t n+1 as follows:

ρ n+1 = ρ n+1,- u n+1 = u n+1,- E n+1 = E n+1,-
Y n+1 = Y n+1,- P n+1 = P (Y n+1 , ρ n+1 , n+1 )

where

: n+1 = E n+1 ρ n+1 - 1 2 (u 2 ) n+1

725876146842836 10 8 P a ( 4 )

 4 These parameters have been computed to be representative of steam-water flows around 81 bars and at saturation (T sat (81 bars) = 569 K).

Figure 1 :

 1 Figure 1: States for the first Riemann problem.

Figure 2 : 5 Figure 3 : 5 Figure 4 :Figure 5 :

 253545 Figure 2: Analytical solution at time t = 2. 10 -4 s for the Riemann problem of section 4.1.

Figure 6 :

 6 Figure6: Comparison of the approximate solutions on a mesh with 800 cells, zoom at the wave locations for VFROENCV (red), VFROENCV-C (black), SEMI-IMP at CF L U +C = 0.5 (blue). The analytical solution is in green.

Remark 4 . 1

 41 The reference state corresponds to the left states of the Riemann problem in section 4.1.

Figure 7 :

 7 Figure 7: Analytical solution for the steady state heating case.

Figure 8 :

 8 Figure8: Comparison of the approximated solutions for a mesh with 400 cells. The CF L U +C is equal to 0.5 for the explicit schemes and to 25 for the SEMI-IMP scheme.

/ 2 Figure 9 :

 29 Figure 9: Convergence curves for the steady-state problem: ln() of the L 1 error norm (5) with respect to ln(1/N ) (N number of cells).The CF L U +C is equal to 0.5 for the explicit schemes and to 25 for the SEMI-IMP scheme.

Figure 12 :

 12 Figure 12: Comparison of the approximated pressure for SEMI-IMP on different meshes and for a time-step corresponding to CF L U +C = 25.

Remark 8 . 1

 81 Actually, for non constant entropy s, the relation is dP = c 2 dρ + βds with β = ∂P (ρ,s)

  Figure 10: Comparison of the approximated pressure for VFROENCV on different meshes. The CF L U +C is equal to 0.5.

	7,7e+06				60 cells		7,7e+06				800 cells
					100 cells							1600 cells
					200 cells							3200 cells
					400 cells							6400 cells
					800 cells							12800 cells
	7,68e+06				sol. exacte	7,68e+06				sol. exacte
	7,66e+06
						7,66e+06
	7,64e+06
						7,64e+06
	7,62e+06
						7,62e+06
	7,6e+06
	0	2	4	6	8	7,6e+06	0	2		4	6	8
	7,7e+06
					60 cells							800 cells
					100 cells						400 cells	1600 cells
	7,9e+06 7,68e+06				200 cells 400 cells 800 cells sol. exacte		7,7e+06				800 cells 1600 cells 3200 cells 6400 cells	3200 cells 6400 cells 12800 cells sol. exacte
											analytical
						7,68e+06
	7,66e+06
	7,8e+06
						7,66e+06
	7,64e+06
	7,7e+06					7,64e+06
	7,62e+06
						7,62e+06
	7,6e+06	0		2			7,6e+06 4	0	2	6	4	6	8 8

Figure 11: Comparison of the approximated pressure for VFROENCV-C on different meshes. The CF L U +C is equal to 0.5.

International Journal on Finite Volumes

Appendix C: The SEMI-IMP scheme

The algorithm detailed below is based on a fractional step method [START_REF] Yanenko | Méthode à pas fractionnaires. Résolution de problèmes polydimensionnels de physique mathématique[END_REF] and follows the classical idea of the so-called pressure correction schemes [START_REF] Chorin | Numerical solutions of the Navier-Stokes equations[END_REF][START_REF] Patankar | A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows[END_REF][START_REF] Temam | On an approximate solution of Navier-Stokes equations by the method of fractional step[END_REF]. Each time step from t n to t n+1 is divided into five steps:

(i) a mass balance step allowing to update the density and to predict an approximation of the pressure;

(ii) a momentum balance step allowing to update the velocity;

(iii) an energy balance step allowing to update the total energy;

(iv) a step for the pure convection of the mixture fractions Y = {α, y, z};

(v) and source term step that accounts for the return to the thermodynamical equilibrium.

For each one of the above five steps (i) to (v), the time derivatives are discretized using an implicit Euler scheme. Moreover, the scheme used for the source term step (v) is exactly the same than that used with all the other schemes of the paper [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF].

The time step is defined as ∆t n = t n+1 -t n . By abuse of notation, for any variable Ψ, Ψ n = Ψ(t n ). Moreover Ψ n+1,-stands for the estimation of the variable Ψ during the fractionnal step, and we set δΨ = Ψ n+1,--Ψ n . All the variables are then updated (from t n+1,-to t n+1 ) once the source term step (v) is achieved.

In the following, we set Q = ρU , E = ρE, and we present the convective part (i.e. steps (i) to (iv)) of the scheme denoted by SEMI-IMP in the paper. It should be mentionned that this algorithm is conservative and it allows then to retrieve the correct shock solutions.

Mass balance step

In this first step, pressure and density are implicit, while the velocity, the mixture fractions Y and the entropy are considered frozen at time t n . The convective mass flux is computed as well, it will be noted Q * in order to distinguish it from the mass flux computed at the next step. The mass balance equation is integrated in time between t n and t n+1 and in space over Ω:

The two terms of the above equation are treated separately. The divergence term is implicit and Green formula is applied:

International Journal on Finite Volumes