Abstract-Potentiometry with ion-selective electrodes (ISEs) provides a simple and cheap approach for estimating ionic activities. However, a well-known shortcoming of ISEs regards their lack of selectivity. Recent works have suggested that smart sensor arrays equipped with a blind source separation (BSS) algorithm offer a promising solution to the interference problem. In fact, the use of blind methods eases the time-demanding calibration stages needed in the typical approaches. In this work, we develop a Bayesian source separation method for processing the outputs of an ISE array. The major benefit brought by the Bayesian framework is the possibility of taking into account some prior information, which can result in more realistic solutions. Concerning the inference stage, it is conducted by means of Markov chain Monte Carlo (MCMC) methods. The validity of our approach is supported by experiments with artificial data and also in a scenario with real data.

Index Terms-Blind source separation, ion-selective electrode, chemical sensor array, Bayesian approach.

I. INTRODUCTION

I ON-selective electrodes (ISEs) [2], [3], [4] are devices used for estimating the ionic activity, a measure of effective concentration of an ion in aqueous solution. The transducer mechanism in an ISE is based on a sensitive membrane in which the electrochemical potential depends on the concentration of a given charged chemical specie. A typical example of ISE is the glass electrode used for measuring the pH value. Moreover, ISEs for ions such as ammonium, potassium, and calcium are of great interest in applications like water quality control [START_REF] Sakai | Evaluation of water quality and pollution using multichannel sensors[END_REF] and biomedical monitoring [START_REF] Oesch | Ion-selective membrane electrodes for clinical use[END_REF], [START_REF] Guenat | Development of an array of ion-selective microelectrodes aimed for the monitoring of extracellular ionic activities[END_REF].

The most relevant drawback of an ISE concerns its lack of selectivity, that is, an ISE may respond to interfering ions other than the target one. While in some situations the interference process can be negligible, there are other cases in which this problem seems to be more accentuated, especially if the ions under analysis have similar physicochemical properties. If no care is taken in these cases, the measurements provided by the ISE may become uncertain.

A possible way to overcome the interference problem relies on the concept of smart sensor array (SSA). In a SSA, the data acquisition is done through several ISEs. Then, the desired information is extracted by a data processing block that tries to exploit the diversity provided by the array. Typically, the Leonardo T. Duarte and C. Jutten are with the GIPSA-lab (UMR CNRS 5216), Institut Polytechnique de Grenoble, Grenoble, France (e-mail: leonardo.duarte@gipsa-lab.inpg.fr, christian.jutten@gipsa-lab.inpg.fr). C. Jutten is also with Institut Universitaire de France.

Saïd Moussaoui is with IRCCyN (UMR CNRS 6597), Ecole Centrale Nantes, Nantes, France (e-mail:said.moussaoui@irccyn.ec-nantes.fr).

A preliminary version of this paper [START_REF] Duarte | Ion-selective electrode array based on a bayesian nonlinear source separation method[END_REF] was presented at the 8th International Conference on Independent Component Analysis and Signal Separation (ICA 2009) and was published in its proceedings. development of a SSA follows a supervised paradigm, i.e. a set of training (or calibration) samples is used to adjust the parameters of the SSA data processing block. Despite the good results obtained in both quantitative and qualitative analysis (for example, see the electronic tongues [START_REF] Vlasov | Electronic tongue: chemical sensor systems for analysis of aquatic media[END_REF]), supervised methods suffers from two practical problems: 1) the acquisition of a sufficient number of training samples may be time-demanding, and 2) the calibration procedure must be performed from time to time due to the sensor's drift. These limitations are the main motivations for the use of unsupervised (or blind) methods. In this new case, the SSA data processing block only uses the array responses and some parametric knowledge on the interference process. Therefore, the calibration stage may be skipped or, at least, quite reduced.

In this work, we are interested in the problem of unsupervised quantitative analysis via an ISE array, which can be expressed as follows: there is a set of source signals, each one representing the activity of a given ion. Due to the interference phenomena in the transduction stage, the signals recorded by the ISE array are mixed versions of the sources. Then, our goal is to retrieve the sources using only the ISE array outputs. This description is an example of a Blind Source Separation problem (BSS) [START_REF] Hyvärinen | Independent component analysis[END_REF], [START_REF] Comon | Séparation de sources 1 : concepts de base et analyse en composantes indépendantes[END_REF], to which the signal processing community has been devoting a great deal of attention. The main difficulty in the application of BSS methods to the problem treated in this work stems from the nonlinear behavior of an ISE. While linear BSS is now a well established matter, the more general problem of nonlinear BSS is still a very challenging one (see [START_REF] Jutten | Advances in blind source separation (BSS) and independent component analysis (ICA) for nonlinear mixtures[END_REF]).

The first studies [START_REF] Bermejo | ISFET source separation: foundations and techniques[END_REF], [START_REF] Bedoya | Improving semiconductor-based chemical sensor arrays using advanced algorithms for blind source separation[END_REF], [START_REF] Duarte | A mutual information minimization approach for a class of nonlinear recurrent separating systems[END_REF] on the design of BSS methods for chemical sensor arrays were based on Independent Component Analysis (ICA) methods [START_REF] Hyvärinen | Independent component analysis[END_REF]. Despite the encouraging results presented in these works, the developed methods have several limitations. Firstly, they do not take into account some prior information that would be useful for improving the separation quality. For instance, the sources in our problem are always non-negative, since they represent ionic activities. Secondly, given that ICA methods works with the assumption that the sources are statistically independent, they cannot be used when the sources are correlated. Finally, all these works were evaluated with synthetic data in noiseless situations.

In this paper, we propose a Bayesian source separation method for processing the outputs of an ISE array. The Bayesian approach has been proved to be attractive in several domains including audio source separation [START_REF] Févotte | A Bayesian approach for blind separation of sparse sources[END_REF] and spectroscopy data analysis [START_REF] Moussaoui | Separation of non-negative mixture of non-negative sources using a bayesian approach and mcmc sampling[END_REF]. Moreover, one of the main advantages of this approach is exactly the possibility of incorporating prior information into our method and of working in noisy scenarios. Finally, the advent of sampling techniques [START_REF] Robert | Monte Carlo Statistical Methods[END_REF] such as Monte Carlo Markov Chain (MCMC) methods has provided an efficient solution to the calculation of the complex integrals that are usually associated with a Bayesian framework.

This work is an extended version of the conference paper [START_REF] Duarte | Ion-selective electrode array based on a bayesian nonlinear source separation method[END_REF], which only contained the basics of our method and some preliminary results. In the present version, we detail important aspects of our method and widely expand the results section, including a comparison with a supervised technique. Concerning the paper's organization, we start with a description of the mixing model associated with the application in mind. After that, in Section III, we describe the building blocks of our Bayesian source separation scheme. In Section IV, we conduct some experiments with artificial and real data to assess the performance of our proposal. Finally, our conclusions are presented in Section V.

II. PROBLEM STATEMENT

Let us consider a solution containing n s ions and an array composed of n c electrodes. Due to the interference problem, the response of the i-th ISE within the array at the instant t, represented by y it , is dependent not only on the activity of its target ion, s it , but also on the activities of the other ions in the solution, s jt . According to the Nicolsky-Eisenman (NE) equation [START_REF]Microcapteurs chimiques et biologiques[END_REF], this phenomenon can be modeled as1

y it = e i + d i log   s it + n s j=1,j =i a ij s z i /z j jt   , (1)
where e i is a constant, z j is the valence of the j-th ion and a ij denotes the selectivity coefficients. Finally,

d i = RT /z i F ,
where R is the gas constant, T the temperature in Kelvin, and F the Faraday constant. When the valences of the ions under analysis are the same, the model (1) becomes a particular case of the class of post-nonlinear (PNL) models [START_REF] Bedoya | Improving semiconductor-based chemical sensor arrays using advanced algorithms for blind source separation[END_REF]. Nonetheless, when the valences are different and the ISEs within the array have different ions as targets, the resulting model becomes tougher due to the appearance of a nonlinearity inside the logarithm term [START_REF] Duarte | A mutual information minimization approach for a class of nonlinear recurrent separating systems[END_REF]. By arranging the response of each ISE within the array, the mixing process can be expressed as follows:

Y = e • 1 1×n d + diag(d) log (A ⊗ z S) , ( 2
)
where n d corresponds to the number of available samples,

Y ∈ R n c ×n d , e = [e 1 , . . . , e n c ] T , d = [d 1 , . . . , d n c ] T , A ∈ R n c ×n s + and S ∈ R n s ×n d + .
The element it of Y corresponds to y it and therefore the i-th row of Y corresponds to the time response of the i-th ISE within the array. Analogously, the element jt of S denotes the activity of the j-th ion at the instant t, i.e. s jt . Matrix A contains the selectivity coefficients. The vector of valences is denoted by z = [z 1 , . . . , z n s ] T , and the operator ⊗ z describes the nonlinear transformation inside the logarithm function present in the NE model (see Eq. ( 1)). If the valences z i are equal, then ⊗ z results in a simple matrix multiplication. Finally, 1 1×n d corresponds to a vector (dimension n d ) with all elements equal to one.

In view of a possible model inaccuracy and/or of the errors introduced by the measurement system, a more realistic description for the array outputs is given by:

X = Y + N, (3)
where N ∈ R n c ×n d represents the noise terms. We assume a zero mean additive white Gaussian noise (AWGN) with covariance matrix C n = diag(σ2 ), where the elements of the vector σ 2 , the variances σ 2 1 , . . . , σ 2 nc , are unknown. With equations ( 2) and (3) in mind, we may formulate the source separation problem treated in this work: given the array response X and assuming that the valences z i are known, our goal is to estimate the elements of S (ionic activities). Since we envisage a blind method, the other parameters related to the mixing model (except z) and the noise variance at each electrode are also unknown and, thus, should be estimated. Furthermore, as it will become clear later, there are other unknown parameters, denoted by φ, which are related to the prior distributions assigned to the sources. Henceforth, all these unknown parameters will be represented by the vector θ = [S, A, d, e, σ, φ] and the following notation will be adopted: θ -θq represents the vector containing all elements of θ except θ q .

III. BAYESIAN SOURCE SEPARATION METHOD

The first step in Bayesian estimation is to assign an a priori probability distribution function, or simply prior, for each unknown parameter. This is done based on some available knowledge or constraint related to the unknown parameters; in Section III-A, we will discuss how this step can be accomplished in our problem. Later, in Section III-B, we will present the likelihood function p(X|θ) resulting from the mixing model (Eqs. ( 2) and ( 3)). Then, in Section III-C, we shall derive the posterior distribution, which can be seen as an update of the priors using the information brought by the likelihood function. According to the Bayes' rule, the posterior is given by p(θ|X) = p(X|θ)p(θ)/p(X). Given that p(X) is not a function of the unknown parameters θ, we are mainly interested in the following expression

p(θ|X) ∝ p(X|θ) p(θ). ( 4
)
Finally, in Section III-D, we shall discuss how an inference scheme based on (4) can be set to resolve our nonlinear source separation problem.

A. Defining the prior distributions 1) Prior distribution of the sources S: We adopted lognormal distributions to model the sources. This distribution is always non-negative which seems natural in our problem as the sources represent ionic activities. Moreover, there are two other aspects that motivate our choice. Firstly, the estimation of the parameters of a log-normal distribution is not complicated. As it will be checked later, it is possible to define conjugate priors in this case 2 . Also, there is a practical argument behind this choice. Ionic activities are expected to have a small variation in the logarithmic scale. This can be taken into account by the log-normal distribution, since such a distribution is nothing but a Gaussian distribution in the logarithmic scale. In mathematical terms, the prior assigned to s jt is given by

p(s jt ) ∝ 1 s jt exp - (ln(s jt ) -µ j ) 2 2σ 2 j 1 [0,+∞[ (s jt ), (5)
where µ j , σ j are the distribution parameters and 1 [0,+∞[ (s jt ) is the indicator function. We assume that the samples are independent and identically-distributed (i.i.d) and also that the sources are mutually statistically independent 3 , that is

p(S) = ns j=1 n d t=1 p(s jt ).
2) Prior distribution of the sources hyperparameters φ: The parameters φ j = [µ j σ j ] in Eq. ( 5) are also unknown in our problem and, thus, priors must be defined for them. We adopted a Gaussian prior distribution for µ j , that is:

p(µ j ) = 1 2πσ 2 µ j exp - µ j -µ µ j 2 2σ 2 µ j , ( 6
)
where µ µ j and σ µ j correspond to the hyperparameters. For r j = 1/σ 2 j , a Gamma distribution is considered:

p(r j ) = r αr j -1 j Γ(α r j )β αr j r j exp -r j β r j 1 [0,+∞[ (r j ), (7)
where α r j and β r j are the hyperparameters. In Section III-D, we will show that these two priors lead to conjugate pairs in the estimation of µ j and r j = 1/σ 2 j .

3) Prior distribution of the selectivity coefficients A:

According to the literature in potentiometric sensors (see [START_REF]Microcapteurs chimiques et biologiques[END_REF], [START_REF] Umezawa | Potentiometric selectivity coefficients of ion-selective electrodes[END_REF]), the selectivity coefficients are also non-negative. Moreover, it is rare to find a sensor whose response depends more on the interfering ion than on the target one, that is, a ij usually lies in the interval [0, A max ] where A max ∈ [0, 1]. Thus, an uniform distribution can be assumed for each 4 a ij , that is

p(a ij ) ∝ 1 [0,A max ] (a ij ). (8)
If no additional information is available, one can set A max = 1. However, it is possible to refine this information by using databases such as [START_REF] Umezawa | Potentiometric selectivity coefficients of ion-selective electrodes[END_REF] or by incorporating some information acquired during the fabrication process. Finally, the coefficients of A are assumed mutually independent, i.e.

p(A) = nc i=1 ns j=1 p(a ij ).
3 If some information concerning a possible dependency between the sources is available, it can be used for defining the priors. In this case, however, one goes toward a less general approach that is useful only to the modeled situation. Besides, the inference problem becomes more difficult in this case. 4 In view of Eq. ( 1), we consider a ij = 1 when i = j.

4) Prior distribution of the Nernstian slopes d: As already discussed, d i is related to physical parameters. For a room temperature, it takes approximately 0.059/z i , and the electrodes with such sensibility are said to have a Nernstian response. However, due to the sensor fabrication process and aging, a deviation from this theoretical value is usually observed. Furthermore, even the injection scheme can play a role on the value of this parameter [START_REF] Dorneanu | Computercontrolled system for ISEs automatic calibration[END_REF]. This possible deviation from the Nernstian value can be taken into account by setting a Gaussian prior of mean µ di = 0.059/z i V, i.e.:

p(d i ) ∝ exp - (d i -µ di ) 2 2σ 2 d i . ( 9
)
σ 2 d i must be high enough to correctly model the derivations from the theoretical value. For example, for an electrode whose target ion is monovalent, one obtains a coherent distribution by setting σ d i = 0.01. Again, this knowledge could be refined if additional information is available, or if the measurements are conducted in different temperatures. Finally, we assume that the elements of d are statistically independent.

5) Prior distribution of e: In contrast to the parameters d i , there is no theoretical value for e i . Still, as can be observed in [START_REF] Gründler | Chemical sensors: an introduction for scientists and engineers[END_REF], [START_REF] Hartnett | Potentometric nonlinear multivariate calibration with genetic algorithm and simplex optmization[END_REF], e i usually lies in the interval [0.05, 0.35]V. Hence, the following Gaussian prior may be adopted

p(e i ) ∝ exp - (e i -µ ei ) 2 2σ 2 e i , ( 10
)
where the mean is given by µ e i = 0.20V and σ 2 e i must be defined so [START_REF] Comon | Séparation de sources 1 : concepts de base et analyse en composantes indépendantes[END_REF] goes toward a flat prior in the interval [0.05, 0.35]. The elements of e are assumed mutually independent.

6) Prior distribution of the noise variance σ 2 i : A very common approach [START_REF] Moussaoui | Separation of non-negative mixture of non-negative sources using a bayesian approach and mcmc sampling[END_REF], [START_REF] Doucet | Monte Carlo methods for signal processing[END_REF] is to assign inverse Gamma priors for the noise variances; that is, γ i = 1/σ 2 i is modeled through a Gamma distribution with parameters α σ i and

β σ i p (γ i ) = γ α σ i -1 i Γ(α σ i )β ασ i σ i exp -γ i β σi 1 [0,+∞[ (γ i ) . ( 11
)
The motivation behind this modeling comes from the fact that it results in a conjugate pair. Moreover, it is possible to set α σi and β σi to obtain a non-informative prior [START_REF] Doucet | Monte Carlo methods for signal processing[END_REF].

B. Probabilistic modeling of the mixing process

Based on the mixing model defined in Section II and on the assumption of i.i.d. Gaussian noise which is also spatially uncorrelated, the likelihood is defined as

p(X|θ) = n d t=1 nc i=1 N x it e i + d i log ns j=1 a ij s z i /z j jt , σ 2 i , ( 12
)
where N x ik (µ, σ 2 ) denotes a Gaussian random variable in x ik with mean µ and variance σ 2 .

C. Bayesian inference

Having defined the priors p(θ) and the likelihood p(X|θ), we may now apply the Bayes' rule to obtain the posterior distribution p(θ|X). Given that the unknown variables of our problem are mutually independent (except S and φ), the prior distribution factorizes and, as a consequence, the posterior distribution becomes p(θ|X) ∝ p(X|θ) p(S|φ) p(φ) p(A) p(e) p(d) p(σ). [START_REF] Bedoya | Improving semiconductor-based chemical sensor arrays using advanced algorithms for blind source separation[END_REF] This a posteriori distribution can be used in the derivation of an inference scheme to estimate the parameters of interest. For instance, one may adopt the minimum mean square error (MMSE) estimator [START_REF] Kay | Fundamentals of statistical signal processing: estimation theory[END_REF], which is given by

θ M M SE = θ p(θ|X)dθ. ( 14
)
The exact implementation of this estimator is difficult because the analytical evaluation of this integral is hard in the problem considered here. Nonetheless, a good approximation of ( 14) can be provided by sampling methods. The idea here is to approximate the MMSE estimator using samples obtained from the posterior distribution p(θ|X). If, for instance, the generated samples are represented by θ (1) , θ (2) , . . . , θ (M ) , then expression ( 14) can be approximated by

θ M M SE = 1 M M m=1 θ (m) . ( 15
)
According to the law of large numbers, θ M M SE = θ M M SE as M → +∞. This important result gives the theoretical foundation for the above-described methodology, which is referred as Monte Carlo integration [START_REF] Robert | Monte Carlo Statistical Methods[END_REF].

D. Gibbs sampling scheme

The implementation of the MMSE estimator boils down to the task of generating samples from p(θ|X). In this work, this is done by the Gibbs' sampler, a Markov Chain Monte Carlo (MCMC) method tailored for simulating joint distributions. The idea in MCMC methods is to generate a Markov chain that admits the desired distribution (p(θ|X) in our case) as stationary distribution. Assuming that x ∼ p(x) stands for the sampling operation, i.e. x is a sample from the distribution p(x), then the Gibbs' sampler can be summarized as follows:

1) Set initial samples θ

(0) 1 , θ (0) 2 , . . . , θ (0) N ; 2) For m = 1 to M , do θ (m) 1 ∼ p(θ 1 |θ (m-1) 2 , θ (m-1) 3 , , . . . , θ (m-1) N , X) θ (m) 2 ∼ p(θ 2 |θ (m) 1 , θ (m-1) 3 , . . . , θ (m-1) N , X) . . . θ (m) N ∼ p(θ N |θ (m) 1 , θ (m) 2 , . . . , θ (m) N -1 , X) end.
Note that the Gibbs' sampler simulates a high-dimensional joint distribution by sequentially sampling from the conditional distribution of each variable 5 . Therefore, aiming at the implementation of the Gibbs sampler for our problem, we should obtain the conditional density for each element of θ. In the sequel, this will be accomplished by observing that p(θ q |θ -θ q , X) ∝ p(X|θ)p(θ q ) [START_REF] Moussaoui | Separation of non-negative mixture of non-negative sources using a bayesian approach and mcmc sampling[END_REF] and by considering the likelihood function [START_REF] Bermejo | ISFET source separation: foundations and techniques[END_REF] and the prior distributions (defined in Section III-A). 1) Conditional distributions for the sources: Substituting expressions ( 12) and ( 5) into ( 16), one has

p(s jt |θ -sjt , X) ∝ exp n c i=1 - 1 2σ 2 i x it -e i -d i log a ij s zi/zj jt + n s b=1,b =j a ib s zi/z b bt 2 - (ln(s jt ) -µ j ) 2 2σ 2 j 1 s jt 1 [0,+∞[ (s jt ). ( 17
)
2) Conditional distributions for the sources hyperparameters φ: Because p(X|θ) is not a function of the parameters φ j = [µ j σ j ], the conditional density of µ j is given by p µ j |σ j , S (j,:) ∝ p S (j,:) |µ j , σ j p(µ j ),

where S (j,:) denotes all the elements of the j-th row. By substituting ( 5) and ( 6) into [START_REF] Robert | The Bayesian Choice[END_REF], we obtain:

p µ j |σ j , S (j,:) ∝ n d t=1 exp - (ln(s jt ) -µ j ) 2 2σ 2 j exp - µ j -µ µ j 2 2σ 2 µ j . ( 19
)
The first term of ( 19) can be rewritten as a Gaussian function in µ j , with mean µ L µ j = 1/n d n d t=1 ln(s jt ) and variance σ 2 L µ j = σ 2 j /n d . Thus, [START_REF] Umezawa | Potentiometric selectivity coefficients of ion-selective electrodes[END_REF] becomes a product of two Gaussian distributions which is also a Gaussian distribution, i.e.:

p µ j |σ j , S (j,:) ∝ exp   - µ j -µ P ost µ j 2 2σ 2 P ostµ j    , ( 20
)
where 6 σ P ost µ j = σ 2

L µ j σ 2 µ µ j /(σ 2 L µ j + σ 2 µ µ j
) and µ P ost

µ j = (µ L µ j σ 2 µ µ j + µ µ µ j σ 2 L µ j )/(σ 2 L µ j + σ 2 µ µ j ).
Similarly to the case of µ j , the conditional density of r j = 1/σ 2 j is given by p r j |µ j , S (j,:) ∝ p S (j,:) |µ j , r j p(r j ). Thus, by considering ( 5) and ( 7), one has p r j |µ j , S (j,:) ∝

n d t=1 √ r j exp -0.5r j (ln(s jt ) -µ j ) 2 r α r j -1 j exp -r j β r j 1 [0,+∞[ (r j ). ( 21
)
We can rewrite p(r j |µ j , S (j,:) ) as a Gamma distribution, i.e. p r j |µ j , S (j,:) ∝ r

α P ostr j j exp - r j β P ost r j 1 [0,+∞[ (r j ), (22)
where α P ost r j = α r j + n d /2 and β -1

P ostr j = n d t=1 (ln(s jt ) -µ j ) 2 /2β r j
It is now clear that we obtain a conjugate pair and, thus, we can sample from p r j |µ j , S (j,:) by sampling from a Gamma distribution.

3) Conditional distribution of a ij : The derivation of the conditional distribution of p(a ij |θ -a ij , X) is close to the one conducted for s jt . Indeed, by considering expression ( 12), [START_REF] Moussaoui | Separation of non-negative mixture of non-negative sources using a bayesian approach and mcmc sampling[END_REF], and (8), it turns out that

p(a ij |θ -aij , X) ∝ exp - 1 2σ 2 i n d t=1
x it -e i -

d i log a ij s z i /z j jt + ns b=1,b =j a ib s zi/z b bt 2 1 [0,1] (a ij ). ( 23
)
4) Conditional distribution of d i : Using Eqs. ( 12) and ( 9), and after some straightforward calculations, p(d i |θ -d i ,X ) can be written as the following Gaussian distribution

p(d i |θ -d i , X) ∝ exp   - d i -µ P ost d i 2 2σ 2 P ost d i    , (24)
where σ P ost

d i = σ 2 L d i σ 2 µ d i /(σ 2 L d i + σ 2 µ d i ), µ P ost d i = (µ L d i σ 2 µ d i + µ µ d i σ 2 L d i )/(σ 2 L d i + σ 2 µ d i
), and

µ L d i = ( n d t=1 (x it -e i )) log n s b=1 a ib s z i /z b bt log ns b=1 a ib s z i /z b bt 2 , ( 25
)
σ 2 L d i = σ 2 µ d i log n s b=1 a ib s zi/z b bt 2 .
(26)

5) Conditional distribution of e i :

The development of this expression is similar to the one performed for d i . After some calculations, it is not difficult to show that

p(e i |θ -ei , X) ∝ exp - e i -µ P oste i 2 2σ 2 P ost e i , ( 27
)
where σ P ost e i = σ 2

L e i σ 2 µ e i /(σ 2 L e i + σ 2 µ e i
), µ P ost e i = (µ L e i σ 2

µ e i + µ µ e i σ 2 L e i )/(σ 2 L e i + σ 2 µ e i
), and

µ Le i = n d t=1 x it -d i log ns b=1 a ib s z i /z b bt n d (28) σ 2 L e i = σ e i /n d . ( 29
)
6) Conditional distribution of the noise variance σ i : As discussed before, the attribution of a Gamma prior for γ i = 1/σ 2 i culminates in a conjugate pair. Indeed, by considering expressions ( 12) and [START_REF] Jutten | Advances in blind source separation (BSS) and independent component analysis (ICA) for nonlinear mixtures[END_REF], one can show that p(γ i |θ -γ i , X) can be reduced to the following Gamma distribution:

p(γ i |θ -γi , X) ∝ γ α P ost σ i -1 i exp -γ i β P ost σ i 1 [0,+∞[ , ( 30
)
where α P ost σ i = α σ i + n d /2 and β -1

P ost σ i = 1 2 n d t=1 x it -e i -d i log ns b=1 a ib s z i /z b bt 2 + β -1 σi .

E. Algorithm Description

Let us make some remarks on the final algorithm, which is summarized in Tab. I. In a first step, we must define hyperparameters that lead to non-informative distributions [START_REF] Robert | The Bayesian Choice[END_REF]. For instance, in the experiments described in Section IV, this strategy was implemented by setting the following values for the variances of the Gaussian priors: σ 2 d i = 0.01 and σ 2 ei = 0.03. Moreover, a high value of the parameter σ 2 µj , which is related to the sources prior, was defined (σ 2 µ j = 100). Concerning the Gibbs' sampler, the resulting conditionals for almost all parameters are given by standard distributions (normal and Gamma) and, thus, sampling is straightforward in these cases. However, p(s jt |θ -s jt , X) and p(a ij |θ -a ij , X) result in non-standard distributions, thus requiring more sophisticated sampling methods. In this work, this task is accomplished by the Metropolis-Hasting (MH) algorithm [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]. Based on an instrumental distribution g(x; y), the MH algorithm simulates a given distribution p(x) in an iterative fashion; given x (t) , the current sample of p(x), the following iteration is done for obtaining the next sample x (t+1) :

1)

x * ∼ g(x; x (t) ) (generation of a candidate sample); 2) Calculate a = min 1, p(x * )g(x (t) ;x * ) p(x (t) )g(x * ;x (t) ) ; 3) u ← sample from a uniform distribution in [0, 1]; 4) if u ≤ a, x (t+1) = x * (accept the proposed sample)
else x (t+1) = x (t) (reject the proposed sample). In this work, truncated Gaussian are considered as instrumental distributions; in order to obtain good acceptance rate, we conducted, for each situation, preliminary simulations to adjust the variances of these distributions.

As the resulting Markov chain of the Gibbs' sampler takes some time to converge to its stationary distribution, the samples obtained in a first period (burn-in period) are not taken into account by the estimator (see Eq. ( 31)). A simple visual inspection of the evolution of some channel states was done for determining the burn-in period. However, more sophisticated strategies for convergence monitoring can be envisaged.

A last step in our algorithm concerns a post-processing stage for dealing with the scale and translation ambiguities inherent in blind methods [START_REF] Comon | Séparation de sources 1 : concepts de base et analyse en composantes indépendantes[END_REF]. For example, in PNL mixtures, the best we can do is to obtain an estimation sjt given by sjt = a j s jt + b j , where a j and b j are unknown, and s jt is the actual source. Given that this ambiguity in not acceptable in a sensing problem, we are forced to use at least two calibration points in order to retrieve the scale and translation parameters. This can be done by a simple linear regression, i.e. we need to find a j and b j that minimize the mean squared error 2 , where s (c) jn denotes the calibration points and s(c) jt the corresponding estimations.

J = 1/N cal N cal n=1 (s (c) jn -s (c) jn )

IV. EXPERIMENTAL RESULTS

To assess the proposed solution, we conducted in a first moment a set of experiments considering artificially generated ,α p j ,β p j ,α σ i ,β σ i ;

2) Random initialization of the current samples θ 0 ; 3) Run Gibbs sampler For m = 1 to M do 20)) σ m j ∼ p(σj |θ-σ j , X) (Eq. ( 22)) 24)) e m i ∼ p(ei|θ-e i , X) (Eq. ( 27)) σ m i ∼ p(σi|θ-σ i , X) (Eq. ( 30)) end 4) Infer the sources through the Bayesian MMSE estimation

• For j = 1, • • • , ns,t = 1, • • • , n d s m jt ∼ p(sjt|θ-s jt , X) (Eq. ( 17
)) • For j = 1, • • • , n s µ m j ∼ p(µj |θ-µ j , X) (Eq. (
• For i = 1, • • • , n c ,j = 1, • • • , n s a m ij ∼ p(aij |θ-a ij , X) (Eq. ( 23
)) • For i = 1, • • • , n c d m i ∼ p(d i |θ -d i , X) (Eq. (
sjt = 1 M -B M m=B+1 s m jt , ∀j, t, ( 31
)
where B denotes the number of iterations of the burn-in period. 5) Retrieve the source scale using N cal calibration points (N cal ≥ 2).

data. Then, in subsection IV-B, we test our proposal in a real situation. In both cases, the performance of each estimated source was inferred according to the following index:

SIR i = 10 log E{s 2 i } E{(s i -ŝi ) 2 } , (32)
where s i denotes the actual source and ŝi its respective estimation (after correct scaling). A global index is defined as SIR = 1/n s n s i=1 SIR i .

A. Experiments with artificial data

We made use of the database of selectivity coefficients presented in [START_REF] Umezawa | Potentiometric selectivity coefficients of ion-selective electrodes[END_REF] to define the following testing scenarios: T and e = [0.100 0.090] T . In the first scenario we have a PNL model as mixing system (see Eq. ( 1)). In the second one, the mixing system is composed of a nonlinear mixing mapping followed by componentwise logarithm functions. Finally, we consider noisy mixtures with a signal-to-noise ratio of SNR = 18 dB in both cases.

We tested our method in a situation where the sources are given by log-normal distributions (the number of samples was n d = 500). For the first scenario, we considered P = 50000 iterations for the Gibbs sampler with a burn-in period of B = 30000. The number of calibration points used in the postprocessing stage was N cal = 5. The performance indices in this situation were SIR 1 = 12.83 dB, SIR 2 = 19.0 dB, SIR 3 = Artificial data (first scenario): actual sources (gray) and their estimation (black). 17.4 dB, and SIR = 16.4 dB, which indicates that our proposal was able to achieve a good source separation. To illustrate that, we show in Fig. 1 the mixed signals, and in Fig. 2, the actual sources and their respective estimations.

In the second scenario, we considered P = 25000 iterations for the Gibbs sampler with a burn-in period of B = 18000, and N cal = 3 calibration points were used. The performance indices in this case were SIR 1 = 17.4 dB, SIR 2 = 16.2 dB, and SIR = 16.8 dB. Again, our method was able to provide fair estimations of the sources.

B. Experiments with real data

A well-known example of interference in ion sensing concerns the estimation of ammonium and potassium. These ions play an important role in applications such as water quality monitoring and food industry [START_REF]Microcapteurs chimiques et biologiques[END_REF]. A brief look at [START_REF] Umezawa | Potentiometric selectivity coefficients of ion-selective electrodes[END_REF] confirms that in many ammonium and potassium electrodes one has selectivity coefficients that are not negligible. Moreover, in some situations, these coefficients can reach quite high values, which represents a strong interference process.

We conducted a set of experiments where a solution containing K + and NH + 4 was analyzed through an ISE array composed of one K + -ISE and one NH + 4 -ISE. The activities of both ions were varied according to the following procedure. Initially, a solution containing a high concentration of potassium (10 -1 M) was prepared. Then, by injecting a solution containing ammonium chloride (NH + 4 Cl), the concentration of NH + 4 was increased whereas the concentration of K + was decreased. We repeat the same experiment but starting with a different initial concentration of K + (almost 10 -4 M). After joining these two experiments in a single dataset, we obtained the sources (the variation in time of ionic activities 7 ) shown in Fig. 5. The total number of samples of the two experiments altogether was n d = 170.

In Fig. 3, we present the responses of the ISE array. Since we have access to the inputs and to the outputs of the electrode array, it is possible to analyze the fitness of the NE model for this case, that is, we can have an idea about the amount of noise in the mixing model. Concerning the potassium electrode, we measured a signal-to-noise ratio (SNR) of SNR K + = 24 dB. For the ammonium electrode, this value was given by SNR N H + 4 = 20 dB. We assumed a Gaussian modeling for the noise in each electrode, and, to have an insight into the pertinence of such assumption, the distribution of the regression errors is plotted in Fig. 4. After scale and translation normalization with N cal = 4 calibration points, the performance indices (average of 30 experiments) obtained by the Bayesian method in this scenario were SIR 1 = 24.0 dB, SIR 2 = 22.5 and SIR = 23.2. 7 The ionic activities were estimated by using the actual concentrations and the Debye-Hückel formalism [START_REF] Gründler | Chemical sensors: an introduction for scientists and engineers[END_REF]. Concerning the parameters of the Gibbs sampler, the number of performed iterations was P = 10000 and the burn-in period was B = 7000. The obtained signals are shown in Fig. 5. Despite a small residual interference, mainly for the K + activity, the method was able to provide good estimations of the sources. On the other hand, the application of a PNL source separation method based on ICA [START_REF] Duarte | Blind source separation of post-nonlinear mixtures using evolutionary computation and order statistics[END_REF] provided poor approximations of the sources (SIR 1 = 7.6 dB, SIR 2 = -0.3 and SIR = 3.6).

C. On the number of calibration points used in the postprocessing stage

We argued in the introductory section that the main benefit brought by BSS methods is the reduction of the calibration step. However, as mentioned before, our BSS method needs at least two calibration points for retrieving the correct sources' scales. From this apparent paradox, we may ask ourselves why not simply use the available training points in a supervised context. For example, we could define a separating system (inverse of NE model) where the estimate ŝjt is given by ŝjt

= n c i=1 a * ij 10 x it -e * i d * i . ( 33
)
Then, we could estimate a * ji , d * i , and e * i based on a supervised approach, i.e. by minimizing the MSE error J =

n s j=1 n d t=1 (ŝ jt -s jt )
2 . We will refer to such strategy as the supervised-NE algorithm.

In fact, the argument presented in the last paragraph is somewhat naive since it implicitly assumes that the performance of a supervised method does not depend on the number of training samples, which is not the case. Indeed, in Fig. 6, we present the evolution of the training error of the supervised-NE method as the number of calibrations points N cal grows. In this figure, which represents an average of 100 runs, one can note that supervised-NE method requires, at least, N cal = 25 for providing good estimations. Conversely, our proposal requires only N cal = 5 to provide equivalent estimations 8 . This is achieved because the "real" data processing in our proposal is done by the Bayesian BSS algorithm, whereas the calibration points are only used to circumvent the signal ambiguities. Hence, although the complete suppression of the calibration stage is not possible, our method can work even if only a few calibration samples are available.

D. Discussion

The experiments with artificial data attested that the proposal can achieve good estimations even in the difficult case where the valences are different. In the scenario with real data, the Bayesian algorithm achieved a much better performance than the ICA-based PNL algorithm. In fact, in this case the sources were clearly dependent and, thus, they violated the fundamental assumption of any ICA method. Evidently, a scenario with dependent sources also poses a problem to a Bayesian method, since there is no guarantee that the provided data representation is unique in such a case. Nonetheless, in contrast to ICA, a Bayesian method does not optimize a functional associated with the statistical independence. Rather, it searches for a "good" data representation given a set of prior information. Thus, in the Bayesian approach, the independence should be seen as a simplifying assumption, that is, we are just omitting an additional prior information to the inference machine, which may still work.

Concerning the algorithm's convergence, we observed that the Gibbs sampler may get trapped in local minima, thus leading to poor estimations. For the situations with two sources the percentage of poor convergences9 was 3% for artificial data and 7% for the real case; in the case with three sources, this value attained 17%. Although we do not have access to the source, it is still possible to identify a bad solution: we observed that poor convergence usually implies in a significant mismatch between the representation provided by the Bayesian algorithm and the actual array response.

Another important point concerns the execution time of our method. For the experiments with real data, the MCMC algorithm took about10 26s to perform 10000 iterations. Nev-ertheless, the algorithm took 380s to perform 50000 iterations in the situation with three sources. This points out a wellknown drawback of MCMC-based Bayesian methods: their computational burden can become quite large as the number of sources and samples grows.

Finally, let us make a remark concerning the choice of the priors. In this work, we tried to define priors 1) that ease the resulting inference problem and 2) that, based on the available information, limit the range of the unknown parameters. Evidently, there is no guarantee that our choices are optimum and, as we mentioned before, a large dataset of measurements would permit to refine the priors' definition. An interesting issue in this context would be to compare the selected priors with alternative models. Because the possibilities are nonexhaustive, we consider a single example where, instead of log-normal priors for the sources, Gamma distributions are considered. This distribution also describes non-negative variables and, moreover, provides a flexible solution as it can model from sparse to almost uniform sources [START_REF] Moussaoui | Separation of non-negative mixture of non-negative sources using a bayesian approach and mcmc sampling[END_REF].

After applying the solution with the Gamma modeling to process the real data, the obtained performance (average of 30 experiments with N cal = 4 calibration points) -SIR 1 = 22.4 dB and SIR 2 = 17.5 dB -was inferior to that of the log-normal prior. Furthermore, unlike the log-normal case, it is not possible to find a conjugate prior in the estimation of the Gamma distribution parameters. Thus, it becomes necessary to incorporate an additional MH algorithm into the Gibbs' sampler, which increases the algorithm's complexity and has the inconvenient of requiring the definition of instrumental distributions, which is not an easy task.

V. CONCLUSIONS

In this work, we proposed a Bayesian source separation method for unsupervised quantitative analysis via an ISE array. A MCMC algorithm, the Gibbs sampler, was used in the inference stage. We defined the prior distributions based on informations that are usually available in a chemical sensing problem, such as the non-negativity of the sources. The resulting separation method did well in different scenarios, including in a real problem where the estimation of the activities of ammonium and potassium was desired. The results shown that ISE arrays equipped with blind source separation can operate even if only a reduced number of calibration points is available. This nice feature can pave the way for alternatives to supervised methods. For example, a less demanding calibration step could ease analysis in the field.

Despite the encouraging results, there are still some questions that demand future work. A first one concerns the possibility of increasing the precision of the developed source separation method. Indeed, in our experiments with real data, the remaining small interference would not be acceptable in very-high-precision applications. We believe that a more precise mixing model could increase the estimation quality. In this spirit, a first step is to search for alternatives to the NE model that take into account, for example, dynamics aspects of the interference modeling. Furthermore, it could be equally interesting to develop a Bayesian framework based on a temporal modeling for the sources, since chemical sources do have a time-structure.

•

  First scenario (n s = 3 and n c = 3): array of three electrodes (one K + -ISE, one NH + 4 -ISE and one Na + -ISE) to estimate K + , NH + 4 and Na + . Mixing system parameters: A = [1 0.16 0.40; 0.25 1 0.19; 0.40 0.13 1], d = [0.059 0.050 0.055] T and e = [0.095 0.105 0.110] T . • Second scenario (n s = 2 and n c = 2): array of two electrodes (one Ca 2+ -ISE and one Na + -ISE) to estimate the activities of Ca 2+ and Na + . Mixing system parameters: A = [1 0.39; 1 0.20], d = [0.026 0.046]

Fig. 1 .

 1 Fig. 1. Artificial data (first scenario): ISE array outputs (mixtures).

  Fig. 2.Artificial data (first scenario): actual sources (gray) and their estimation (black).

Fig. 3 .

 3 Fig. 3. Real data: responses provided by the ISE array.

Fig. 4 .

 4 Fig. 4. Histograms of the regression errors resulting from a fitting with the NE model. The black curves correspond to the fitted Gaussian distribution.

Fig. 5 .

 5 Fig. 5. Real data: Retrieved signals (black) and actual sources (gray).

Fig. 6 .

 6 Fig.6. Performance index as the number of calibration points grows.

  A Bayesian Nonlinear Source Separation Method for Smart Ion-selective Electrode Arrays Leonardo T. Duarte, Student Member, IEEE, Christian Jutten, Fellow, IEEE, and Saïd Moussaoui

In this work, log and ln stands for the logarithm of base 10 and the natural logarithm, respectively.

From the Bayes' rule p(X/Y ) ∝ p(Y /X)p(X). If p(X/Y ) and p(X) belong to the same family, then p(X) is said to be conjugate with respect to the likelihood p(Y /X). When the prior is conjugate, the Bayesian inference only needs the update of the posterior distribution (see[START_REF] Robert | The Bayesian Choice[END_REF], for instance).

Actually, the Gibbs' sampler requires only the conditional distributions up to a proportional gain.

The derivation of the mean and variance of a product of two Gaussian can be found in[START_REF] Kay | Fundamentals of statistical signal processing: estimation theory[END_REF].

In Fig.6, if the test error were plotted, then the difference between the supervised NE method and the Bayesian method would become even smaller.

These percentages were obtained through a visual inspection after 30 executions.

The method was implemented in Matlab (Windows XP) and the simulations were performed in a Intel Core 2 duo 3GHz, 2048MB RAM.

ACKNOWLEDGMENTS L. T. Duarte thanks the CNPq (Brazil) for funding his PhD research. The authors are grateful to Pierre Temple-Boyer, Jérôme Launay and Ahmed Benyahia (LAAS-CNRS) for the support in the acquisition of the dataset used in this work.

Finally, the authors also would like to thank the anonymous reviewers for their insightful comments.


G Barbut
email: barbut@imft.fr

M Braza

Y Hoarau

G Barakos

A Sévrain

J B Vos

Prediction of transonic buffet around a wing with flap

The present study presents numerical simulations and turbulence modelling of the flow around a NACA0012 airfoil including a deflected aileron. The results are compared with experiments that have been performed in the N-3 wind tunnel of the Institute of Aviation (IoA), Warsaw. The experiment focused on unsteady flow characteristics and buffet phenomena arising as the result of the transonic shock wave / boundary layer interaction (SWBLI). The transonic buffet is a natural and self-sustaining oscillation of the shock wave and separated flow region, caused by pressure fluctuation. The first objective is to capture the transonic buffet unsteadiness by means of URANS and DES turbulence modelling approaches. Secondly, the periodic flap oscillation has been used to modify the oscillation amplitudes towards an outlook of attenuation of the transonic buffet.

Introduction

Within the UFAST (Unsteady Effects in Shock Wave induced separation) program, it has been shown through its main test-cases, that URANS approaches were not successful to capture unsteady shock wave boundary layer interaction, SWBLI for internal geometries, whereas they were more successful for external flow configurations. These facts are reported in the UFAST books, published also by Springer in 2010. The DES approaches were more successful than URANS for internal flow configurations, whereas they are comparable to URANS for external flow configurations (example: the INCAS-biconvex test-case), where also transonic buffet occurs. In the present study, the transonic buffet over the NACA0012 wing with deflected aileron is investigated. First, the natural unsteadiness is examined as a self-sustaining oscillation of the shock wave and of the separated flow region, caused by pressure fluctuation (Seegmiller et al, 1978). A first objective of the present study is to capture the transonic buffet unsteadiness by means of URANS and DES turbulence modelling approaches. A second objective is to perform a sensitivity study concerning the effect of the periodic flap oscillation, that is used to modify the oscillation amplitudes. This goal is related to current and prospective studies for attenuation of the transonic buffet by deflected ailerons. In the present study, the URANS-Spalart-Allmaras (SA), the URANS-k-the DES-SA and the DDES -SA turbulence modelling approaches have been used.

Test-case Description 1 Experimental set-up

The NACA0012 airfoil model equipped with 22.6% of airfoil chord flap and modified gap as shown in Fig. 1 was tested. The model's chord is equal to 180 mm and flap length is 40.69 mm. The axis of rotation of the main part of the model is at 35% of the chord. There is a smaller, 0.3 mm gap between airfoil main part and flap leading edge. The flap leading edge radius is 4.69 mm. The axis of the aileron/ flap rotation is at 80% of the airfoil model chord. The span of the model is 596 mm. The computations have been carried by using the real wind-tunnel 3D geometry (fig. 2). The angle of attack is 5°. The Mach number is M = 0.73 and Reynolds number Re = 2.69•10 6 . The amplitude of the flap deflection angle up to   = ±2°.

The oscillation frequency of the flap, ff = 100 Hz and 30 Hz (close to the buffet frequency and one third of the buffet frequency, respectively).

Advanced signal processing on experimental signal

The experimental signals display a strong chaotic character, interspersed by peaks of more organised events in respect to the buffet phenomenon (figure 3). This complex behaviour was enhanced by the fact that the upper and lower walls are not porous but solid, fixed walls. IMFT has applied advanced signal processing techniques, as autoregressive modelling (AR), Burg algorithm as described in [START_REF] Marple | Digital Spectral Analysis[END_REF] and Morlet's wavelets analysis (Grossmann & Morlet-1984), to extract the organised part and to assess the contribution of the organised and random parts in the physical process. The AR model indicates organisation of the flow according to a bi-modal pattern with a predominant frequency of order 95-100 Hz, corresponding to the buffet phenomenon and a first sub-harmonic of order 50 Hz. The third sub-figure in fig. 3 shows the FFT analysis that provides a multitude of frequencies in the present case. A similar behaviour is obtained from the Morlet's wavelet analysis, according to the previously mentioned bi-modal pattern. In case of forced oscillation of the flap by a frequency f=35 Hz, the organisation of the signal shown by AR and by Morlet's wavelets indicates a shift of the organised motion around the forcing frequency (35Hz), as well as the appearance of a first harmonic of this frequency, of order 70 Hz (Doerffer et al, UFAST research program, 2010).

Numerical simulation of the IoA test-case

The transonic flow case IoA was studied for clean and flapped wing configurations. From the very start, the influence of the wind tunnel walls was identified as a key problem and a special mesh was put together to allow for modelling of the walls and easy comparison with free-stream conditions. For this Freq.

(Hz) Spectrum CT Burg, order=30, duration=0.05s, step = 0.0031 Duration=0.05s

Freq.

Freq.

Third Symposium on Hybrid RANS-LES Methods, 10-12 June 2009, Gdansk, Poland

The AR model indicates organisation of the flow according to a bi-modal pattern with a predominant frequency of order 95-100 Hz, corresponding to the buffet phenomenon and a first sub-harmonic of order 50 Hz. The lower sub-figure in fig. 3 shows the FFT analysis that provides a multitude of frequencies in the present case. A similar behaviour is obtained from the Morlet's wavelet analysis, according to the previously mentioned bi-modal pattern. In case of forced oscillation of the flap by a frequency f=35 Hz in the IoA experiment [START_REF] Miller | Parameters for aerofoil and aileron oscillations[END_REF], the organisation of the signal shown by AR and by Morlet's wavelets indicates a shift of the organised motion around the forcing frequency (35Hz), as well as the appearance of a first harmonic of this frequency, of order 70 Hz. The signal processing of this forced case is not presented in the present paper.

Numerical simulation of the IoA test-case

The transonic flow case of the IoA was studied for clean and flapped wing configurations. From the very start, the influence of the wind tunnel walls was identified as a key problem and a special mesh was put together to allow for modelling of the walls and easy comparison with free-stream conditions. At the time this report was compiled very few experimental data were available for comparisons and the exact flow case conditions were not specified. For this purpose, the University of Liverpool (ULIV) and IMFT focussed on a set of cases and comparisons with test data were made. The numerical grids are provided by ULIV (fig. 4 up) with 2.65 million cells including 18 blocks and by IMFT that performed refinement on the upper and lower walls as well as an optimised grid around the gap (fig. 4 down). The refinement in the upper and lower walls was done in respect of the y+ value (less than 1). IMFT final grids are of order 3.3 million and 3.75 million cells. A special attention has been attributed to correctly mesh the gap region and to provide CFD on a computational domain with the real wind-tunnel walls. The grid spacing in the mid-span is of order 0.03. The DES computation has been active in the flow regions beyond the near-wall ones.

The resolution of the wing surface was higher than the side walls which were put at exactly the same location as in the experiment. Care was taken to allow adequate grid resolution near the front 20% of the NACA-0012 wing section. This would allow for shocks to be resolved and capture any separation due to shock/boundary-layer interaction. The work focussed at relatively high incidence angles (up to 6.8 degrees) for which buffet clearly appears and for this reason, the upper surface of the wing was also meshed at a finer level. In total 2.65 million cells were used and the baseline grid included 18 blocks at each span-wise station. The time step size used was 10 -7 for the URANS and DES computations. Statistics were evaluated following the elimination of initial transient phases.. The statistics obtained for the DES would need to be checked for convergence by using longer sample sizes in a future study.

The numerical code used by ULIV is the PMB code [START_REF] Barakos | 48-month IoA Case Contribution[END_REF] and by IMFT the NSMB, (Navier-Stokes MultiBlock code) [START_REF] Vos | Recent advances in aerodynamics inside the NSMB (Navier-Stokes Multi-Block) consortium[END_REF] in which IMFT, being part of the NSMB consortium, has implemented higher-order turbulence models.

In the present study, the employed turbulence models concerning URANS are the Spalart Almaras (SA) (1992) and the k- model [START_REF] Wilcox | Turbulence Modelling for CFD[END_REF]. Concerning Detached Eddy Simulation, the DES-Spalart Allmaras (DES-SA) (1997), the DDES-Spalart-Allmaras (DDES-SA) ( 2006) have been used. This suggests that a zero gradient condition could be imposed in this direction and a slice of the wing near the centre would be taken into account for the simulation. This approach would result in lower computational cost. However, the top and lower walls have a stronger effect and as was investigated, their effect must be accounted for. The flow reversal near the mid-plane of the wing and the walls is clearly visible in Figure 5, as well as the corner vortices. The effect of the tunnel walls was studied. Computations with the upper and lower walls and with free-stream conditions result in different surface pressure coefficients Overall, the presence of the walls resulted in higher levels of surface pressure coefficient comparing with the non-confined case. Without the presence of the walls, the buffet onset is shifted to higher Mach numbers and angle-ofattack. This significant influence suggested that the upper and lower wind tunnel walls must be taken into account for the flow simulation.

A study of the spanwise conditions is presented on Figure 6. It can be seen that the full 3D configuration compares better with the experimental data than the computation where spanwise symmetry was imposed. It is therefore considered The buffet frequency obtained by the above IMFT results is of the order of 95 Hz, in good agreement with ULIV results and with the physical experiment. It is recalled that the experimental data processing in respect of the buffet frequency has been performed later than the CFD evaluation of this frequency. Therefore, the present test-case was a 'blind' one in respect of the predominant frequency value. The results presented in Fig. 10 Furthermore, comparisons of the CFD pressure coefficient with the experiment are presented on fig. 13. ULIV have used the DES-SA modelling, where IMFT has used the DDES-SA (Delayed DES). A quite good comparison of the mean wall pressure coefficient is achieved, especially by ULIV. The difference between the two DES could however be a consequence of the short preliminary statistical sample computed. Figure 14 shows URANS-SA and DDES computations at 50% of chord, for the suction and pressure sides. The URANS field has been used as initial condition for the DDES computations. The two signals are practically in phase opposition by 180°. The buffet frequency is clearly pointed out by both approaches. Longer computational times should be used for the DDES approach to capture a more rich statistical content.

Flap actuation for flow control

The effect of the flap actuation has been studied in two cases in the UFAST program, with actuation frequency of 30 Hz, close to one third of the buffet frequency and a frequency of 100 Hz, closed to the buffet frequency respectively. The flap was oscillating at an amplitude of 2 degrees. After an initial time where the flow adapts from a partially-converged steady-state to buffet, the loads show an almost periodic response. Figure 15 shows the lift and drag time-dependent coefficients with and without flap deflection. It is noticeable that the present actuation does not considerably attenuate the buffet frequency.

The results of Figure 16 present comparison of the controlled and uncontrolled cases with blended and slotted flaps. In the case of actuation with a frequency very close to the buffet natural frequency, it is found that the amplitudes of the oscillation increase. In the case of actuation with a third of the buffet frequency, these amplitudes decrease. Furthermore, the flap motion appears introducing extra frequencies in the flow and induce loads modulation. Both cases lead to an increase of lift coefficient amplitudes, as was pointed out by the computations of ULIV. This can be explained as a 'resonance' lock-in between the buffet frequency and the forcing frequency (100 Hz), created by the flap motion. Figure 18 shows the 3D modulation of the shock area along the span in the first case (flap oscillations starting at an instant of maximum Clift). Furthermore, this figure shows the influence of the phase concerning the starting of the flap motion, at minimum ("Phase 1") or at maximum ("Phase 2") of lift coefficient respectively. It is found that the buffet mechanism finally follows the phase of the flap motion oscillations.

Conclusions

A numerical study has been carried out to analyse the transonic interaction around a NACA0012 wing of constant section, including a trailing edge aileron (flap). First, the static case has been analysed and secondly, a forced flap motion case with frequencies equal to the natural buffet and to one third of the buffet frequency. The results are compared with the experimental results of the Institute of Aviation, IoA, provided in the UFAST research program.

A detailed signal processing of the experimental results has been carried out by IMFT by using especially wavelet analysis and Autoregressive (AR) modelling. This treatment was useful to extract the organised frequencies of the buffet phenomenon from the complex chaotic behaviour due to turbulence. A special attention has been attributed to examine the influence of a low amplitude (2°) flap motion.

The numerical simulations have been carried out respecting the exact wind tunnel geometry, by the University of Liverpool (ULIV) and the Institut de Mécanique des Fluides de Toulouse (IMFT). The approaches that have been used in the context of URANS and preliminary DES, show ability of capturing the buffet frequency, as well as the wall pressure coefficient, in a quite good agreement with the experiment. Comparisons were first carried out between the two sets of simulations without knowing the experimental values and afterwards, detailed comparisons with the experiment have been performed. The buffet frequency has been found around 95 Hz. Concerning the flap motion, both simulations have shown an increase of the amplitudes at forcing equal to the natural buffet frequency. A slight attenuation of the aerodynamic coefficients has been achieved in the case of forcing by one third. This study has shown that the phase response of the system is independent on the initial starting phase of flap motion. The present study opens a promising way towards a further significant control of the buffet, by using this kind of simulation approaches within an optimisation kernel. In this way, it may be possible to attenuate rather than amplify the oscillations arising from transonic buffet.

Fig. 1 .

 1 Fig. 1. The geometry of the NACA 0012 model equipped with the flap.

Fig. 2 .

 2 Fig. 2. Geometry of the test section, IoA N-3 wind tunnel with the NACA 0012 flapped model.

Fig. 3 .

 3 Fig. 3. Experimental signal, (up); autoregressive modeling (AR)-(middle), FFT spectrum (down).

Fig. 3 .

 3 Fig. 3. Experimental signal, surface pressure at 50% of chord versus time, (up); autoregressive modeling (AR)-(middle), FFT spectrum (down), M=0.7, =6.8, Re=2.63x10 6

Fig. 4 .

 4 Fig. 4. View of the baseline mesh constructed by ULIV with ICEM-CFD Hexa. Up: overall domain configuration and block boundaries; down: adapted mesh by IMFT around the gap.

Fig. 5 .

 5 Fig. 5. Overview of flow configuration: flow recirculation and corner vortices near the side walls of the tunnel. Streamlines and iso contours of pressure coefficicient (left). Iso-(U=0) velocity contours (right). IMFT: URANS-SA

  also the lateral walls in the simulation. An interesting observation is that the flow near the flap may be influenced by the exact geometry of the gap between the main element and the flap of the wing. The slotted flap configuration has a small effect near the trailing edge, fig.6(right). To enhance the possibility of using a simplified geometry for actuation of the flap, slotted and blended flap configurations were used. The conditions were Mach number of 0.75 and Reynolds number of 2.81 million. The angle of attach of the wing was changed between 0 and five degrees (where the onset of the buffet was close). At the same time the flap was deflected by three degrees. For all cases, blended and slotted flaps were used. Flow visualisation is shown in Fig.7for the flow region near the slot. The conditions were selected so that the pressure difference between the lower and upper wing surfaces drives the flow through the gap. As a result a small flow distortion is shown. This flow distortion is also evident in the surface pressure coefficient plots shown in Fig 6 (right).

Fig. 6 .

 6 Fig. 6. Cp wall pressure coefficient compared with the experiment; axisymmetric assumption providing lower Cp, versus full 3D WT geometry computation. Slotted mesh including the gap between main and flap elements (left). Comparison between blended and slotted flap configurations (right), Re= 2.81 x 10 6

Fig . 7 .Fig. 8 .Fig. 9 .

 789 Fig . 7. Iso-streamlines and iso-Mach, blended (on the left) and slotted (on the right) flap configurations. Flap deflected by 3°, URANS-k-ULIV).

  (left) show flow visualisation by means of instantaneous iso-U contours at the median plane. The SWBLI is clearly illustrated and compared with the experimental Schlieren visualisation (fig. 11 left) of the IoA. Figure 10 (right) shows rms of the pressure field at the suction side. This illustrates the region affected by the shock motion. Figure 11 (right) shows an instantaneous field of the ratio (grad/),  being the density. This compares quite well with the instantaneous Schlieren visualisation.

Fig. 10 .

 10 Fig.10. Left: Instantaneous Iso-U contours at the medium span. Right: rms fluctuation of midspan pressure illustrating the compressibility effects in the shock region contours; URANS-SA.

Fig. 11 .

 11 Fig.11. Left : IoA shadowgraph showing SWBLI on the NACA 0012 airfoil with undeflected flap at α = 5° , Mach number M = 0.70 and Reynolds number Re = 2.69•10 6 . Right: Instantaneous iso-contours of the ratio (grad/),  being the density.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Left: Time-averaged wall pressure coefficient at 30% spanwise position (CFD). Right: experimental wall pressure coefficient (right), URANS -SA model.

Fig. 14 .

 14 Fig. 14. Surface pressure signal at 50% of chord, URANS-SA and DDES-SA; up: suction side; down: pressure side.

Fig. 15 :

 15 Fig. 15: Flap oscillation of 2 degrees at the buffet frequency; left: lift coefficient and flap angle, f, versus time. Right: drag coefficient and f. M= 0.7, Re= 2.63 x10 6 (ULIV).

Fig. 16 .

 16 Fig. 16. Comparison of controlled and uncontrolled cases (flow parameters as in previous figure), ULIV, with two oscillation frequencies.

Fig. 17 .

 17 Fig. 17. Left: mid span iso-Mach lines at two different phases of the natural buffet oscillations (maximum and minimum of the Clift), corresponding to the starting of flap actuation. Right: Flap actuation at forcing frequency of 100 Hz and a deflection of ± 2°. In legend, "steady" means fixed flap. IMFT -URANS-SA.

Figure 17

 17 Figure 17 shows the influence of forcing at 100Hz (very close to the natural buffet frequency), actuated at two different phases (maximum and minimum of the Clift

Fig. 18 .

 18 Fig. 18. Instantaneous iso-Mach surface (M=1), illustrating the supersonic 3D area along the span (right), M=0.73, a=5°, R=2.63 x 10 6

Fig. 19 :

 19 Fig. 19: Wall pressure coefficient oscillations as a function of time at 75% chord position and at mid-span distance. "Upper surface" means suction side. Motion 1=starting flap motion at minimum of Clift oscillation. Up: Static reference case. Middle: comparison of

Acknowledgement

This work was sponsored by the UFAST research program, "Unsteady Effects in Shock Wave Induced separation", coordinated by IMP-PAN (P. Doerffer). The authors are grateful for the use of the French supercomputing facilities of CINES, IDRIS, GRID-5000.