Nassim Razaaly
email: nassim.razaaly@inria.fr

Giacomo Persico

Pietro Marco Congedo

Impact of Geometric, Operational, and Model Uncertainties on the Non-ideal Flow Through a Supersonic ORC Turbine Cascade

Keywords: Geometric Variability, Uncertainty Quantification, ORC Turbine, SU2, Non-Ideal Compressible Fluid Dynamics

Typical energy sources for Organic Rankine Cycle (ORC) power systems feature variable heat load and turbine inlet/outlet thermodynamic conditions. The use of organic compounds with heavy molecular weight introduces uncertainties in the fluid thermodynamic modeling. In addition, the peculiarities of organic fluids typically leads to supersonic turbine configurations featuring supersonic flows and shocks, which grow in relevance in the aforementioned off-design conditions; these features also depends strongly on the local blade shape, which can be influenced by the geometric tolerances of the blade manufacturing. This study presents an Uncertainty Quantification (UQ) analysis on a typical supersonic nozzle cascade for ORC applications, by considering a two-dimensional high-fidelity turbulent Computational Fluid Dynamic (CFD) model. Kriging-based techniques are used in order to take into account at a low computational cost, the combined effect of uncertainties associated to operating conditions, fluid parameters, and geometric tolerances. The geometric variability is described by a finite Karhunen-Loeve expansion representing a non-stationary Gaussian random field, entirely defined by a null mean and its autocorrelation function. Several results are illustrated about the ANOVA decomposition of several quantities of interest for different operating conditions, showing the importance of geometric uncertainties on the turbine performances.

Introduction

Organic Fluids are characterized by complex molecules and moderate to large molecular weights, resulting in critical point temperature significantly lower than that of water. By virtue of these features, the implementation of organic fluids in Rankine cycles, referred to as Organic Rankine Cycle (ORC) power systems, has prompted the development of widely distributed, thermal energy conversion devices suitable to exploit heat sources featuring low temperature level such as, for example, geothermal energy, concentrated solar energy, biomass or waste heat from industrial processes. These sources also share a certain degree of variability in the heat load, which might combine with variations in the condenser temperature due to the inherent annual variability of ambient temperature; as a result, the thermodynamic condition at the inlet and, in some cases, at the outlet of the ORC turbine undergoes significant variation, resulting in a variable expansion ratio and aerodynamic loading on each of the cascades of the machine (see [START_REF] Colonna | Real-gas effects in organic rankine cycle turbine nozzles[END_REF]). To improve the reliability of this technology, the robustness to variation in input conditions must be taken into account at an early stage of the development process.

On the other hand, just a few highly uncertain thermodynamic data are available for many ORC working fluids, especially the heaviest and most complex ones [START_REF] Cinnella | Quantification of thermodynamic uncertainties in real gas flows[END_REF][START_REF] Congedo | Backward uncertainty propagation method in flow problems: Application to the prediction of rarefaction shock waves[END_REF][START_REF] Cinnella | Sensitivity analysis of dense gas flow simulations to thermodynamic uncertainties[END_REF]. This complicates the development of high-fidelity Equations of State describing the thermodynamic behavior of these fluids, and leaves a certain margin for uncertainty in the parameters of the model equations (see [START_REF] Congedo | Numerical investigation of dense-gas effects in turbomachinery[END_REF] for a discussion). Furthermore, the need for a reduction in production costs of ORC systems may push ORC companies towards large-scale manufacturing techniques based on low-quality machining and assembling technologies, which might result in significant geometric tolerances on the final product [START_REF] Congedo | Shape optimization of an airfoil in a BZT flow with multiple-source uncertainties[END_REF].

On the basis of the aforementioned considerations, three areas of uncertainties can be identified when analyzing ORC turbines. Theoretical analysis of ORC turbines is nowadays performed by resorting to Computational Fluid Dynamics (CFD); in the classical approach, CFD simulations are run considering a pre-defined geometric configurations and by assigning a set of deterministic input data and model parameters. In order to evaluate the effects of geometry and parameter variation through a CFD model, several stochastic approaches have been developed in recent years [START_REF] Le Maitre | Stochastic Spectral Methods for Uncertainty Quantication with Applications to Computational Fluid Dynamics[END_REF] [8], which treat the problem solution as a vector of random variables depending, via the governing equations, on one or more stochastic input parameters. These may affect the solution via the equations themselves, via the initial and/or boundary conditions and via the geometry of the computational domain. The results of these stochastic simulations is no longer a single deterministic solution, but rather the probability distribution function of the solution or at least its lower-order statistical moments, such as the expectancy and the variance. By virtue of this class of simulations under uncertainty, the mean system performance can be evaluated along with its variability, and this information can be taken into account for system optimization. For instance, one can think of optimizing the mean system performance while minimizing some variability parameter to get a more robust design and avoid overfitting problems. The interested reader may refer, e.g., to [START_REF] Congedo | Shape optimization of an airfoil in a BZT flow with multiple-source uncertainties[END_REF] and [START_REF] Congedo | Tsi metamodels-based multi-objective robust optimization[END_REF] for a discussion of optimization strategies under uncertainties in real-gas scenario.

The present work, for the first time in the open Literature on ORC power systems, investigates the application of some advanced stochastic techniques to the analysis of a typical supersonic turbine cascade for ORC applications by considering a high-fidelity turbulent solver and including a detailed characterization of the geometric uncertainties. In fact, some recent works dealt with ANOVA-based analysis applied to the simulation of non-ideal gas flows, without an accurate characterization of geometric uncertainties: i) using Euler solvers and considering uncertainties on operating conditions and thermodynamics models [START_REF] Cinnella | Robust optimization of dense gas flows under uncertain operating conditions[END_REF][START_REF] Geraci | High-order statistics in global sensitivity analysis: Decomposition and model reduction[END_REF][START_REF] Congedo | Tsi metamodels-based multi-objective robust optimization[END_REF], ii) using turbulent solvers [START_REF] Bufi | Efficient uncertainty quantification of turbulent flows through supersonic orc nozzle blades[END_REF].

In particular, in this paper, a high-fidelity two-dimensional CFD model is used, coupled with the k -ω-SST model to introduce the effects of turbulence and with the Peng-Robinson-Stryjek-Vera Equation of State to treat the non-ideal thermodynamics of the organic vapour under consideration. Multiple sources of uncertainty are taken into account, related to the geometric tolerances of the blade and to the operating conditions, and also considering those of some parameters of the thermodynamic model. This paper is organized as follows. Section 2 provides a description of the CFD model used in this paper. How uncertainties are modeled is illustrated in Section 3. Then, Section 4 is devoted to the presentation of the UQ methods used for the sensitivity analysis. Section 5 illustrates the main outcomes of this paper and a physical analysis. Finally, some conclusions and perspectives are drawn in Section 6.

Computational Flow Model

In this section, the flow solver and the flow model used in this work are recalled. The turbine configuration and the operating conditions of interest are also introduced.

CFD Solver

The code SU2 [START_REF] Palacios | Stanford University Unstructured (SU2): An open-source integrated computational environment for multi-physics simulation and design[END_REF][START_REF] Economon | Performance optimizations for scalable implicit {RANS} calculations with {SU2}[END_REF] has been employed to carry out the numerical analysis presented in this work, in particular the solver version developed for non-ideal Computational Fluid-Dynamics (SU2 NI-CFD) [START_REF] Vitale | Extension of the SU2 Open Source CFD code to the simulation of turpinibulent flows of fluids modelled with complex thermophysical laws[END_REF]. The reliability of solver predictions was assessed in [START_REF] Gori | Experimental assessment of the open-source su2 cfd suite for orc applications[END_REF]. The SU2 NI-CFD solver relies on an embedded thermodynamic library which includes several Equation of States (EoS) and advanced turbulence models. In the following, a brief overview of the main features of SU2 are recalled. The three-dimensional Navier-Stokes equations for compressible flows read

∂u ∂t + ∇ • f(u) = ∇ • d(u) (1)
The vector of the conserved variables u = (ρ, m, E t ) T includes the mass density ρ ∈ R + , the three component vector related to momentum density m ∈ R 3 , and the specific total energy E t ∈ R, E t = ρ e + 1 2 ||v|| 2 (where e is the specific internal energy while v = m/ρ is the velocity vector).

The functions f

(u) = [m, (m ⊗ m)/ρ + PI, m(E t + P)/ρ] T and d(u) = [0, Π, v T • Π -q] T ∈ R 5 × R 3 are
, respectively, the advection and pressure fluxes and the viscous and thermal fluxes. In these latter functions, P = P(u) is the pressure, I is the identity matrix of dimension 3, Π = Π(v) is the viscous stress tensor and q is the thermal flux.

The system of equations is supplemented by constitutive relations that bound the state of the fluid to the thermodynamic and to the transport quantities. SU2 NI-CFD solves, in fact, the Reynolds-Averaged Navier-Stokes (RANS) equations, in which the effects of turbulence are introduced via the so-called turbulent or eddy viscosity and of the turbulent thermal conductivity, which are evaluated using well-established turbulence models, as reported in [START_REF] Economon | Performance optimizations for scalable implicit {RANS} calculations with {SU2}[END_REF].

The SU2 suite is equipped with Non-Reflecting Boundary Conditions (NRBC) [START_REF] Giles | Nonreflecting boundary conditions for euler equation calculations[END_REF], designed to avoid spurious pressure oscillations due to the reflection of acoustic waves at the domain boundaries and to perform computations on a truncated domain, therefore lowering the computational burden; this is particularly relevant in the present case, which features shock waves generated in the rear sucton side of the blades and propagating downstream up to the outflow boundary.

Flow Model and Turbine Configuration

In the analysis presented hereinafter, the inviscid fluxes are discretized using a second order approximate Riemann solver (ARS) of Roe upwind type [START_REF] Roe | Approximate riemann solvers, parameter vectors, and difference schemes[END_REF][START_REF] Vinokur | Generalized flux-vector splitting and roe average for an equilibrium real gas[END_REF][START_REF] Guardone | Roe linearization for the van der Waals gas[END_REF] along with the slope limiter proposed by van Albada; central differences are used for the viscous terms. The combination of such schemes guarantees high resolution in the numerical solution of the system of equations. The k -ω Shear Stress Transport (SST) model [START_REF] Menter | Zonal two equation kw turbulence models for aerodynamic flows[END_REF] is used to retrieve the effects of turbulence in RANS simulations. The organic vapour of interest is considered a polytropic non-ideal gas, whose volumetric behavior is described by the Peng-Robinson-Stryjek-Vera Equation of State. Uniform thermo-physical quantities were also assigned, estimated as representative values for the entire transformation by resorting to the RefProp library. The configuration of interest is a well known axial-flow supersonic nozzle cascade operating with siloxane MDM, first presented in [START_REF] Colonna | Real-gas effects in organic rankine cycle turbine nozzles[END_REF]. The properties of the MDM are reported in Table 1. In this work, blade-to-blade effects are investigated, considering a two-dimensional flow around the blade profiles at the midspan section of the cascade. Total Pressure P t in , total Temperature T t in , and axial flow direction are assigned at the inlet, while static pressure P s out is given at the outlet. A representation of the blade shape, alongside the computational mesh, is given in Figure 1.

Following [START_REF] Colonna | Real-gas effects in organic rankine cycle turbine nozzles[END_REF], we consider first a nominal (or full-load) operating condition for this nozzle cascade, characterized by an inlet thermodynamic state close to the saturation curve, and a high expansion ratio equal to ≈ 6. As well known, ORC power systems are requested to operate at part-load for long periods during their technical life, due to changes in the thermal power made available by the heat source and in the condenser temperature; this variation implies a large change in the turbine pressure ratio, resulting in a variation of aerodynamic loading on each cascade. In [START_REF] Colonna | Real-gas effects in organic rankine cycle turbine nozzles[END_REF], the implication of part-load operation for this cascades was estimated so to reduce the pressure ratio to ≈ 4, by an increase of cascade outlet pressure. Both full-load and part-load operating conditions are considered in this study, both of them reported in Table 2.2.

Grid Generation and Convergence

The computational meshes are generated by resorting to an advancing-front/Delaunay algorithm. To create hybrid grids suitable for viscous simulations, quadrilateral elements are first added over the solid walls, to build a boundary layer mesh. Afterwards, the advancing front algorithm triangulates the remaining portion of the computational domain. Since the cascade configuration features a typical fish-tail shock pattern at the blade trailing edge and the onset of a main shock in the unguided turning region on the rear suction side [START_REF] Colonna | Real-gas effects in organic rankine cycle turbine nozzles[END_REF], in these areas the grids were refined so to enhance the quality of the solution. A representative mesh composed by 180 kcells is reported in Figure 1, alongside two enlargements corresponding to the diverging channel and to the trailing edge. As no wall functions are used, the near-wall region of the boundary layer is solved entirely by the simulation. To ensure a proper resolution in the near-wall region, the height of the first cell of the boundary layer grid was set to equal to ≈ 2 • 10 -7 , corresponding to a y + value lower than 1. The y + value was evaluated considering the dynamic viscosity value assumed for the whole domain, and given in Table 1; the considered values of stream velocity and density are instead the ones resulting from simulations at the exit of the cascade, the reference length being the blade chord.

A dedicated grid dependence analysis was performed at the full-load condition. Two quantities were used as metric for evaluating the mesh convergence, namely the total pressure loss Y = P t in -P t out P t in -P s out (where P t out denotes the mass-averaged total pressure at the outlet), and the mass flow rate per unit span ṁ. Figure 2 reports the relative percentage error on the performance estimators. In this study, grid independence was considered achieved with the mesh composed by 590k cells; on the basis of the quantitative results, the 180k cells mesh (depicted in detail in Figure 1) was considered as the optimal compromise between accuracy and computational cost (error of 2%) and, hence, it was employed to compute the blade performance in the UQ framework. We perform a posteriori assessment that the variability induced by the uncertainty quantification study is mainly predominant with respect to the mesh error.

SU2 Solver/CPU

The RANS simulations were performed on a cluster equipped with Intel(R) Xeon(R) CPU X5650 at 2.67GHz. An adaptative CFL is selected to between 10 and 100. Simulations are assumed to be converged when residuals on the density are decreased by six orders of magnitude, or after 15000 iterations. Ill-converged simulations are re-run automatically with a minimal CFL of 1. Simulations run with 6 cores take about 8 hours for meshes composed by 180 kcells.

Cascade Flow Field at Full-and Part-Load Conditions

Calculations of the flow in the cascade for the full-and part-load conditions are now reported to properly show the most relevant flow features, in view of the UQ analysis which is the core of this study. Figure 3 reports the distribution of Mach number for the two conditions, and show that high supersonic flows (with peak Mach number values of 2) are established in the rear sides of the cascade and in the cascade-exit flow. As a result of the high Mach number of the free-streams (i.e., out from the boundary layer) on both the suction and pressure sides of the blade, classical fish-tail shock systems are generated at the blade trailing edge. The shock generated on the pressure side of the trailing edge impinges on the suction side of the adjacent blade, and it is reflected; this shock wave grows in strength when propagating downstream, as it merges with the compression occurring on the rear suction side of the blade, and due to the local concave shape of the blade. This leads to the onset of a strong shock wave in the cascade-exit flow field. A similar shock pattern can be observed for the two conditions, even though the main shock is stronger and slightly less inclined (with respect to the axial direction) in the part-load condition. The flow angle is also clearly different in the two conditions, as marked by the higher inclination of the wakes at part-load operation. This feature can be interpreted by considering that, in this latter condition, the converging-diverging nozzle becomes significantly overexpanded with respect to the outlet pressure; as a result, a post-compression takes place that, in cascade configuration, induces a flow turning in tangential direction (so to reduce the 'virtual' cross-section normal to the flow at the exit). The post-compression also explains the higher shock strength at part-load operation. The blade wakes and, especially, the shocks contribute to the loss generation, which is quantified in this study in terms of the aforementioned total pressure loss Y. Y and the mass flow rate per unit span ṁ are the Quantities of Interest (QoI) considered in that work.

Source and Modeling of Uncertainties

In this work, three classes of uncertainties are taken into account. At first, the geometric variability along the blade due to manufacturing tolerances was considered, discussed in subsection 3.1 and modeled through the realization of a random vector of practical dimension (N σ = 7 here), stemming from a infinite dimension gaussian process. Moreover, uncertainties in operating conditions as well as in two parameters of the thermodynamic model were also considered and are their quantification is discussed in section 3.2.

Modeling Geometric Variability

Various approaches have been introduced to model the geometric variations due to blade turbine manufacturing. Based on measurement data, Principal Component Analysis (PCA) can be used to build a probabilistic model of variability from the empirical mean and covariance of surface deviations at different locations on the blade [START_REF] Garzon | Probabilistic aerothermal design of compressor airfoils[END_REF][START_REF] Häcker | Statistical analysis of manufacturing deviations and classification methods for probabilistic aerothermal design of turbine blades[END_REF][START_REF] Lange | Principal component analysis on 3d scanned compressor blades for probabilistic cfd simulation[END_REF]. Following [START_REF] Dow | The implications of tolerance optimization on compressor blade design[END_REF][START_REF] Dow | Optimal design and tolerancing of compressor blades subject to manufacturing variability[END_REF], we assume that the geometric variability in manufactured turbine blades can be accurately described as a non-stationary Gaussian Random Field e(s, ω), ω being a coordinate in the sample space Ω, and (Ω, F , P) a complete probability space. The arclength s ∈ [0, 1] parametrizes the location on the blade surface, starting at the trailing edge (s = 0), going around the leading edge (s = 1

2 ), and continuing back to the trailing edge on the opposite side of the blade (s = 1). This approach is somehow similar to the PCA [START_REF] Garzon | Probabilistic aerothermal design of compressor airfoils[END_REF], except that in the present work, mean and covariance are provided in closed form rather than estimated from data.

Random Field. The Random field e(s, ω) represents the error between the manufactured surface and the nominal (perfect) one in the normal direction at the point parametrized by s. It is fully defined by its mean ē(s) (null here) and autocovariance function C(s, t): it captures the correlation between manufacturing errors at locations s and t along the blade surface, and describes the smoothness and correlation length of the random field. It is written as [START_REF] Dow | Optimal design and tolerancing of compressor blades subject to manufacturing variability[END_REF]:

C(s, t) = σ(s) σ(t) ρ(s, t) (2)
where σ(s) is the standard deviation of the random field at location s, considered constant here: σ(s) = σ 0 . σ 0 quantifies le level of manufacturing variability. The non-stationary autocorrelation function ρ is defined by [START_REF] Dow | Optimal design and tolerancing of compressor blades subject to manufacturing variability[END_REF]:

ρ(s, t) = exp(- |s -t| 2 L 2 (s, t) ) ( 3
)
where the square correlation length L 2 (s, t) is defined by:

L 2 (s, t) =L(s)L(t) L(s) =L 0 + (L LE -L 0 )exp(- |s -1 2 | 2 ω 2 ) (4)
The values L 0 = 0.1, ω = 0.1 [START_REF] Dow | Optimal design and tolerancing of compressor blades subject to manufacturing variability[END_REF] and L LE = 1.0 × 10 -2 , all normalized by the blade half-arclength were used.

In that study, trailing edge is modelled as a circular arc, thus, the impact of manufacturing variability at the trailing edge is not addressed. Discretization/Simulation. The Karhunen-Loeve (KL) expansion is used to represent the random field: it is based on a spectral decomposition of the autocovariance functions. The random field can be represented exactly by [START_REF] Betz | Numerical methods for the discretization of random fields by means of the karhunen-loève expansion[END_REF][START_REF] Karhunen | Über lineare Methoden in der Wahrscheinlichkeitsrechnung[END_REF][START_REF] Loeve | Functions aleatoires du second ordre[END_REF]]:

e(s, ω) = ē(s) + ∞ i=1 λ i Φ i (s)Z i (ω), (5)
where the λ i are arranged in descending order. Z i are standard uncorrelated normal random variables; λ i ∈ R + and Φ i : [0, 1] → R are respectively the eigenvalues and orthonormal eigenfunctions of the autocovariance function C, also referred to as the kernel function in this context, obtained by solving the homogeneous Fredholm integral equation of the second kind:

[0,1] C(s, t)Φ i (t)dt = λ i Φ i (s), ∀s ∈ [0, 1]. ( 6
)
This integral equation is solved numerically on the discretized blade surface using the Nyström method [START_REF] Betz | Numerical methods for the discretization of random fields by means of the karhunen-loève expansion[END_REF][START_REF] Nyström | Über die praktische auflösung von integralgleichungen mit anwendungen auf randwertaufgaben[END_REF], with N g Gauss quadrature points (N g = 500 here) to approximate the integral. Discretized Eigenmodes are consequently evaluated solving a discrete eigenproblem, using for instance a Singular Value Decomposition (SVD). It results in a finite expansion further truncated to give a stochastic dimension of N KL , determined by the decay of the eigenvalues. The normalized cumulative energy from the first p modes is defined by

E(p) = p i=1 λ i N g i=1 λ i . ( 7
)
It is used to select the number of modes retained in the KL expansion:

N KL = arg min{k ∈ N * s.t. E(k) < }, (8)
where is a user defined parameter driving the accuracy of the expansion. = 99% is chosen here. Figure 4 illustrates the eigenmodes decay and selection. By analyzing the effects of the modes, it is evident that the odd modes alter, with different 'wavenumber', the blade thickness: mode 1 appears as a simple omothetic scaling of the blade, mode 3 increases the thickness in the trailing edge region and makes the blade more slender in the leading edge, mode 5 makes the blade more slender in the central part of the blade and increases the thickness both at leading edge and trailing edge, and so on for mode 7. Conversely, the even modes act mainly as blade rotation: mode 2 appears indeed as a rigid blade rotation centered close to the leading edge, mode 4 and mode 6 are more complex deformations in which the leading edge, the trailing edge and the central part of the blade 'deflect' (by a different amount) with respect to the original blade.

Figure 5 shows samples of geometric perturbation with the proposed method. Three levels of geometric variability are considered to study its impact on the QOI:

• No geometric variability: σ 0 = 0.

• Moderate geometric variability: σ 0 = 3 × 10 -5 m.

• Large geometric variability: σ 0 = 6 × 10 -5 m.

In the following, the unity of σ 0 will be skipped for a sake of brevity.

Mesh Deformation

The introduction of geometric uncertainties involves issues related to the computational mesh for the flow simulation. It requires either to modify the initial mesh, or to resort to an automated mesh tool; the former approach was used here.

We rely on the use of a grid deformation in-house tool in order to accurately deform the grid for an assigned boundary displacement. The present approach follows the work of [START_REF] De Boer | Mesh deformation based on radial basis function interpolation[END_REF], successfully applied in [START_REF] Pini | Adjoint method for shape optimization in real-gas flow applications[END_REF], to achieve a highly flexible and robust deformation tool for unstructured meshes, based on the interpolation of boundary nodes displacements to the whole mesh with Radial Basis Function (RBF's). A linear system of equations only involving boundary nodes has to be solved, and no grid connectivity information is needed.

In this subsection, the following notations are. d denotes the dimension space (here d = 2), x ∈ R d a node in the mesh, x new ∈ R d its new location, n b the number of boundary nodes, {x b j } j∈ 1,n b the set of nodes at the boundary, p a polynomial.

The interpolation function s : R d → R d , describing the displacement in the whole domain, can be written as a sum of basis functions:

s(x) = n b j=1 α j Φ( x -x b j ) + p(x), (9)
where Φ is a given RBF. Coefficients of the linear polynomial p and α j are determined by the interpolation conditions:

s(x b j ) = d b j (10)
n b j=1 α j q(x b j ) = 0 ( 11
)
d b j
is the imposed displacement of the boundary node x b j . Equation 11 has to satisfied for all polynomials q with a degree less or equal than that of polynomial p. In our case, the displacement d b j is nullified for all boundary nodes, except the ones belonging to the blade. Independently for each spatial direction, the coefficients of the polynomial p and α j are recovered solving a (n b + 4) × (n b + 4) symmetric positive definite linear system, using a Cholesky decomposition for instance.

The new position x new of a node in the interior domain initially located in x is then directly derived:

x new = x + s(x) (12)
Each point is moved individually involving that no mesh-connectivity information is needed. The RBF retained in this work is the so-called Volume function, with global support, defined by Φ(r) = r.

Operating Conditions and Thermodynamic Model Uncertainties

Uncertainties associated to the thermodynamic conditions are considered in this work, namely on P t in , T t in and P s out , for both full-load and part-load turbine operation. The values of the uncertainties are reported in Tables 3 and4. The uncertainties were estimated considering typical measurement errors in the monitoring devices of power systems; in particular, ±5kPa of uncertainty in the pressure measurements and ±0.35K of uncertainty in temperature measurements are assumed.

Two coefficients associated to the thermodynamic model of the fluid are also considered uncertain, namely the acentric factor ω and the heat capacity ratio γ of the fluid (which is considered polytropic in this study). A specific quantification of the variability of these parameters not being available, a uniform uncertainty of 5% is assumed for ω and a uniform range ±0.004 ensuring thermodynamical consistency for γ. This is consistent with respect to other studies performed in literature over the coefficients of the thermodynamic models (see for example [START_REF] Cinnella | Sensitivity analysis of dense gas flow simulations to thermodynamic uncertainties[END_REF][START_REF] Cinnella | Quantification of thermodynamic uncertainties in real gas flows[END_REF][START_REF] Congedo | Tsi metamodels-based multi-objective robust optimization[END_REF]).

Tables 3 and4 summarizes the scenarios corresponding respectively to full-and part-load, each considering three different geometric uncertainty levels.

Globally, the random vector is then denoted as θ = (P

UQ methods

This section introduces the techniques employed for the quantification of the uncertainties on the quantities of interest.

UQ Methods for scalar QoI

Let us introduce the generic model Y(θ) = g(θ), where g is a scalar response function representing the QoI associated to a physical model depending of a vector θ. This is a set of d independent random variables, where each θ i features a Probability Density Function (PDF) f θ i . In the following, we provide a basic description about the methods for evaluating the statistics, and the Sobol Indices (ANOVA).

Note that in practice, the generic model Y is replaced by a Kriging-based surrogate model [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF][START_REF] Dubourg | Adaptive surrogate models for reliability analysis and reliability-based design optimization[END_REF] Ỹ. An anisotropic Matern 5 2 kernel is chosen for the covariance function, along with a null trend. The hyperparameters are tuned maximizing the likelihood. The latter is built based on Latin Hypercube Sampling (LHS) evaluations.

UQ Surrogated-based: Scalar Statistics Evaluation

In this study, the following statistics of Y(θ) are considered: mean

µ Y = E θ [Y], standard deviation σ Y = V 1 2
θ (Y) and α-quantile = q Y α . We recall below their definition. The expectation of f (θ), a.k.a. its mean, is defined as:

E θ [Y] = Y(θ) f θ (θ)dθ. ( 13
)
where f θ designs the PDF of θ, defined as d i=1 f θ i here, due to the independence assumption made above. The so-called variance of f (θ) is defined as:

V θ [Y] = E θ [Y 2 ] -E θ [Y] 2 . ( 14
)
For α ∈ [0, 1], the α-quantile, is defined by:

q Y α = inf{q ∈ R s.t. P θ (Y < q) α}. (15)
Otherwise speaking, it measures the probability that Y, due to the variability of θ, is lower than q Y α is α, i.e. P θ (Y < q Y α ) = α. The Coefficient Of Variation (CoV) defined by CoV = σ Y µ Y is also used to quantify the normalized variability of

Y.

Corresponding empirical estimators based on Monte-Carlo technique are computed with the following algorithm:

• Build a DOE (LHS) and Compute a CFD computation for each sample, yielding: {θ j , Y(θ j )} j∈ 1,N LHS

• Build a surrogate from the DOE: Ỹ(θ)

• Sample a large MC set:

{θ k } k∈ 1,N MC
• Compute MC empirical statistics μ Ỹ , σ2 Ỹ and q Ỹ α based on the surrogate Ỹ using Equations 16, considering the set { Ỹ(θ k )} k sorted:

μ Ỹ = 1 N MC N MC k=1 Ỹ(θ k ) σ2 Ỹ = 1 N MC N MC k=1 ( Ỹ(θ k ) -μ Ỹ ) 2 q Ỹ α = Ỹ(θ αN MC ) (16)

ANOVA: Sobol Indices

Variance-cased sensitivity indices are common tools in the analysis of complex physical phenomena. More precisely, we are interested in the total Sensitivity Index (SI) [START_REF] Saltelli | Sensitivity Analysis in Practice[END_REF][START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates[END_REF].

The variance decomposition of the response can be written as follows:

V(Y) = u⊆U σ 2 u (θ u ), (17)
where U = (1, 2, ..., d) is the set of random variables indexes and σ 2 u is the variance introduced by interactions of random variables θ u ⊆ θ U . The associated sensitivity measure of θ u , named Sobol SI, is written as the correlation ratio:

S u = V(E[Y|θ u ]) V(Y) = v u σ 2 v V(Y) . ( 18
)
For a single variable θ i , the first order Sobol SI S i is given by:

S i = V(E[Y|θ i ]) V(Y) . ( 19
)
It quantifies the contribution to V(Y) of the main effect of θ i varying alone, but averaged over variations in other input parameters. Note that d i S i = 1. The total order Sobol SI S T i is defined by:

S T i = E[V(Y|θ (i) )] V(Y) = 1 - V(E[Y|θ (i) ]) V(Y) , ( 20
)
where the θ (i) notation indicates the set of all variables except θ i . It measures the contribution to V(Y) of θ i caused by its interactions with all other input variables, with d i S T i 1. The Sobol SI are estimated using a Monte-Carlo algorithm [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates[END_REF].

UQ methods for High Dimensional QoI

We are interested in evaluating efficiently statistics (mean and variance) where the quantity of interest is a vector, possibly with a large dimension (say 10 5 ), for example to determine the distribution of statistical quantities over the computational domain. As a consequence, the method presented in Subsubsection 4.1 is not tractable, since it would require to build a large number of surrogate models. Let us consider the generic random variable G(X) where G : R d → R n is a multidimensional function depending on a random vector X characterized by its PDF f X . In the following, we describe a stochastic method combining a PCA (Principal Component Analysis) decomposition to reduce the dimension space n, surrogate modeling and MC sampling. Below, a summary of the method:

1. Set an initial DOE: Data generation using LHS 2. PCA: compute the eigenmodes, and select the most energetic ones. 3. Build a surrogate of each scalar coefficient in the reduced basis. 4. MC sampling on the scalar coefficients. 5. Recover the statistics of G 1) Initial DOE. N LHS samples using LHS are generated, then stored in the mean subtracted data matrix Y = [G(x i ) ( j)µ ( j) G ] i j ∈ R N LHS ×n . G(x i ) ∈ R n denotes the i-th sample. µ G ∈ R n denotes the data empirical mean, defined as ( j) . Mean subtraction is an integral part of the solution towards finding a principal component basis minimizing the data approximation mean square error [START_REF] Miranda | New routes from minimal approximation error to principal components[END_REF].

µ ( j) G = 1 N LHS N LHS k=1 G(x i )
2) PCA and modes selection. The semi-definite positive empirical covariance matrix is built:

C = Y T Y N LHS -1 (21) C ∈ R n×n is then decomposed into an orthogonal basis of eigenvectors V = (v 1 , ..., v n ), involving v T i v j = δ i j , as C = VΛV T (22)
Λ = diag(λ 1 , ..., λ n ) is a diagonal matrix where the positive eigenvalues are sorted in decreasing order. Similarly to Equation 8, the l most energetic modes are retained, with = 99.99%.

Note that in the cases that we consider, N LHS < n, so the covariance matrix C has a rank N LHS at maximum so l ≤ N LHS . The l eigenvectors {v k } k∈ 1,l with v k ∈ R n are set as basis vectors, the vector solution G(x) at an unknown location x is sought in the form

G(x) -µ G = l i=1 h i (x)v i ( 23
)
where h i (x) are the scalar coefficients in the reduced basis, which need to be approximated by a surrogate.

3) Build a surrogate.

For each i ∈ 1, l , the DOE {x k , v T i (G(x k ) -µ G )} k∈ 1,N LHS is used to build the surrogate hi for h i . A surrogate for G consequently reads G(x) = l i=1 hi (x)v i + µ G ( 24
)
4-5) MC sampling and statistics of G. The decomposition in the reduced basis (Equation 23) is used to approximate the random vector. From the linearity of the operator E X , we obtain:

E X [G(X)] = µ G + l i=1 µ h i v i ( 25
)
where

µ h i = E X [h i (X)
] denotes the mean of the coefficient h i . Please note that this quantities differs from the mean µ G only for the inherent error of using a limited number of modes, hence it is of the order of 1 -. The variance of G(X) is derived by component j ∈ 1, n , since they are uncorrelated w.r.t. the random vector X:

V X [G(X) ( j) ] = V X [ l i=1 h i (X)v ( j) i ] = l i=1 v ( j) i v ( j) i σ 2 h i + 2 n i=1 n k=i+1 v ( j) i v ( j) k ρ ik ( 26
)
where

σ 2 h i = V X [h i (X)] and ρ ik = E X [h i (X)h k (X)] -µ h i µ h k
denote respectively the variance of h i and the covariance between h i and h k . This derivation (Equation 26) takes into account that the coefficients h i can be correlated.

µ h i , σ 2 h i , ρ ik are evaluated numerically by means of MC on the metamodels hi so the mean (resp. variance) of G(X) are easily reconstructed using Equation 25(resp. Equation 26).

Computational aspects. Note that in practice, the eigen-decomposition (Λ, V) is not obtained by evaluating the covariance matrix C and then perform the diagonalization as suggested by Equation 22. Indeed, C ∈ R n×n can be very large and ill-conditioned since cond(C) = cond(Y) 2 . A Singular Value Decomposition (SVD) on the mean subtracted data matrix Y is preferred, yielding

Y = UΣV T (27)
U ∈ R N LHS ×N LHS is an orthonormal matrix with the property

U T U = I N LHS . Σ ∈ R N LHS ×n is a diagonal matrix of eigenvalues {σ k } k . Equation 27 is justified by C = (UΣV T ) T UΣV T N LHS -1 = V Σ T Σ N LHS -1 V T (28)
Equation 28 also shows the underlying relation between eigenvalues of matrices C and Y:

λ i = σ 2 i N LHS -1 (29)

Stochastic Convergence

An analysis of the convergence of statistics is performed considering full load and σ 0 = 6 × 10 -5 in order to justify the number of samples necessary to conduct the UQ analysis based on the tools presented above. Three sets of LHS DOE of respective size 100, 200 and 400 are considered. Those sizes are chosen on a basis of similar studies conducted with a inviscid model, not presented here for a sake of brevety.

Based on the three LHS sets, Figure 6 plots the PDFs of ṁ and Y, while Figure 7 summarizes their total Sobol SI of input parameters. PDF's and total Sobol SI are in very good agreement for the different LHS sets; however, a slight error in the PDF of Y is observed when using 100 samples, hence all the UQ studies were conducted based on a LHS set of size N LHS = 200.

Results

This section presents a comprehensive uncertainty quantification on the basis of the stochastic framework formulated above. Based on a LHS set of size N LHS = 200 for each case, a comprehensive UQ analysis is performed. In Subsection 5.1, a sensitivity analysis (ANOVA) is carried out. In Subsection 5.2, marginal and joint PDF's of the QoI are studied, their statistics in Subsection 5.3. In Subsection 5.4, the impact of the uncertain input on the flow field is studied.

Sensitivity Analysis

A Sensitivity analysis was performed for the QOI Y and ṁ, and is presented in Figures 8 and9 respectively, for both full-and part-load. First considering the influence of thermodynamic conditions in absence of geometric variability, the uncertainty in the outlet static pressure has the highest impact on the loss coefficient, especially in fullload condition; this is because the uncertainty in the outlet static pressure is relatively much higher than that on the inlet total pressure, and hence it dominates on the variability of pressure ratio; the variability of pressure ratio, in turns, induces a variation of shock strength and, hence, aerodynamic losses. In part-load condition, the inlet total pressure grows in relevance on Y; as a matter of fact, for this off-design condition the difference between the inlet and the outlet pressures reduces, while the uncertainty remains the same, so both the uncertainties contribute in a measurable way to the variability of pressure ratio. The outlet static pressure does not have any impact on the flow rate, as expected since the cascade is chocked in both full-and part-load conditions. Conversely, the uncertainty in inlet total pressure dominates the variability of flow rate in both the conditions, as the inlet total pressure determines the value of fluid density; this latter, in turn, is only marginally influenced by the inlet total temperature in the conditions of interest.

The uncertainties in the two parameters of the thermodynamic model have a minimal impact on the variability of Y. Conversely, the ratio of specific heat capacities has a relevant impact on the flow rate, as it influences the relationship between the mass flow rate and the expansion ratio in chocked-flow conditions. The acentric factor does not have a measurable effects on both the QoIs. The relevance of operational and model uncertainties is drastically reduced when manufacturing tolerances are plugged into the flow model, even though the mutual hierarchy of the formers remains the same, and the considerations reported above still apply. If geometric variability is large, manufacturing tolerances dominate the variability of the QoIs. In particular, it is observed that mode 1 and mode 4 are responsible for most of the variability of flow rate, and also have the highest impact on the variability of loss coefficient. An analysis of the geometric modes reveals that Mode 1 induces a net variation of geometric throat and also mode 4 has a large impact on the throat, as it is associated to a rotation which alters significantly both the leading edge and the rear side of the blade, which define the geometric throat of the cascade. These considerations justify the relevance of mode 1 and 4 on the flow rate; however, in a supersonic configuration a variation of the geometric throat without a corresponding change in expansion ratio makes the nozzle to operate in off-design conditions, leading to the onset of stronger shocks and, hence, modifying the cascade performance. This may explain why mode 1 and mode 4 also dominate the variability of Y. This latter, however, is influenced by all the modes; as a matter of fact, the development of boundary layer and the onset of the leading shock depend on the detailed shape of the blade, both in the diverging part of the bladed duct and in the region of un-guided turning, and all the modes alter the blade in the rear part. As a further consideration, it is well known in cascade aerodynamics that trailing edge thickness has a relevant impact on the wake width and, hence, on the viscous loss; this contributes to explain the relative relevance of all the odd modes on Y, as all of them alter the trailing edge thickness; this is particularly visible in full-load condition, for that viscous effects contribute more than shocks to the aerodynamic loss, as at full-load the shocks have lower strength.

PDF Comparisons

The PDFs of Y and ṁ are provided in Figure 10 for the three different geometric variations, in both full-load and part-load cases. In an alternative two-dimensional representation, the joint PDFs are also provided in Figure 11.

In absence of geometric variability, the PDFs of Y have a complex non-gaussian shape with different mean values and support; in particular, the support is narrow enough to make the full-load conditions always outperforming the part-load one. Conversely, the PFDs of ṁ are of gaussian type with very narrow support; this result had been already observed and discussed in detail in [START_REF] Romei | Assessment of deterministic shape optimizations within a stochastic framework for supersonic organic rankine cycle nozzle cascades[END_REF], and is motivated by the combination of inlet pressure and temperature uncertainties. As the cascade is chocked in the two conditions, and the uncertainties in the inlet conditions are the same, the PDFs appear identical. The influence of geometric tolerances is, again, very significant from both the qualitative and quantitative point of view. At first, the introduction of geometric variability makes the PDFs of Y to become of gaussian type. This is probably because the geometric variability itself is modeled as gaussian process; however, this indicate the dominant effect of manufacturing tolerances, when they are plugged into the uncertainty framework. On the quantitative ground, the geometric variability clearly enlarges the support of the PDFs in all the conditions and for both the QoIs. In particular, the variability in loss coefficient increases from about 2% in absence of tolerances to about 4-6% for the largest tolerance. Despite this, the mean of the PDFs undergoes a slight change (only visible for part-load condition) which does not alter the ranking between full-load and part-load conditions in terms of mean performances; however, when considering the actual realizations of the flow process in presence of the largest geometric variability, a significant overlapping between the PDFs in the two conditions exists (between 21.5% and 24.5%) making the two operating conditions equivalent from the performance perspective for a relatively large percentage of realizations.

The impact of manufacturing tolerances on the flow rate is identical for the two conditions and simply enlarges by a factor ≈ 5 the support of the PDF.

Statistics

On the basis of the PDFs reported in the previous section, the most relevant statistics (mean, variance, 5-quantile, and 95-quantile) were computed and are reported in Table 5. They show, in more quantitative terms, the features already observed when commenting the PDFs in terms of mean, variance, and support (which is roughly given by the difference between the two quantiles). The interference between the PDFs of Y at large geometric variability is confirmed by the 5-quantiles and 95-quantiles, thus confirming that it is not just a matter of tail probabilities. The analysis also confirms that the quantitative equivalence between the flow rate in the full-and part-load conditions to be gaussian processes whose mean correspond to the original blade shape. However, the same uncertainties produce significant perturbations, that grow significantly when activating the geometric variability. CoV values slightly exceeding 2% are found, with the highest values concentrated in the wakes and in the shock fronts.

In case of no geometric variability, almost null CoV is found in the bladed channel up to the shock on rear suction side of the blade; after that, high CoV is observed in the whole rear suction side up to the trailing edge, where local peaks of CoV occur in the wake and in the fishtail shock system. Moving downstream, the wake width increases due to turbulent and viscous diffusion, resulting in a cascade of wide regions of moderate CoV; in correspondence to the main shock front, the CoV remains locally high and increases moving downstream, as the leading shock grows in strength. As already commented, the main shock grows in magnitude moving downstream as the compression waves generated in the rear suction side, and induced by the local concave shape of the blade, coalesce in a shock front at a certain distance from the blade, enforcing the weak reflected shock coming from the trailing edge of the adjacent blade. The variability in operating conditions (in outlet pressure, in particular) and in the specific heat capacity ratio influence the shock strength justifies the local high values of CoV.

When geometric variability is plugged into the uncertainty framework, the character of the CoV distribution does not change but some evident quantitative effects appear; similar features take place for the two operating conditions. At first, progressively higher CoV values are found upstream of the blade when the manufacturing tolerance increases; the CoV reduces in the accelerating region of the bladed nozzle, and especially in the supersonic divergent, where the mean Mach number grows significantly. This suggests that manufacturing tolerances induce a nearly uniform variance in this region, so that the CoV distribution is dominated by that of the mean. On the rear suction side of the blade and in the cascade-exit flow field high CoV values are observed in the wakes and in the shock regions, as already observed in absence of geometric variability; however, wider areas of high CoV appear, and, at a certain distance from the blade, moderate-to-high CoV cover the entire flow field. The larger extension of high CoV regions downstream of the cascade can be explained by considering that the modes alter significantly the blade shape in the rear part. In particular, the compression wave generated on the concave part of the rear suction side is perturbed by the local variability in the blade shape, altering both the shock strength and the shock inclination; this contributes to explain the wider high CoV region observed across the shock in presence of manufacturing tolerances. Beside that, the wake appears even more sensitive than the shock to geometric uncertainty; this is probably motivated by the fact that the geometric variability, and in particular the odd modes, alter the trailing edge thickness, which has a direct impact on the wake width.

Conclusion

This paper has proposed an uncertainty quantification analysis applied to the simulation of a turbulent non-ideal flow within a supersonic turbine cascade for Organic Rankine Cycle applications. In particular, a Kriging-based method has been coupled with a Computational Fluid Dynamic solver permitting to consider multiple sources of uncertainties associated with operating conditions, fluid parameters and an accurate representation of geometric tolerances. Two main scenarios are considered, i.e. a full-load (a high expansion ratio equal to ≈ 6) and a part-load (pressure ratio to ≈ 4) operating conditions.

The primary outcome of the paper is an accurate statistical study of two quantities of interest, namely the total pressure loss and the mass flow rate. In particular, the sensitivity analysis, the probability density functions, and the statistical moments are computed. Concerning the sensitivity analysis, results confirm the findings of previous work; the uncertainty in the outlet static pressure has the highest impact on the loss coefficient and is predominant with respect to the thermodynamic model uncertainties; the inlet total conditions, instead, dominate on the mass flow rate. With respect to the present state-of-the art, a novel result has been observed, i.e. the impact of operational and model uncertainties is drastically reduced when manufacturing tolerances are considered. The modal representation of manufacturing uncertainties and the related statistical analysis allow to highlight the most relevant modes, suggesting that specific areas of the blade have the highest impact on the overall uncertainty. The support of the probability density functions of the quantities of interest enlarges significantly when the manufacturing tolerances are plugged into the uncertainty analysis, and the functions take a regular gaussian shape (differently from what observed if only operative and model uncertainties are considered). On a distributed basis, the flow released by the cascade exhibits relevant sensitivity to all sources of uncertainty, but again its statistical moments grow significantly in presence of manufacturing tolerances; as the dominating modes alter the blade shape in the rear suction side and the trailing edge thickness, the strength and inclination of the main rear shock undergoes significant variability, so as the width and defect of the blade wake.

The behaviour characterized in this work could be a very relevant statement for design. Conclusions drawn will constitute the basis for future studies about robust optimization of supersonic turbine cascades, e.g. minimizing one of the statistics of the Quantity of Interest (such as its mean or quantile). Moreover, considering the potential interest in using low-fidelity computations, we aim to include multi-fidelity approaches for controlling the loss of accuracy generated by the use of a low-fidelity model in the optimization loop.

  8.0 bars 543.65 K 1.33 bars Part-Load 8.0 bars 543.65 K 2.0 bars

Figure 1 :

 1 Figure 1: 180k cells mesh.

Figure 2 :

 2 Figure 2: Grid Analysis on the baseline configuration (RANS), for the two QoI ṁ and Y.

Figure 4 :

 4 Figure 4: Eigenmodes Illustration, with scale=50, σ 0 = 6 × 10 -5 .

  scale=10, σ 0 = 6 × 10 -5 .

Figure 5 :

 5 Figure 5: Geometric Variability Modeling: Samples.

Figure 6 :

 6 Figure 6: PDF Convergence analysis, based on LHS sets of size 100, 200 and 400. Stars indicate 5% and 95% quantiles, 'plus' indicate µ ± σ and diamonds indicate µ.

Figure 7 :

 7 Figure 7: Total Sobol SI Convergence analysis, based on LHS sets of size 100, 200 and 400.

  (a) Full Load: σ 0 = 0 (b) Full Load: σ 0 = 3 × 10 -5 (c) Full Load: σ 0 = 6 × 10 -5 (d) Part Load: σ 0 = 0 (e) Part Load: σ 0 = 3 × 10 -5 (f) Part Load: σ 0 = 6 × 10 -5

Figure 8 :

 8 Figure 8: Sensitivity Analysis Y: Total Sobol SI

  (a) Full Load: σ 0 = 0 (b) Full Load: σ 0 = 3 × 10 -5 (c) Full Load: σ 0 = 6 × 10 -5 (d) Part Load: σ 0 = 0 (e) Part Load: σ 0 = 3 × 10 -5 (f) Part Load: σ 0 = 6 × 10 -5

Figure 9 :

 9 Figure 9: Sensitivity Analysis ṁ: Total Sobol SI.

  (a) Full Load: σ 0 = 0 (b) Full Load: σ 0 = 3 × 10 -5 (c) Full Load: σ 0 = 6 × 10 -5 (d) Part Load: σ 0 = 0 (e) Part Load: σ 0 = 3 × 10 -5 (f) Part Load: σ 0 = 6 × 10 -5

Figure 11 :

 11 Figure 11: Joint PDF ( ṁ,Y)

  (a) Mean: σ 0 = 0 (b) Mean: σ 0 = 3 × 10 -5 (c) Mean: σ 0 = 6 × 10 -5 (d) CoV: σ 0 = 0 (e) CoV: σ 0 = 3 × 10 -5 (f) CoV: σ 0 = 6 × 10 -5

Figure 12 :

 12 Figure 12: Mean and CoV[%] Mach contours for full load.

  (a) Mean: σ 0 = 0 (b) Mean: σ 0 = 3 × 10 -5 (c) Mean: σ 0 = 6 × 10 -5 (d) CoV: σ 0 = 0 (e) CoV: σ 0 = 3 × 10 -5 (f) CoV: σ 0 = 6 × 10 -5

Figure 13 :

 13 Figure 13: Mean and CoV[%] Mach contours for part load.

Table 1 :

 1 Gas properties of the siloxane MDM

	Critical pressure	14.152 bar
	Critical temperature	564.1 K
	Critical density	256.82 kg.m -3
	γ	1.0165
	Acentric factor ω	0.529
	Gas constant	35.152 J/kg/K
	µ	1.1517 × 10 -5 Pa.s
	k	0.03799 W/(m.K)

Table 2 :

 2 Operating Conditions

Table 3 :

 3 Full Load: Uncertain Input, for different geometric variability levels

	Geometric Variations	No	Moderate	Large
		σ 0 = 0	σ 0 = 3 × 10 -5	σ 0 = 6 × 10 -5
	P s out [bars] P t in [bars] T t in [K]	U[1.28, 1.38] U[7.95, 8.05] U[543.3, 544.0] U[543.3, 544.0] U[543.3, 544.0] U[1.28, 1.38] U[1.28, 1.38] U[7.95, 8.05] U[7.95, 8.05]
	γ [-]	U[1.012, 1.020] U[1.012, 1.020] U[1.012, 1.020]
	ω [-]	U[0.518, 0.540] U[0.518, 0.540] U[0.518, 0.540]
	ξ	-	N(0, I 7 )	N(0, I 7 )
	Dimension θ	5	12	12
	Geometric Variations	No	Moderate	Large
		σ 0 = 0	σ 0 = 3 × 10 -5	σ 0 = 6 × 10 -5
	P s out [bars] P t in [bars] T t in [K]	U[1.95, 2.05] U[7.95, 8.05] U[543.3, 544.0] U[543.3, 544.0] U[543.3, 544.0] U[1.95, 2.05] U[1.95, 2.05] U[7.95, 8.05] U[7.95, 8.05]
	γ [-]	U[1.012, 1.020] U[1.012, 1.020] U[1.012, 1.020]
	ω [-]	U[0.518, 0.540] U[0.518, 0.540] U[0.518, 0.540]
	ξ	-	N(0, I 7 )	N(0, I 7 )
	Dimension θ	5	12	12

t in , T t in , P s out , γ, ω, ξ 1 , ..., ξ N σ ), or θ = (P t in , T t in , P s out , γ, ω) if geometric variations are not taken into account.

Table 4 :

 4 Part Load: Uncertain Input, for different geometric variability levels

(a) Full Load (b) Part LoadFigure 3: Mach contours: Nominal Conditions

does not hold only for the mean properties but it is valid for the whole stochastic framework. This result, which is physically consistent, might be considered as a further proof of the validity of the proposed statistic methodology based on surrogates.

Full Load (µ,σ,q 5 ,q 95 ) Y

Physical Analysis

In this section, QoI of high dimension variability is conducted using the UQ tools described in Subsection 4.2 to perform a physical analysis of the flow. The distributions of Mach number over the entire flow domain and its variability are now discussed, by resorting to the mean and the Coefficient of Variance.

The QoI here shown refer to Mach number at node location of the initial mesh, the latter allowed to be deformed conserving the initial nodes connectivity. Prescribing the QoI to assigned spatial location would require a curbersome interpolation/extrapolation, adding numerical inaccuracies of the QoI, in the vicinity of the blade wall, due geometric blade perturbations.

The distributions of mean and CoV of the Mach number on the computational domain is shown Figures 12 and13, corresponding respectively to the full-and part-load conditions. The mean Mach number distributions are actually very similar each other and to the deterministic ones. This is probably because the uncertainties are relatively small and are symmetric with respect to the deterministic values; moreover, the manufacturing tolerances are also assumed


Xiaofei Li
email: xiaofei.li@inria.fr

Laurent Girin

Radu Horaud

Expectation-Maximization for Speech Source Separation using Convolutive Transfer Function

This paper addresses the problem of under-determinded speech source separation from multichannel microphone singals, i.e. the convolutive mixtures of multiple sources. The time-domain signals are first transformed to the short-time Fourier transform (STFT) domain. To represent the room filters in the STFT domain, instead of the widely-used narrowband assumption, we propose to use a more accurate model, i.e. the convolutive transfer function (CTF). At each frequency band, the CTF coefficients of the mixing filters and the STFT coefficients of the sources are jointly estimated by maximizing the likelihood of the microphone signals, which is resolved by an Expectation-Maximization (EM) algorithm. Experiments show that the proposed method provides very satisfactory performance under highly reverberant environments.

Introduction

Most of speech source separation techniques are designed in the short time Fourier transform (STFT) domain where the narrowband assumption is generally used, e.g. [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF][START_REF] Mandel | Model-based expectation-maximization source separation and localization[END_REF][START_REF] Winter | MAPbased underdetermined blind source separation of convolutive mixtures by hierarchical clustering and 1 -norm minimization[END_REF][START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF]. In the narrowband assumption, at each frequency band, the time-domain filter is represented by the acoustic transfer function (ATF), and the time-domain convolutive process is transformed to a product between the ATF and the STFT coefficients of the source signal. This assumption is also referred to as the multiplicative transfer function (MTF) approximation [START_REF] Avargel | On multiplicative transfer function approximation in the short-time Fourier transform domain[END_REF]. Based on the ATF or its variant, e.g. relative transfer functions [START_REF] Gannot | Signal enhancement using beamforming and nonstationarity with applications to speech[END_REF][START_REF] Li | Estimation of relative transfer function in the presence of stationary noise based on segmental power spectral density matrix subtraction[END_REF], beamforming techniques are widely used for multichannel speech source separation and speech enhancement. Popular beamformers for multisource separation include linearly constrained minimum variance/power (LCMV/LCMP) [START_REF] Gannot | Signal enhancement using beamforming and nonstationarity with applications to speech[END_REF][START_REF] Van Trees | Detection, estimation, and modulation theory[END_REF]. Furtherly, because of the spectral sparsity of speech, the microphone signal can be assumed to be dominated by only one speech source in each timefrequency (TF) bin. This is refered to as the W-disjoint orthogonality (WDO) assumption [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF]. The binary masking method [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF][START_REF] Mandel | Model-based expectation-maximization source separation and localization[END_REF] and 1norm minimization method [START_REF] Winter | MAPbased underdetermined blind source separation of convolutive mixtures by hierarchical clustering and 1 -norm minimization[END_REF] exploit such WDO assumption. More examples of narrowband assumption-based techniques can be found in [START_REF] Gannot | A consolidated perspective on multimicrophone speech enhancement and source separation[END_REF] and references therein.

In real scenarios, the time-domain filter, i.e. the room impulse response (RIR), is normally much longer than the STFT window (frame) length, since the latter should be set to be sufficiently short to account for the local stationarity of speech. For this case, the narrowband assumption is no longer valid, and thus leads to unsatisfied speech source separation performance. In the literature, only a few studies had questioned the validity of the narrowband assumption and attempted to tackle this problem. Based on the narrowband assumption, the theoretical covariance matrix of one source image is a rank-one matrix [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF]. To mitigate the invalidity of this matrix in practice, a full-rank spatial covariance matrix was adopted in [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], even if the narrowband assumption is used. To circumvent the inaccurancy of the narrowband assumption, the wideband timedomain convolution model was used in [START_REF] Kowalski | Beyond the narrowband approximation: Wideband convex methods for underdetermined reverberant audio source separation[END_REF][START_REF] Arberet | Sparse reverberant audio source separation via reweighted analysis[END_REF], where the source STFT coefficients are recovered by minimizing the fit cost between the time-domain mixtures and the time-domain convolution model. Meanwhile, based on the Lasso technique, the 1 -norm of the source STFT coefficients is minimized as well to impose the sparsity of speech spectra. This method achieves good performance, but its computational complexity is very large due to the high cost of timedomain convolution operation. In [START_REF] Leglaive | Multichannel audio source separation: variational inference of time-frequency sources from time-domain observations[END_REF][START_REF] Leglaive | Separating timefrequency sources from time-domain convolutive mixtures using non-negative matrix factorization[END_REF], based on the criterion of likelihood maximization, a variational EM algorithm was proposed also using the time-domain convolution model and the STFT domain signal model.

In [START_REF] Avargel | System identification in the short-time Fourier transform domain with crossband filtering[END_REF], the time-domain convolution can be ideally represented in the STFT-domain by the cross-band filters. More precisely, the STFT coefficients of the source image can be computed by summing multiple convolutions (over frequencies) between the STFT coefficients of the source signal and the STTF-domain filter. Note that the convolution is conducted along the frame axis. To simplify the analysis, for each frequency, the band-to-band filter, i.e. the CTF model [START_REF] Talmon | Relative transfer function identification using convolutive transfer function approximation[END_REF], is used, while with the cross-band information omitted. Compared to the narrowband assumption that uses a frequency-wise scalar product to approximate the time-domain convolution, the CTF model uses a frequency-wise convolution and thus is more accurate. Following the principle of the wide-band Lasso [START_REF] Kowalski | Beyond the narrowband approximation: Wideband convex methods for underdetermined reverberant audio source separation[END_REF], based on the CTF model, a subband Lasso technique was proposed in [START_REF] Li | Audio source separation based on convolutive transfer function and frequency-domain lasso optimization[END_REF], which largely reduces the complexity relative to the wide-band Lasso technique. In [START_REF] Li | Multichannel speech separation and enhancement using the convolutive transfer function[END_REF], two CTF inverse filtering methods were proposed based on the multiple-input/output inverse theorem (MINT). In [START_REF] Talmon | Convolutive transfer function generalized sidelobe canceler[END_REF], the CTF was integrated into the generalized sidelobe canceler beamformer. A CTF-based EM algorithm was proposed in [START_REF] Schwartz | Online speech dereverberation using kalman filter and EM algorithm[END_REF] for single-source dereverberation, in which the Kalman filter was exploited to achieve online EM update. The cross-band filters were adopted in [START_REF] Badeau | Multichannel high-resolution NMF for modeling convolutive mixtures of non-stationary signals in the time-frequency domain[END_REF], combined with a non-negative matrix factorization model for the source signal. To estimate the source signals, the likelihood of the microphone signals is maximized via a variational EM algorithm. In [START_REF] Higuchi | Joint audio source separation and dereverberation based on multichannel factorial hidden Markov model[END_REF], also based on likelihood maximization and EM, a STFT-domain convolutive model was used for source separation, combined with an HMM model for source activity. Even though this STFT-domain convolutive model was not named as CTF in [START_REF] Higuchi | Joint audio source separation and dereverberation based on multichannel factorial hidden Markov model[END_REF], it actually plays the same role as CTF.

Due to the high model complexity, the above mentioned source separation techniques that are beyond the narrowband assumption actually can not be performed in a blind manner, in other words, some prior knowledge are required. For example, the RIRs for the wide-band Lasso techniques, and the CTFs for the CTF-Lasso and MINT techniques, are required to be knwon or well estimated. Both the mixing filters and source parameters are required as a good initialization for the (variational) EM techniques [START_REF] Leglaive | Multichannel audio source separation: variational inference of time-frequency sources from time-domain observations[END_REF][START_REF] Leglaive | Separating timefrequency sources from time-domain convolutive mixtures using non-negative matrix factorization[END_REF][START_REF] Badeau | Multichannel high-resolution NMF for modeling convolutive mixtures of non-stationary signals in the time-frequency domain[END_REF][START_REF] Higuchi | Joint audio source separation and dereverberation based on multichannel factorial hidden Markov model[END_REF]. A blind multichannel CTF identification method was proposed in [START_REF] Li | Multichannel identification and nonnegative equalization for dereverberation and noise reduction based on convolutive transfer function[END_REF] and the identified CTF can be fed in the semi-blind methods. However, this CTF identification method was only suitable for the single-source case.

In the present paper, based on the CTF model, we propose a likelihood maximization method for speech source separation. First, the CTF model is presented in a source mixture probabilistic framework. Then, an EM algorithm is proposed for resovling the likelihood maximization problem. The STFT coefficients of the source signals are taken as hidden variables, and are estimated in the expectation step (E-step). The CTF coefficients and source parameters are estimated in the maximization step (M-step). Experiments show that the proposed method performs better than the narrowband assumpation based methods [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF][START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF] and the CTF-lasso method [START_REF] Li | Audio source separation based on convolutive transfer function and frequency-domain lasso optimization[END_REF] within a semiblind setup where the mixing filters are initialized with a perturbed version of the ground-truth CTF.

The rest of this paper is organized as follows. Section 2 presents the CTF formulation, which is plugged in a probabilistic framework in Section 3. The proposed EM is given in Section 4. Experiments are presented in Section 5. Section 6 concludes the paper. This paper is an extension of a conference paper [START_REF] Li | An EM algorithm for audio source separation based on the convolutive transfer function[END_REF]. The main improvements over [START_REF] Li | An EM algorithm for audio source separation based on the convolutive transfer function[END_REF] consists of i) we present the methodology in more detail, such as the two vector/matrix formulations in Section 3, the detailed derivation of the EM algorithm in Section 4 and the execution process of EM in Algorithm 1; ii) in Section 5, we add the experiments with CTF perturbations, and analyze the computational complexity of the proposed method.

CTF Formulation

In an indoor (room) environment, a speech source signal propagates to the receivers (microphones) through the room effect. In the time domain, the received source image y(n) is given by the linear convolution between the speech source signal s(n) and the RIR a(n):

y(n) = a(n) s(n), (1)
where denotes convolution. Applying STFT, the narrowband approximation usually approximates the time domain convolution (1) as the product y(p, k) = a(k)s(p, k), where y(p, k) and s(p, k) are the STFT coefficients of the corresponding signals, and a(k) is the ATF. Let N denote the STFT window (frame) length, then k ∈ [0, N -1] is the frequency bin index, and p ∈ [1, P ] is the frame index. This narrowband approximation is no longer valid when a(n) is long relative to the STFT window.

We use the CTF model to circumvent the inaccurancy of the narrowband assumption, then y(p, k) can be presented as [START_REF] Talmon | Relative transfer function identification using convolutive transfer function approximation[END_REF]:

y(p, k) = a(p, k) s(p, k) = Q p =0 a(p , k)s(p -p , k), (2)
where a(p, k) denotes the CTF, which can be derived from the timedomain filter a(n) by:

a(p, k) = (a(n) ζ k (n))| n=pL . ( 3
)
This equation is the convolution with respect to the time index n evaluated at multiples of the frame step L, with

ζ k (n) = e j 2π N kn +∞ m=-∞ ωa(m) ωs(n + m),
where ωa(n) and ωs(n) are respectively the STFT analysis and synthesis windows. The length of CTF, i.e. Q + 1, approximately equals the length of RIR divided by L.

Mixture Model Formulations

Basic Formuation for Mixture Model

We consider a source separation problem with J sources and I sensors, which could be either underdetermined (I < J) or (over)determined (I ≥ J). Using the CTF formulation (2), in the STFT domain, the microphone signal x i (p, k) is:

x i (p, k) = J j=1 a ij (p, k) s j (p, k) + e i (p, k), (4)
where a ij (p, k) is the CTF from source j, j = 1, . . . , J to sensor i, i = 1, . . . , I, and e i (p, k) denotes the noise signal.

Probabilistic Model

In the literature of source separation, each source signal s j (p, k) is normally assumed to be independent to other sources, and is also independent across STFT frames and frequencies. Each STFT coefficient s j (p, k) is assumed to follow a complex Gaussian distribution with a zero mean and variance v j (p, k) [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF][START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], i.e. its probability density function (pdf) is:

Nc(s j (p, k); 0, v j (p, k)) = 1 πv j (p, k) exp - |s j (p, k)| 2 v j (p, k) .
The noise signal is assumed to be stationary, temporally uncorrelated, and independent to the speech source signals. We define the noise vector accross microphones as e(p, k) = [e 1 (p, k), . . . , e I (p, k)] ∈ C I×1 . This vector is assumed to follow a zero-mean complex Gaussian, with a non-diagonal covariance matrix denoted as Σe(k). This convariance matrix encodes the spatial correlation of the noise signals. The pdf is:

Nc(e(p, k); 0, Σe(k)) = 1 π I |Σe(k)| e -e(p,k) H Σe(k) -1 e(p,k) ,
where H denotes complex transpose, | • | the determinant of matrix.

Since the proposed source separation method is carried out independently at each frequency, hereafter, we omit the frequency index k for notational simplicity.

Vector/matrix Formulation 1

To formulate the mixture model (4) more compactly, we have several different choices to organize the signals and the convolution operation in vector/matrix forms. To facilitate the derivation of the following EM algorithm, we use two different vector/matrix formulations. In this section, Formulation 1 will be presented to enable us to easily derive the M-step, and Formulatioin 2 will be presented used for the derivation of the E-step in the next section. This two formulations are different just in the organization of the variables and parameters, thence transforming from the E-step to the M-step, and vice-versa, will only necessitate reorganizing the vector/matrix elements.

In Formulation 1, we define the source signals in vector form as, for p ∈ [1, P ]:

s j (p) = [s j (p), . . . , s j (p -q), . . . , s j (p -Q)] ∈ C (Q+1)×1 , s(p) = [s 1 (p) , . . . , s j (p) , . . . , s J (p) ] ∈ C J(Q+1)×1 ,
where denotes vector/matrix transpose. If p ≤ q, we set s j (pq) = 0. Define the CTF in vector/matrix form as:

a ij = [a ij (0), . . . , a ij (p), . . . , a ij (Q)] ∈ C (Q+1)×1 , A j = [a 1j ; . . . ; a ij ; . . . ; a Ij ] ∈ C I×(Q+1) , A = [A 1 , . . . , A j , . . . , A J ] ∈ C I×J(Q+1) ,
We already defined e(p) ∈ C I×1 . Similarly, the microphone signal is x(p, k) = [x 1 (p, k), . . . , x I (p, k)] ∈ C I×1 . Finally, we can rewrite (4) as:

x(p) = J j=1 A j s j (p) + e(p) = As(p) + e(p). (5)
In Formulation 1, the source vector s(p) follows a zeromean complex Gaussian distribution with J(Q + 1) × J(Q + 1) diagonal covariance matrix Rs(p) where the first Q + 1 diagonal entries (for the first source) are v 1 (p), . . . , v 1 (P -Q), the next Q + 1 diagonal entries (for the second source) are v 2 (1), . . . , v 2 (P ), and so on. The pdf of the mixture given the sources is Nc(x(p); As(p), Σe).

Vector/matrix Formulation 2

In Formulation 2, let sj = [s j (1), . . . , s j (p), . . . , s j (P )] T ∈ C P ×1 , ẽi = [e i (1), . . . , e i (p), . . . , e j (P )] T ∈ C P ×1 , xi = [x i (1), . . . , x i (p), . . . , x j (P )] T ∈ C P ×1 denote the j-th source vector and the i-th noise and microphone vectors, all involving all P frames. Concatenate them along source/microphone as:

s = [s 1 , . . . , s j , . . . , s J ] ∈ C JP ×1 , ẽ = [ẽ 1 , . . . , ẽ i , . . . , ẽ I ] ∈ C IP ×1 , x = [x 1 , . . . , x i , . . . , x I ] ∈ C IP ×1 .
Define the CTF convolution matrix in C P ×P :

A ij =          a ij (0) 0 • • • • • • • • • 0 . . . . . . . . . . . . . . . . . . a ij (Q) . . . a ij (0) 0 . . . 0 . . . . . . . . . . . . . . . . . . 0 • • • 0 a ij (Q) • • • a ij (0)         
, where the CTF {a ij (p)} is first flipped and then duplicated as the row vectors, with one element shift per row. Concatenate it along source and microphone as:

A i = [A i1 , . . . , A ij , . . . , A iJ ] ∈ C P ×JP , A = [A 1 , . . . , A i , . . . , A I ] ∈ C IP ×JP . (6)
Then we can rewrite (4) as:

xi = J j=1 A ij sj + ẽi = A i s + ẽi , or x = As + ẽ. (7)
In Formulation 2, the pdf of source vector s is a zero-mean complex Gaussian distribution with JP × JP diagonal covariance matrix

Ψs = Diag{[v 1 (1), . . . , v 1 (P ), . . . , v J (1), . . . , v J (P )] } (8)
where Diag{•} denotes diagonal matrix of a vector. The noise vector ẽ follows a zero-mean complex Gaussian distribution with IP × IP covariance matrix Ψe. The entries of Ψe are

Ψe((i 1 -1)P + p 1 , (i 2 -1)P + p 2 ) = Σe(i 1 , i 2 ), if p 1 = p 2 , 0, otherwise. ( 9
)
where the arguments in parentheses denotes the row and column indices;

i 1 , i 2 ∈ [1, I], p 1 , p 2 ∈ [1, P ].
The pdf of the mixture given the sources is Nc(x; As, Ψe).

Expectation-Maximization Algorithm

Collect the source variances for all sources and frames, we have the source variance set V = {v j (p)} j∈[1,J],p∈ [1,P ] . The parameters set in the present problem is Θ = {V, A, Σe} in Formulation 1 or Θ = {Ψs, A, Ψe} in Formulation 2. The likelihood of the mixture will be maximized by an EM algorithm, in which the parameters Θ are the optimization variables. Meanwhile, the the STFT coefficients of the source signals, i.e. {s j (p)} j,p , are taken as hidden variables, whose posterior statistics will be inferred, and the posterior mean is taken as the estimation of the source signals. The proposed EM algorithm is summarized in Algorithm 1.

E-step

The E-step will be derived based on Formulation 2. Using the parameters estimates Θ, given in the preceding M-step by [START_REF] Talmon | Relative transfer function identification using convolutive transfer function approximation[END_REF] in Formulation 1, we construct the CTF convolution matrix A, the source covariance matrix Ψs and the noise covariance matrix Ψe following ( 6), ( 8) and ( 9), respectively.

In Formulation 2, the posterior distribution of the source signals is p(s|x, Θ) ∝ p(x|s, Θ)p(s|Θ). Since both p(x|s, Θ) and p(s|Θ) are Gaussian, p(s|x, Θ) is also Gaussian. Let E s| x,Θ [ . ] denotes the expectation in the sense of the posterior distribution p(s|x, Θ). From the exponent of p(s|x, Θ), i.e.

-(x -As) H Ψe -1 (x -As) -sH Ψs -1 s, the posterior mean s = E s| x,Θ [s] and covariance matrix Σs = E s| x,Θ [(s -s)(ss) H ] can be derived by reorganizing the quadratic and linear forms in s . We obtain:

Σs = A H Ψe -1 A + Ψs -1 -1 , s = ΣsA H Ψe -1 x, (10)
and then the posterior second-order moment matrix Rs = E s| x,Θ [ss H ] can be computed as

Rs = s sH + Σs. (11)
As derived based on the narrowband assumption in [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF], Eq. ( 10) is also in the form of classical Wiener filtering. The difference is, here the interframe elements in the posterior covariance matrix Σs are nonzero, which means the correlation between frames due to the convolution is encoded. As a result, the posterior mean of source, i.e. s, is recovered by deconvoluting the mixture.

M-step

The M-step will be derived based on Nc(s j (p); 0, v j (p)).

Considering only the terms related to the parameters and hidden variables, the corresponding loglikelihood writes:

Denote the auxiliary function for likelihood maximization as

Q(Θ, Θ old ) = E S|X,Θ old [log(p(X, S|Θ))],
where Θ old denotes the parameters estimated at the previous iteration. From the loglikelihood [START_REF] Arberet | Sparse reverberant audio source separation via reweighted analysis[END_REF], the auxiliary function can be derived as:

Q(Θ, Θ old ) = - P p=1 log(|Σe|) + Trace{Σe -1 (A s(p)x(p) H + x(p) s(p) H A H -A Rs(p)A H )} - J j=1 P p=1 log(v j (p)) + v j (p) v j (p) + const, ( 13
)
where Trace{•} denotes matrix trace, and

s(p) = E S|X,Θ old [s(p)], Rs(p) = E S|X,Θ old [s(p)s(p) H ], v j (p) = E S|X,Θ old [|s j (p)| 2 ]
are the posterior statistics of the source signal, namely the posterior mean, the posterior second-order moment matrix and the elementwise posterior second-order moment, respectively. Actually, v j (p) is the ((j -1)(Q + 1) + 1)-th diagonal entry of Rs(p). Reformulation: s(p) and Rs(p) can be obtained by reformulating s and Rs derived in the preceding E-step. The reformulation is mainly to find the elements with the same source and frame indices:

s(p) {j,p-q} = s{j,p-q} Rs(p) {j1,p-q1},{j2,p-q2} = Rs {j1,p-q1},{j2,p-q2} , (14)
where the subscript {j,p} denote "the j-th source at p-th frame". In a vector, or in a row/column of a matrix, {j,p} represents i) the {(j -1)P + p}-th element in s and Rs, based on Formulation 2, and ii) the {(j -1)(Q + 1) + q + 1}-th (q ∈ [0, Q]) element in s(p) and Rs(p), based on Formulation 1.

With respect to A * ( * denotes conjugate), v j (p) and Σe, the (complex) derivative of Q(Θ, Θ old ) are respectively:

∂Q(Θ, Θ old ) ∂A * = -Σe -1 P p=1 x(p) s(p) H -A Rs(p) , ∂Q(Θ, Θ old ) ∂v j (p) = v j (p)v -2 j (p) -v -1 j (p), ∂Q(Θ, Θ old ) ∂Σe = P p=1 Σe -1 (A s(p)x(p) H + x(p) s(p) H A H -A Rs(p)A H )Σe -1 -Σe -1 (15)
To maximize Q(Θ, Θ old ), the three derivatives are set equal to zero, then A, v j (p) and Σe can be estimated as, respectively:

A new = P p=1 x(p) s(p) H P p=1 Rs(p) -1 , v new j (p) = v j (p). Σe new = 1 P P p=1 A new s(p)x(p) H + x(p) s(p) H A new H -A new Rs(p)A new H (16)
5 Experiments

Experimental Configuration

The binaural (two-channel) simulated signals were used to evaluate the proposed EM algorithm. The experiments were conducted under various acoustic conditions, in terms of reverberation time, number of sources and intensity of room filter perturbation.

Simulation set-up:

We use a KEMAR dummy head [START_REF] Gardner | HRTF measurements of a KEMAR dummy-head microphone[END_REF] with one microphone embedded in each ear as the recording system. The head related impulse responses (HRIRs) for a large grid of directions were measured in advance. The ROOMSIM simulator [START_REF] Campbell | The roomsim user guide (v3.3)[END_REF] simulates the binaural room impulse responses using these HRIRs as both the direct-path wave and the reflections. Four reverberant conditions were simulated, with the reverberation time T 60 as 0 s (the anechoic case), 0.22 s, 0.5 s and 0.79 s, respectively. The TIMIT [START_REF] Garofolo | Getting started with the DARPA TIMIT CD-ROM: An acoustic phonetic continuous speech database[END_REF]

EM initialization:

Depends on what types of prior knowledges are available, the EM algorithm can be initialized either from the E-step or from the M-step. For both choices, the accuracy of the initialization is crucial for the EM iterations to converge to a good solution. In this experiment, we initialize EM from the M-step. Due to the difficulty of the blind initialization, we consider a semiblind initialization scheme. The time-domain filters are assumed to be known, from which the CTFs are computed by [START_REF] Winter | MAPbased underdetermined blind source separation of convolutive mixtures by hierarchical clustering and 1 -norm minimization[END_REF]. To fit the realistic situation that the time-domain filters (or CTFs) should actually be blindly estimated and suffer from some estimation error, a proportional Gaussian random noise is added to the time-domain filters to generate the perturbed filters. The normalized projection misalignment (NPM) [START_REF] Morgan | On the evaluation of estimated impulse responses[END_REF] in decibels (dB) is used to measure the intensity of the perturbation. The lower NPM is, the intenser the perturbation will be. To have a good initialization for the source variance, the CTF-Lasso method proposed in [START_REF] Li | Audio source separation based on convolutive transfer function and frequency-domain lasso optimization[END_REF] is first applied, which solves the problem:

min s As -x 2 2 +λ s 1 ,
where • 1 denotes 1 -norm, and is used to impose the spectral spasity of speech sources. For more details, please refer to [START_REF] Li | Audio source separation based on convolutive transfer function and frequency-domain lasso optimization[END_REF]. Then the source variance is initialized as the magnitude square of each source coefficient estimate. It is found that, for most of acoustic conditions and frequency bins, the number of EM iterations required for convergency is less than 10, thence in this experiment we use a constant number, i.e. 7, of EM iterations.

5.1.3

Baseline methods: Three baseline methods were used for comparison: i) The CTF-Lasso method used for initialization;

ii) The binary masking (BM) method [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF], which is based on the narrowband approximation. To make a fair comparison, the narrowband mixing filters are also computed using the knwon perturbed time-domain filters. However, to compute the one-order mixing filters, for the high reverberation time cases, the time-domain filters should be first truncated to have a length being equal to (or less than) the STFT window length. Based on some pilot experiments, we use the HRIRs (without reverberation) as the truncated filters, which achieves the best results. For source separation, each TF bin is assigned to one of the sources based on the mixing filters; iii) the full-rank spatial covariance matrix (FR-SCM) method [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF]. The full-rank spatial covariance matrix for each source was separately estimated using the corresponding source image, and kept fixed during the EM, following the line of the semi-oracle experiments in [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF].

Performance metrics:

The signal-to-distortion ratio (SDR), signal-to-interferences ratio (SIR) and signal-to-artifacts ratio (SAR) [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF], all in decibels (dB), are used as the separation performance metric. In the following, three sets of experiments will be conducted: i) for various reverberation time, ii) for various numbers of sources, i.e. with 2, 3, 4 and 5 sources, iii) with various NPM settings. For each condition, the metric scores are averaged over 20 mixtures.

The computation complexity of each method is measured with the real-time factor, which is the processing time of one method divided by the length of the processed signal. Note that all the methods are implemented in MATLAB.

Results as a Function of Reverberation Time

Fig. 1 plots the performance measures obtained for the 4 reverberation times. The number of sources is 3. NPM is set to -35 dB, for which the filter pertrubation is light, and the CTFs used for CTF-Lasso and the proposed method are very accurate. Thence, for this experiment, the CTFs are fixed during EM iterations. For the anechoic case, all the four methods largely improve the SDR and SIR scores compared to the unprocessed signals, however slightly reduces the SAR score. This indicates that the three methods can efficiently separate the multiple sources, but introduces some artifacts. Especially, BM suffers from more artifacts than the other methods, since the hard assignment of the time-frequency bins to the dominant source largely distorts the less dominant sources. As T 60 increases, the performance measures of BM and FR-SCM dramatically decrease, since the length of RIRs is much larger than the STFT window and the narrowband approximation is not valid any more for the high reverberation cases. FR-SCM outperforms BM due to the use of the full-rank convariance matrix, which however is only suitable for the low reverberation cases. This can be testified by the fact that the most prominent advantage of FR-SCM over BM presents when T 60 is 0.22 s. In contrast to BM and FR-SCM, CTF-Lasso and the proposed EM method achieve good performances. It is a bit surprising that the performance measures actually increases with the increase of reverberation time. The possible reason is that the long filters involve more information to differentiate and separate the multiple sources. Of course, to satisfy this assertion the filters should be accurate enough. Compared to CTF-Lasso whose outputs are taken as the initial point for the EM algorithm, the EM algorithm improves the SDR by about 1.5 dB for every reverberation times, which indicates that the EM iterations are able to refine the quality of the source estimate.

Results as a Function of Number of Sources

Fig. 2 plots the results for various number of sources, for T 60 = 0.5 s. NPM is set to -35 dB, and the CTFs are fixed during EM iterations. As expected, the performance measures of all methods degrade when the number of sources increases. For BM, the WDO assumption for speech sources becomes less valid when more sources present. For FR-SCM, CTF-Lasso and the proposed EM method, the mutual confusion of sources gets larger with increasing number of sources. The performance degradation rate of the four methods are similar to the one for the unprocessed signals. Overall, when the CTFs are properly initialized, CTF-Lasso and the proposed EM method achieve good source separtion performance even for the mixtures with five sources using only two microphones.

Results as a Function of NPM

To evaluate the proposed method under the conditions that the initialized CTFs suffer from large estimation error, we conducted the experiments with various NPM settings, and the results are shown in Fig. 3. Since the initialization is not accurate, in this experiment, the CTFs are updated during EM interations to refine the CTFs. To demonstrate the efficiency of CTF updation, the results with the fixed CTFs are also given. As NPM increases, the performance of CTF-Lasso and the proposed EM method largely degrade. When NPM is larger than -20 dB, the proposed method with fixed CTFs does not improve the performance of CTF-Lasso, which means the source estimation can not be refined based the inaccurate CTFs. In contrast, the proposed method with updating CTFs improves the performance of CTF-Lasso, due to the refining of CTFs in the M-step. The performance measures of CTF-Lasso and the proposed method become close to the ones of BM and FR-SCM, which indicates that the CTFbased methods are more sensitive to the filter perturbations than the narrowband assumption based methods.

Computational Complexity Analysis

Table 1 shows the real-time factor of the four methods for the case with three sources and T 60 of 0.5 s. BM is the fastest, since it is an one-step method. The other three are all iterative methods, whose computational complexity is proportional to the number of iterations. FR-SCM and the proposed method are similar in the sense that they both are based on EM iteration: estimating the source statistics using a Wiener-like filter in the E-step, and estimating the mixing filters in the M-step. The main difference between them is that the proposed method uses the CTF with a length of Q, while FR-SCM uses the mixing filter with a length of 1. As a result, the proposed method has a much larger complexity than FR-SCM. CTF-Lasso also uses the CTF. However, unlike FR-SCM and the proposed method that several matrix inverse operations are performed (as shown in [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF] and ( 16)), Lasso optimization only executes the first-order convolution operation. Thence, the real-time factor of CTF-Lasso is lower than the factor of FR-SCM and the proposed method.

Conclusion

In this work, an EM algorithm has been proposed for speech separation. The subband convolutive model, i.e. CTF model, was adopt. To concisely derive the M-step and E-step, two convolution vector/matrix formulations were used. The CTF model based methods, i.e. CTF-Lasso and the proposed EM method, outperform the narrowband assumption based method, i.e. binary masking and FR-SCM, for the reverberant case. The proposed EM algorithm is capable to refine the CTFs and source estimates starting with the output of CTF-Lasso, and thus improves the source separation performance. Only the semi-blind experiments were conducted in this work due to the difficulty of EM initialization. In the future, a blind CTF identification method could be developed to enable the blind initialization of EM. To this aim, the CTF identification methods proposed in [START_REF] Li | Multichannel identification and nonnegative equalization for dereverberation and noise reduction based on convolutive transfer function[END_REF] and [START_REF] Li | Multiple-speaker localization based on direct-path features and likelihood maximization with spatial sparsity regularization[END_REF] could be combined, which are in the contexts of single-source dereverberation and multisource localization, respectively.

Formulation 2 .

 2 Collect the multichannel mixture vectors and source signal vectors along frames, we have the observation set X = {x(p)} p∈[1,P ] and the source signal set S = {s(p)} p∈[1,P ] . The complete-data (including observations and hidden variables) likelihood function is: p(X,S|Θ) ∝ p(X|S, Θ)p(S|Θ)

  log(p(X, S|Θ)) = -P p=1 log(|Σe|) + (x(p) -As(p)) H Σe -1 (x(p) -As(p))j (p)) + |s j (p)| 2 v j (p) + const.

Fig. 1 :

 1 Fig. 1: Performance measures as a function of T 60 . The number of sources is three. NPM is -35 dB. 'unproc.' represents the unprocessed mixture signals.

  speech signals were used as the speech source signals, and were convolved with the simulated BRIRs to generate microphone signals (mixtures). The sampling rates of source signals and microphone signals are both 16-kHz, and the length of source signals are about 3 s. The speech sources were set to locate at different directions in front of the dummy head. The noisy microphone signals is generated by adding a spatially uncorrelated stationary speech-like noise to the noise-free signals. One SNR (signal-to-noise ratio) condition, i.e. 20 dB, is tested. The STFT uses Hamming window with a length of 1,024 samples (64 ms) and frame step of 256 samples (16ms). In this experiment, the ground-truth noise covariance matrix Σe(k) is used, and fixed during EM iterations.

Fig. 2 :

 2 Fig. 2: Performance measures as a function of the number of sources. T 60 = 0.5 s. NPM is -35 dB.

Fig. 3 :

 3 Fig. 3: Performance measures of the proposed methods as a function of NPM. T 60 = 0.5 s. The number of sources is 3. The solid curves with cross marker represent the proposed method with fixed CTFs, while the dashed curves with cross marker represent the proposed method with CTFs updated during EM iterations.

Table 1

 1 Real-time factor of the four methods.

	BM FR-SCM CTF-Lasso Proposed
	0.01	81.7	54.2	630.0

c The Institution of Engineering and Technology 2015