Massimiliano Todisco
email: todisco@eurecom.fr

Héctor Delgado
email: delgado@eurecom.fr

Kong Aik Lee
email: k-lee@ax.jp.nec.com

Nicholas Evans
email: evans@eurecom.fr

Tomi Kinnunen

Junichi Yamagishi
email: jyamagis@nii.ac.jp

Integrated Presentation Attack Detection and Automatic Speaker Verification: Common Features and Gaussian Back-end Fusion

Keywords: automatic speaker verification, spoofing, countermeasures, presentation attack detection

The vulnerability of automatic speaker verification (ASV) systems to spoofing is widely acknowledged. Recent years have seen an intensification in research efforts to develop spoofing countermeasures, also known as presentation attack detection (PAD) systems. Much of this work has involved the exploration of features that discriminate reliably between bona fide and spoofed speech. While there are grounds to use different frontends for ASV and PAD systems (they are different tasks) the use of a single front-end has obvious benefits, not least convenience and computational efficiency, especially when ASV and PAD are combined. This paper investigates the performance of a variety of different features used previously for both ASV and PAD and assesses their performance when combined for both tasks. The paper also presents a Gaussian back-end fusion approach to system combination. In contrast to cascaded architectures, it relies upon the modelling of the two-dimensional score distribution stemming from the combination of ASV and PAD in parallel. This approach to combination is shown to generalise particularly well across independent ASVspoof 2017 v2.0 development and evaluation datasets.

Introduction

Presentation attack detection (PAD) systems capable of detecting and deflecting so-called spoofing attacks, or presentation attack (PA) in ISO/IEC 30107 1 nomenclature, leveled at automatic speaker verification (ASV) systems have been under development for a number of years. While ASV systems aim to verify the identity claimed by a speaker, PAD systems aim to verify the authenticity of the speech signal itself, namely whether it is bona fide speech or whether, instead, it is artificially created or somehow manipulated, i.e. spoofed.

While early PAD systems used features similar to those used for ASV, being distinctly different tasks, most efforts to develop effective PAD systems have focused on the design of new features tailored to discriminate between bona fide and spoofed speech. While the use of features designed specifically for PAD have been shown to give better performance than systems that 1 https://www.iso.org/standard/67381.html use features designed for ASV, the use of different front-ends augments computational complexity.

It can hence be convenient to use a single front-end. The use of such a single front-end avoids redundant processing and can also simplify the combination of ASV and PAD decisions. The search for features which perform well for a combined ASV and PAD task is the subject of this paper.

A second contribution relates to the manner in which ASV and PAD systems scores can be combined. It extends previous work [START_REF] Sahidullah | Integrated spoofing countermeasures and automatic speaker verification: an evaluation on ASVspoof 2015[END_REF] which proposed cascade and parallel approaches to system combination and is similar in nature to the combination architecture reported in [START_REF] Chingovska | Anti-spoofing in action: Joint operation with a verification system[END_REF]. New to this paper is a twodimensional score modelling technique which avoids the joint optimisation of separate ASV and PAD decision thresholds. The explicit modelling of target and impostor trial scores encompassing genuine, bona fide trials in addition to both zeroeffort and spoofed impostor trials provides for greater flexibility in decision boundaries and hence more reliable decisions. The merits of these two contributions are assessed through experiments with the ASVspoof 2017 database of bona fide and spoofed speech signals and protocols for the assessment of combined ASV and PAD systems.

The remainder of the paper is organised as follows. Section 2 describes the different front-ends used in this work. The approach to system combination is presented in Section 3. Experiments are reported in Section 4 whereas results are reported in Section 5. Conclusions are presented in Section 6.

Front-end processing

This paper aims to determine a common front-end for both ASV and PAD tasks. While ASV calls for features that capture speaker-discriminant information, PAD systems rely on features that capture the tell-tale signs of spoofing. The study includes four different front-ends, each of which is described here.

Mel-frequency cepstral coefficients (MFCCs): MFCCs are used widely in speech and speaker recognition and have been explored extensively as features for spoofing detection [START_REF] Sahidullah | A comparison of features for synthetic speech detection[END_REF]. MFCCs are usually derived from short-time Fourier transform (STFT) decompositions, the application of perceptually motivated Mel-frequency scaled filterbank [START_REF] Davis | Comparison of parametric representations for monosyllabic word recognition in continuously spoken sentences[END_REF] and standard cepstral analysis.

Table 1: Classification of trials in PAD and ASV tasks. By "-" we assume that ASV has zero or no capability to reject spoofing imposter trials, similarly for PAD that cannot differentiate zero-effort imposter and target trials. Alternatively "-" means arbitrary: since the last row corresponds to the logical AND operator to combine the two systems, choosing either +1 or -1 in the place of "-" gives rise to the same trial classifications for the integrated system.

class C1 C2 C3 system/trial target zero-effort nontarget spoof nontarget PAD +1 - -1 ASV +1 -1 - ASV + PAD +1 -1 -1
Linear frequency cepstral coefficients (LFCCs): LFCCs are similar to MFCC except for the use of a linear-scaled in place of a Mel-scaled filterbank, thereby giving a constant spectral resolution. LFCCs have also been applied to both speech and speaker recognition, in addition to spoofing detection [START_REF] Sahidullah | A comparison of features for synthetic speech detection[END_REF].

Infinite impulse response constant Q Mel cepstral coefficients (ICMCs): ICMCs have been applied successfully to ASV [START_REF] Delgado | Further optimisations of constant Q cepstral processing for integrated utterance verification and textdependent speaker verification[END_REF], utterance verification (UV) [START_REF] Delgado | Further optimisations of constant Q cepstral processing for integrated utterance verification and textdependent speaker verification[END_REF] and speaker diarization [START_REF] Patino | EURECOM submission to the Albayzin 2016 Speaker Diarization Evaluation[END_REF]. Features are based upon the perceptually-motivated infinite impulse response constant Q transform (IIR-CQT) approach to spectro-temporal decomposition [START_REF] Cancela | An efficient multiresolution spectral transform for music analysis[END_REF]. In contrast to the STFT, the spectral resolution has a constant Q factor which reflects filter selectivity, defined as the ratio between the centre frequency and bandwidth. Efficient feature extraction is obtained from the IIR filtering of the fast Fourier transform (FFT) giving a variable-resolution decomposition with greater frequency resolution at low frequencies and greater time resolution at higher frequencies. ICMCs are then obtained from Mel-scaling and standard cepstral analysis.

Constant Q cepstral coefficients (CQCCs): CQCC features were designed specifically for spoofing detection [START_REF] Todisco | A new feature for automatic speaker verification anti-spoofing: Constant Q cepstral coefficients[END_REF][START_REF] Todisco | Constant Q cepstral coefficients: A spoofing countermeasure for automatic speaker verification[END_REF] and applied subsequently to ASV [START_REF] Todisco | Articulation rate filtering of CQCC features for automatic speaker verification[END_REF]. CQCCs rely on the same CQT approach used in ICMCs but are extracted according to the more computationally demanding approach described in [START_REF] Brown | Calculation of a constant Q spectral transform[END_REF][START_REF] Brown | An efficient algorithm for the calculation of a constant q transform[END_REF]. Resampling [START_REF] Todisco | Constant Q cepstral coefficients: A spoofing countermeasure for automatic speaker verification[END_REF] is applied to warp the geometric scale of the CQT to the linear scale of the discrete Cosine transform (DCT).

Integration of PAD with ASV

The integration of ASV and PAD can be achieved at the model/feature level [START_REF] Sizov | Joint speaker verification and antispoofing in the i-vector space[END_REF] or at the score level [START_REF] Sahidullah | Integrated spoofing countermeasures and automatic speaker verification: an evaluation on ASVspoof 2015[END_REF]. This paper focuses on the latter. Dedicated classifiers are developed for ASV and PAD, and scores produced by each system are combined (i.e., late fusion). Even at this score level, there are different approaches to combination including both cascaded and parallel combinations [START_REF] Sahidullah | Integrated spoofing countermeasures and automatic speaker verification: an evaluation on ASVspoof 2015[END_REF]. We describe below the cascade approach and then the proposed Gaussian back-end fusion approach.

Task and trial definitions

Whereas ASV and PAD systems are both binary classifiers, they tackle different tasks. and trial combination. PAD systems aim to distinguish bona fide speech from spoofed speech, while ASV systems aim to verify a claimed identity. This same idea is illustrated in the two-dimensional score space of Fig. 1. For a threshold θ2 applied to PAD scores sPAD, all trials for which the score is greater than the threshold, sPAD > θ2, are classified as bona fide speech. Those for which the score is below the threshold are classified as presentation attack speech. Similarly, for a threshold θ1 applied to ASV scores sASV, all trials for which the score is greater than the threshold, sASV > θ1, are classified as genuine speaker trials. Those for which the score is below the threshold are classified as impostor trials. The green-colored region to the upper-right-hand corner of Fig. 1, where {sPAD > θ2} ∩ {sASV > θ1}, corresponds to target trials (class C1). Class C2 (zero-effort nontarget) is represented by the blue-colored region, which is the resultant of the impostor and bona fide axes. Similarly, class C3 (spoof nontarget) is represented by the red-colored region, which is the resultant of the target and spoofed axes. As shown in Table 1, only target trials should be positively verified. Both forms of nontarget trial should be rejected.

Cascaded/tandem combination

ASV and PAD systems can be cascaded in either order -PAD followed by ASV, or ASV followed by PAD. ASV and PAD systems can be optimised independently or jointly, e.g. considering an architecture whereby PAD precedes ASV, then the PAD threshold θ2 can be optimised, for instance, to minimise the equal error rate (EER) of the ASV system. In order to estimate the performance of the integrated system, trials classified as spoofs are assigned arbitrarily -∞ scores and are thereby rejected automatically by the ASV system that follows.

The cascaded approach relies on two thresholds, θ2 and θ1, applied to PAD and ASV scores, respectively. The cascading of ASV and PAD in such a way is equivalent to the partitioning of the score space into rigid decision regions with the vertical and horizontal decision boundaries being those illustrated in Fig. 1.

Gaussian back-end fusion

In contrast to the cascaded approach to combination involving the optimisation of two different thresholds for PAD and ASV, the parallel approach to combination requires the optimisation of only a single threshold. Recall the N = 3 classes of trial illustrated in Table 1, namely target (C1), zero-effort nontarget (C2) and spoof nontarget (C3). Let s = [sPAD, sASV] T ∈ R2 be a two-dimensional score vector of PAD and ASV scores corresponding to a single ASV trial.

We treat s as a 2D feature vector and model the classconditional probability density of s using a Gaussian

p(s|C l ) = N (s|µ l , Σ l ),
where µ l ∈ R 2 and Σ l ∈ R 2×2 are the mean vector and the covariance matrix corresponding to class C l where l ∈ {1, 2, 3}. For each trial t, a fused score is computed as the log-likelihood ratio according to:

st = log p(st|H0) p(st|H1) = log p(st|C1) α × p(st|C2) + (1 -α) × p(st|C3) , (1)
where the null hypothesis H0 represents the likelihood that st is a target speaker trial (from class C1) while the alternative hypothesis H1 signifies that st belongs to either C2 or C3. Therefore, (1) converts the 2D detection score vector st into a loglikelihood ratio, a scalar where higher relative values are associated with a stronger support for the null hypothesis. Classification decisions corresponding to the last line of Table 1 may then be made upon the application of a single threshold θ to Note the alternative hypothesis likelihood in the denominator of (1) represents, in fact, a 2-component GMM with mixing weight α. The hyper-parameter α ∈ [0, 1], determines the weight of C2 and C3 in the denominator of Eq. 1 (the alternative hypothesis). Parameters µ l , Σ l and threshold θ can be learned from development data.

Modeling capacity of the two integration approaches

The cascaded approach has two parameters that require optimization -the thresholds θ1 and θ2. The Gaussian back-end approach, in turn, has 3 × 2 = 6 parameters to specify the three mean vectors, 3 × 3 = 9 parameters to specify the three covariance matrices (they are symmetric so 3 parameters are sufficient), one combination parameter α, plus the final decision threshold, i.e. a total of 6 + 9 + 1 = 16 parameters.

While the number of parameters in the Gaussian back-end fusion approach is greater than that for the cascaded approach, generative score modelling offer better potential for generalization. In specific, the decision regions in the cascaded approach have a limited modeling capacity: as illustrated in Fig. 1, the decision region for the target class is an infinite rectangular region, with the lower left corner at (θ1, θ2). Since the joint distribution of the PAD and ASV scores is certainly not defined by such regions, the cascaded approach has limited modeling power. In the Gaussian back-end fusion, in turn, we use three Gaussians with arbitrary (non-shared) covariance matrices, which yields generally a more complex nonlinear decision boundary [START_REF] Duda | Pattern Classification[END_REF] and makes the model better adaptable to arbitrary scores.

Similar 2D back-end fusions have been used earlier in the context of general biometrics and joint operation of PAD and biometric modalities, e.g. [START_REF] Chingovska | Anti-spoofing in action: Joint operation with a verification system[END_REF] and [START_REF] Korshunov | Joint operation of voice biometrics and presentation attack detection[END_REF]. The authors of [START_REF] Chingovska | Anti-spoofing in action: Joint operation with a verification system[END_REF][START_REF] Korshunov | Joint operation of voice biometrics and presentation attack detection[END_REF] used a linear classifier whose parameters (i.e. slope and intercept) were trained using a logistic loss.

Experimental setup

This section describes the database used in this work with implementation details for the individual PAD and ASV systems.

ASVspoof 2017 v2.0 database and protocols

Experiments relate to the ASVspoof 2017 v2.0 database [START_REF] Delgado | ASVspoof 2017 Version 2.0: metadata analysis and baseline enhancements[END_REF] of bona-fide and replayed, spoofed short utterances of about 1 to 5 seconds each. PAD systems are trained using the official training partition containing 1507 bona fide and 1507 replayed speech segments. Joint PAD and ASV experiments were performed using a specifically designed protocol encompassing target segments and both zero-effort and spoofed speech segments 2 . The number of speakers and trials in the development and evaluation sets are illustrated in Table 3.

Front-ends

Experiments were conducted using the four front-ends described in Section 2 using pre-and post-processing comprising the addition of log-energy parameters, cepstral mean and variance normalization (CMVN) [START_REF] Auckenthaler | Score normalization for text-independent speaker verification systems[END_REF], relative spectral (RASTA) filtering [START_REF] Ellis | PLP and RASTA (and MFCC, and inversion) in MAT-LAB[END_REF] and articulation rate (ARTE) filtering [START_REF] Todisco | Articulation rate filtering of CQCC features for automatic speaker verification[END_REF]. Dynamic coefficients up to double deltas are also considered.

The configuration of MFCCs and LFCCs is standard: 19 (S)tatic coefficients (excluding the 0-th), RASTA filtering with appended (D)elta and (A)cceleration coefficients. The ICMC configuration is that reported in [START_REF] Delgado | Further optimisations of constant q cepstral processing for integrated utterance and text-dependent speaker verification[END_REF] and is the same as for MFCCs. Based on configurations used previously for textdependent ASV [START_REF] Todisco | Articulation rate filtering of CQCC features for automatic speaker verification[END_REF] and PAD [START_REF] Delgado | ASVspoof 2017 Version 2.0: metadata analysis and baseline enhancements[END_REF], CQCCs includes 29 S coefficients with appended D coefficients, ARTE filtering and log-energy coefficients. None of the experiments reported here use speech activity detection (SAD). Any single experiment reported in this paper involve PAD and ASV systems that use the exact same front-end configurations.

PAD and ASV systems

Both the PAD and the ASV classifiers are conventional Gaussian mixture models (GMMs). The PAD classifier uses models of 512 components. Models are learned for bona fide and spoofed speech with an expectation-maximisation (EM) algorithm with random initialisation. Classifier scores for a given test utterance are computed as the log-likelihood ratio between the GMMs for bona fide and spoofed speech. The ASV classifier also uses models of 512 components and learns speaker specific models from the maximum a posteriori (MAP) adaptation of a universal background model (UBM) trained on the RSR2015 database [START_REF] Larcher | Text-dependent Speaker Verification: Classifiers, Databases and RSR[END_REF]. Scores are the log-likelihood ratio given the target model and the UBM.

Integration of PAD and ASV

Concerning the Gaussian back-end fusion described in Section 3.3, we use maximum likelihood to obtain the means and covariances of all the three classes (C1, C2 and C3). The value of α is set empirically to 0.96 using a grid search on the development set.

Experimental results and discussion

Results are presented in Table 2 for development (D) and evaluation (E) partitions of the joint PAD-ASV protocol (see Table 3). Results are presented for each front-end and for the cascaded combination and the proposed Gaussian backend fusion (bottom part of Table 2). The performance of two alternative methods, namely linear regression (LR) and polynomial linear regression (PLR) fusion approaches [START_REF] Chingovska | Anti-spoofing in action: Joint operation with a verification system[END_REF], is also reported (top part of Table 2). These contrastive methods also aim to split the 2D score space into two classes. Results are also presented separately for target trials combined with zero-effort and spoofed impostor trials, and the average.

LFCC, ICMC and CQCC features perform marginally better than MFCC features. LFCC features generalise better across development and evaluation subsets, giving the best average (zero-effort and spoofing impostor) results of 13% and 10% EER for each approach to combination. This observation shows that other features that give better performance in terms of spoofing detection do not necessarily give the best performance when ASV and PAD are combined. This is probably due to use of single features which avoids redundant processing and simplifies the combination of ASV-PAD decisions.

These same results also show that the Gaussian back-end fusion approach proposed in this paper outperforms the cascade/tandem combination and the LR and PLR approaches. The improvement in performance is attributed to use of a single, flexible and jointly-optimised, instead of independently optimised, rigid thresholds. The former gives better capacity to reject spoofed trials with less impact on the rejection of target trials. The Gaussian back-end fusion approach reported in this paper therefore offers better robustness to spoofing and better usability. Results for LR and PLR fusion approaches are globally worse with respect to the other two approaches. Moreover, they exhibit similar performance for development and evaluation sets, but with a lack of generalisation.

Finally, Fig. 2 shows the 2D score space representation for the LFCC features for development and evaluation set. By analysing the data distribution, it is clear that a decision can not be taken using approaches with rigid thresholds optimised on development set.

Conclusions

This paper presents a comparative study of the performance of different front-ends when used for both automatic speaker verification (ASV) and presentation attack detection (PAD). Performance is assessed with difference approaches to system integration including cascaded combination, linear and polynomial logistic regression and a Gaussian back-end fusion approach. The use of a single front-end for both ASV and PAD systems simplifies integration in terms of convenience and efficiency; computational effort is reduced by avoiding redundant processing.

Performance is assessed using the ASVspoof 2017 v2.0 database. Results show that feature that achieve the best results for independent ASV or PAD tasks do not give the best performance when systems are combined. The cascaded approach to ASV and PAD combination, as well as the logistic and polynomial logistic regression approach, improve reliability in the case where the nature of non-target trials is known. When faced with unknown, or previously unseen forms of non-target trials and spoofing attacks, then performance degrades significantly; these approaches to ASV and PAD combination fail to generalise. In contrast, the Gaussian back-end approach to integration is shown to generalise well and gives the lowest equal error rate for the independent evaluation set.

Figure 1 :

 1 Figure 1: A two-dimensional view of the score space formed by PAD and ASV scores.

Figure 2 :

 2 Figure 2: Two-dimensional PAD-ASV score representation corresponding to LFCC features for development (left) and evaluation (right) set. Green dots represent target trials, while blue stars represent zero-effort trials and red crosses are spoof trials.

Table 1

 1

defines the three different types of trial that ASV and PAD systems may encounter: (1) target, (2) zero-effort nontarget and (3) spoof nontarget trials. Also illustrated in Table

1

are the ground-truth labels for each task

Table 2 :

 2 Speaker verification performance in terms of EER using linear regression fusion, polynomial logistic regression fusion, cascade/tandem combination and proposed Gaussian back-end fusion of PAD and ASV scores for the ASVspoof 2017 v2.0 database. D: development set, E: evaluation set. The best average results for development set are shown in boldface.

	impostor type	zero-effort	spoof		average		zero-effort	spoof		average
	feat. config / tested on	D	E	D	E	D	E	D	E	D	E	D	E
			Logistic regression fusion [2]			Polynomial logistic regression fusion [2]
	MFCC	3.78	2.42	42.72 31.02	23.25 16.72	3.86	2.50	43.81 35.14	23.84 18.82
	LFCC	5.72	2.11	46.41 35.71	26.06 18.91	5.47	2.20	37.64 26.99	21.55 14.60
	ICMC	2.67	2.16	43.60 33.59	23.14 17.88	2.60	2.08	37.58 29.31	20.09 15.69
	CQCC	6.02	3.52	38.76 33.17	22.39 18.34	6.02	7.93	42.67 47.96	24.34 27.94
			Cascaded/tandem combination			Gaussian back-end fusion
	MFCC	5.19	5.36	23.07 24.65	14.13 15.00	3.99	3.26	21.02 24.35	12.51 13.81
	LFCC	7.28	4.96	22.41 21.28	14.84 13.12	5.71	2.90	21.26 17.98	13.48 10.43
	ICMC	4.30	4.92	27.08 27.82	15.69 16.37	3.06	3.71	20.90 22.51	11.98 13.11
	CQCC	7.31	8.30	15.71 25.26	11.51 16.78	6.04	4.71	17.34 18.11	11.69 11.41

Table 3 :

 3 Statistics of the ASVspoof 2017 joint PAD+ASV protocol.

		#spk target zero-effort spoof
	Dev.	8	742	5186	940
	Eval.	24	1106	18624	10878

Note to the Interspeech reviewers: this custom protocol will be made public later on, to ensure reproducibility.


Andrea Gadda
email: andrea.gadda@polimi.it

Luca Mangani

Giulio Romanelli

Paolo Mantegazza

Ernesto Casartelli

S Y Mpo

O Tatin

A GPU-accelerated compressible RANS solver for Fluid-Structure Interaction simulations in turbomachinery

Keywords:

INTRODUCTION

When dealing with aerodynamic components, especially in the turbomachinery field, the usual approach adopted with CFD relies on the negligible structural deformation due to aerodynamic loads. However, if performance and power efficiency are the ultimate goals (e.g. propulsion), a correct prediction of the structure-flow interaction could improve results accuracy. This is particularly important in a typical geometry optimization loop, where hundreds of slightly different geometries have to be analyzed. Aerolastic simulations, i.e trim and flutter analysis, are nowadays a fundamental step in aircrafts [START_REF] Romanelli | Computational Aeroservoelasticity of Free-Flying Deformable Aircraft[END_REF] and helicopters blades design. In fact, with these aerodynamic components, tip displacements can be easily appreciated as they usually are multiple of the blade thickness. With turbomachinery blades however, deformations due to aerodynamic loads are usually less noticeable. As an example, in our trim simulations of the NASA's Rotor 67 test case, tip displacements are smaller than the blade thickness. Nonetheless, an accurate prediction of this kind of behavior is particularly important due to the presence of the shroud wall near the blade tip. Turbomachinery trim and flutter predictions still represent a challenge due to complex phenomena like rotor-stator interactions, separations and shock waves. The usual time-linearized, frequency-domain strategies can be inadequate when this kind of strong non-linear phenomena occur in the flow, making necessary time-domain simulations [START_REF] Zheng | Coupled fluid-structure flutter analysis of a transonic fan[END_REF][START_REF] Volker Carstens | Coupled simulation of flow-structure interaction in turbomachinery[END_REF] or the harmonic balance technique [START_REF] Earl | Modeling of fluidstructure interaction[END_REF]. The scenario is further complicated by the periodic and multi-stage nature typical of axial compressors and turbines. In flutter and aerodynamic damping simulations, a fundamental role is also played by the Inter-Blade Phase Angle (IBPA) [START_REF] Donini | Aeroelasticity of turbomachines linearized flutter analysis[END_REF] when solution spatial periodicity is not related to a single blade anymore. However, in order to save computational time by reducing the computational domain to a single passage, as usually done with periodic boundary conditions, time-delayed boundary conditions [START_REF] Romanelli | Implementation of a cfd-based aeroelastic analysis toolbox for turbomachinery applications[END_REF] can be used.

Besides unsteady analysis, another important aspect, not yet adequately investigated, is represented by the trim analysis, which is fundamental for an accurate steady analysis that aims to consider static blade elasticity for the performance evaluation of turbomachinery. Trim analysis can be used to provide more accurate initial conditions for unsteady simulations or to improve the characteristic curve computations.

With the aim of an optimization process, another requirement of totally different nature arises: computational efficiency. This problem is typically addressed by implementing an efficient implicit formulation [START_REF] Mangani | Implementation and Evaluation of a Fully Coupled Solver in Openfoam: Steady State Incompressible Turbulent Flows[END_REF][START_REF] Mangani | Implementation and Evaluation of a Fully Coupled Solver in Openfoam: Steady State Incompressible Turbulent Flows in Rotational Reference Frames[END_REF] or by accelerating explicit formulations with multiple combined convergence acceleration techniques [START_REF] Romanelli | Implementation of explicit density-based unstructured cfd solver for turbomachinery applications on graphical processing units[END_REF]. Nowadays, GPUs provide a cheap and relatively easy way to obtain great SIMD (Single Instruction Multiple Data) floating point performance at the cost of a limited amount of available memory with respect to the usual amount of system RAM [START_REF] Favale | A gpu parallelized two fields full potential formulation for real gases[END_REF][START_REF] Romanelli | Implementation of explicit density-based unstructured cfd solver for turbomachinery applications on graphical processing units[END_REF][START_REF] Mangani | Comparison of acceleration techniques on cfd open-source software for aerospace applications[END_REF]. Thus, the idea is to exploit GPUs architecture to accelerate an explicit FSI solver, keeping a low memory usage. Furthermore, the structural FEM model is reduced to a modal representation [START_REF] Romanelli | Coupled csd/cfd non-linear aeroelastic trim of free-flying flexible aircraft[END_REF][START_REF] Romanelli | Computational Aeroservoelasticity of Free-Flying Deformable Aircraft[END_REF] in order obtain both an accurate and efficient FSI formulation.

The purpose of this article is to describe the architecture and the validation of a GPU-accelerated Fluid-Structure Interaction (FSI) solver for compressible viscous (using RANS models) flows. A pioneer work describing the impact of static deformation on the turbomachinery performances will be provided. In particular the well known NASA's Rotor 67 case [START_REF] Arnone | Viscous analysis of three-dimensional rotor flow using a multigrid method[END_REF][START_REF] Grosvenor | Rans prediction of transonic compressive rotor performance near stall[END_REF] is chosen for the trim analysis. This is a typical axial compressor rotor blade. Effects of trimmed solutions on the most important integral quantities (i.e. massflow, characteristic curves) are investigated and a comparison with pure aerodynamic results is provided. Numerical results are also compared with two sets of experimental data available in literature. The GPU solver is written in OpenCL in order to be compatible with both AMD and NVIDIA GPUs. This choice allows the solver to be natively executed also on CPUs, without source code modifications, improving the solver hardware compatibility and allowing to easily perform comparisons with multi-core CPUs. The proposed approach is validated with typical industrial cases.

METHODS

Aerodynamic numerical formulations

The numerical aerodynamic formulations are chosen with the aim to obtain a computationally efficient general-purpose compressible RANS solver. The purpose of this section is to introduce the aerodynamic numerical formulations implemented in the solver. For a better description of the adopted formulations, the reader is referred to [START_REF] Romanelli | Implementation of explicit density-based unstructured cfd solver for turbomachinery applications on graphical processing units[END_REF].

The Reynolds Averaged Navier-Stokes (RANS) equations governing the dynamics of a compressible, viscous and conductive fluid are written within an Arbitrary Lagrangian Eulerian (ALE) framework in integral conservative form as follows:

d dt V U dV + S f(U)-g(U)-U w •n dS = V h(U) dV (1)
where V (t) is the time-varying spatial domain delimited by the boundary S(t) with normal unit vector n(x, t). In the special case of local velocity of the moving boundaries w(x, t) equal to zero, the classical Eulerian description of the continuum is recovered.

The arrays of the conservative variables U(x, t), convective fluxes f(U) and viscous fluxes g(U) are defined as follows:

U =            ρ ρ u ρ E t            f =            ρ u ρ u u T + p I ρ E t u + p u            g =            0 τ τ • u + q            (2)
where ρ(x, t), u(x, t) and E t (x, t) are the density, the velocity and the total specific energy. p(x, t), τ(x, t) and q(x, t) are the pressure, the viscous stresses tensor and the power exchanged by conduction. The problem is closed when the equation of state for the pressure (Polytropic Ideal Gas model), the constitutive equations for the viscous stresses tensor (Newtonian fluid model) and the power exchanged by conduction (Fourier model) are specified. The source term h(U) allows for modeling gravitational forces, Coriolis acceleration, porous media.

The source term, together with the ALE framework, is also used when dealing with Multiple Reference of Frame (MRF) simulations, such as for multi-stage compressors and turbines [START_REF] Robert | Recent enhancements to the FUN3D flow solver for moving-mesh applications[END_REF]. Given Ω as the rotor row angular velocity, by enforcing the ALE velocity as w = Ω × r and the source term as h(U) = {0, -Ω × ρu, 0} T it is possible to solve the Navier-Stokes equations (1) for rotor cases without actually moving the mesh.

The solver is capable to handle hybrid unstructured meshes, i.e. meshes where different cells could have different number of faces and where there are no simple index-based strategies to identify cell neighbors. This feature offers significant advantages in terms of handling complex geometries, thanks to numerous unstructured mesh generation tools available. However, this comes at the price of a reduced efficiency in memory access patterns, especially for GPU executions. Therefore the code must be designed to maximize coalesced memory access and minimize branch divergence [START_REF] Favale | A gpu parallelized two fields full potential formulation for real gases[END_REF].

Within an unstructured grid framework, let us denote with Γ i = {Γ i1 , Γ i2 , . . . , Γ in } the set of n interfaces delimiting the i-th finite volume ∆V i , the FV cell-centered discretized form of the RANS equations corresponds to the following system of ordinary differential equations:

d(U i ∆V i ) dt = Q i = - n j=1 |Γ i j | (F i j -G i j ) +H i ∆V i ∀ i, (3)
where U i (t) is the vector of unknown variables averaged on the i-th cell ∆V i and collocated in its center. F i j (U i , U j , W i j ) and G i j (U i , U j ) are the convective and viscous numerical fluxes through interface Γ i j . Finally H i (U i ) is the source term evaluated in i-th cell center.

Concerning the space discretization operators, in order to obtain a monotone but sharp solution near shock waves and contemporarily to achieve a second order accuracy in space in smooth flow regions, a combination of the monotone first order accurate Roe's upwind flux and a second order accurate centered flux is used. The blending of the two is automatically controlled by means of the flux limiter by van Leer, requiring the numerical solution on the extended cells i * and j * to be available [START_REF] Romanelli | A Free Approach to Modern Computational Aeroelasticity[END_REF]. This operation is performed using a connectivity data structure built in the pre-processing stage. As alternative options AUSM+ [START_REF] Liou | A sequel to ausm: Ausm+[END_REF] and CUSP [START_REF] Tatsumi | Flux-limited schemes for the compressible navier-stokes equations[END_REF] fluxes are implemented. In this case second order accuracy is achieved by means of min-mod limited solution reconstruction.

The viscous numerical fluxes are assembled using a celllimited Gauss algorithm for the evaluation of the velocity and temperature gradients. As alternative option a thin layer approximation is also implemented, which is much more efficient since it does not require to assembly the gradient fields but has a limited validity interval. Turbulence is modeled using a 1-equation Spalart-Allmaras model in which the additional transport-diffusion-reaction equation is solved for the new variable ν [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF]. Automatic wall treatment is implemented to handle meshes at all nondimensional wall distance (y + ), raging from the viscous sublayer to the log layer, thanks to a blended formulation [START_REF] Popovac | Compound wall treatment for RANS computation of complex turbulent flows/heat transfer[END_REF][START_REF] Kalitzin | Near-wall behavior of RANS turbulence models and implications for wall functions[END_REF] between the solutions provided by the two layers. Alternatively, a simple algebraic Mixing Length (ML) model is also implemented, where the kinematic viscosity ν t = l 2 B is a function of the strain rate tensor magnitude B and the length l = k min(y, ∆) [START_REF] Chima | An algebraic turbulence model for three-dimensional viscous flows[END_REF].

Besides the usual inflow, outflow and wall boundary conditions other types required for turbomachinery applications are implemented [START_REF] Romanelli | Implementation of explicit density-based unstructured cfd solver for turbomachinery applications on graphical processing units[END_REF], such as Mixing Plane (MP) [START_REF] Hanimann | Development of a novel mixing plane interface using a fully implicit averaging for stage analysis[END_REF][START_REF] Dzung | Konsistente mittelwerte in der theorie der turbomaschinen fuer kompressible medien[END_REF], AMI periodic boundary conditions, total and mass-flow inlet boundary conditions.

Since the explicit nature of the solver, convergence acceleration techniques like Multi-Grid (MG), Residual Smoothing (RS) and Local Time Stepping (LTS) [START_REF] Brandvik | An accelerated 3d navier-stokes solver for flows in turbomachines[END_REF][START_REF] Klostermeier | Investigation into the capability of large eddy simulation for turbomachinery design[END_REF][START_REF] Romanelli | Implementation of explicit density-based unstructured cfd solver for turbomachinery applications on graphical processing units[END_REF] must be used to damp residuals and to obtain convergence rates that can be comparable with the ones provided by implicit solvers. LTS allows different cells to advance with different time-steps based on their geometrical and flow properties, accelerating in particular the convergence of the bigger cells far from the wall. A simplified MG strategy is implemented, an adaption to unstructured polyhedral grids of the formulation originally proposed by Denton and better explained in [START_REF] Romanelli | Implementation of explicit density-based unstructured cfd solver for turbomachinery applications on graphical processing units[END_REF]. Basically the original mesh is agglomerated in multiple levels of coarser meshes. Transients are dispersed more quickly on the coarser levels, as larger time steps are allowed, while retaining the spacial accuracy of the finest [START_REF] Brandvik | An accelerated 3d navier-stokes solver for flows in turbomachines[END_REF][START_REF] Klostermeier | Investigation into the capability of large eddy simulation for turbomachinery design[END_REF][START_REF] Romanelli | Implementation of explicit density-based unstructured cfd solver for turbomachinery applications on graphical processing units[END_REF]. These techniques accelerate the convergence of steady solutions, but can't be directly used for unsteady simulations since the time (often called "pseudo-time" at this point) has no longer any physical meaning. The problem of performing unsteady simulations by maintaining the convergence boost provided by these techniques is addressed using the Dual Time Stepping (DTS) formulation [START_REF] Blazek | Computational Fluid Dynamics: Principles and Applications[END_REF].

FSI numerical formulations

The following subsections are aimed to introduce the numerical formulations that, alongside with the previously introduced CFD strategies, allow the solver to perform FSI simulations. The reader is also referred to [START_REF] Romanelli | Computational Aeroservoelasticity of Free-Flying Deformable Aircraft[END_REF][START_REF] Romanelli | Coupled csd/cfd non-linear aeroelastic trim of free-flying flexible aircraft[END_REF] for a more detailed description of the implemented strategies.

Modal reduction of the structural model

With the aim of obtaining a solver capable of flutter and aerodynamic damping computations, the structural model is reduced using the modal approach. Neglecting the structural damping, the free response equation of the structural system can be represented in the matrix notation as follows:

M { üs (t)} + K {u s (t)} = {0} (4)
where {u s (t)} is the vector of nodal displacements, M is the mass matrix in structural nodal coordinates, K is the stiffness matrix in structural nodal coordinates. From this equation we can obtain modal shapes and frequencies, basically eigenvectors and eigenvalues of the free deformable system.

We can represent the transformation from modal coordinates {q(t)} to nodal coordinates as follows:

{u s (t)} = [U] {q(t)} (5)
Matrix [U] represents modal shapes and its transposed is used to obtain modal displacements from structural nodal displacements:

{q(t)} = [U] T {u s (t)} (6)
This way, using the transformation matrix [U] it is possible to rewrite equation 4 in modal coordinates:

[M] { q(t)} + [K] {q(t)} = {0} (7)
where [M] is the generalized mass matrix and [K] is the generalized stiffness matrix, the matrices in modal coordinates.

Modes can be normalized in different ways, e.g. to obtain a unitary generalized mass, in which case the generalized mass matrix becomes the identity. The main advantage of the modal representation of the structural behavior consists in the possibility to build an accurate and efficient reduced model using just few modes. An important aspect is that when computing modes of rotating components like rotor blades, centrifugal effects must be taken into account. Basically those effects are translated into stiffness contributions by the FEM solver. When the structural system is loaded with aerodynamic forces, equation 7 can be rewritten adding the generalized aerodynamic forces vector to the right hand side:

[M] { q(t)} + [K] {q(t)} = {Q(t)} (8)
Matrix [U] T can be used to obtain generalized aerodynamic forces from aerodynamic forces known on structural nodes:

{Q(t)} = [U] T { f s (t)} (9)
Thus, thanks to relations 5, 9 and equation 8, knowing aerodynamic forces on structural nodes allows the solver to compute displacements on structural nodes.

Aeroelastic Interface

Here we briefly describe the adopted aeroelastic interface. For a more detailed description, the reader is referred to [START_REF] Romanelli | Coupled csd/cfd non-linear aeroelastic trim of free-flying flexible aircraft[END_REF]. Usually aerodynamic and structural grids are independent. This means that aerodynamic nodes {x a } are usually different from structural nodes {x s }. Thus we need a way to project aerodynamic loads from aerodynamic nodes ({ f a (t)}) to structural nodes ({ f s (t)}) in order to compute structural displacements, and a way to project structural displacements from structural nodes {u s (t)} to aerodynamic nodes {u a (t)} in order to provide the FSI solver the updated geometry from which re-compute the flow.

The aeroelastic interface can be represented by means of a linear operator [I]:

{u a (t)} = [I] {u s (t)} (10)
As explained in [START_REF] Romanelli | Coupled csd/cfd non-linear aeroelastic trim of free-flying flexible aircraft[END_REF], it is fundamental for the aeroelastic interface the conservation of momentum and energy exchanged between the aerodynamic and the structural subsystems. This means that the virtual work made by the aerodynamic forces { f a (t)} for the structural displacements interpolated on the aerodynamic nodes {u a (t)} must be equivalent to the virtual work made by the aerodynamic forces interpolated on structural nodes { f s (t)} for the structural displacements {u s (t)}.

The consequence of this requirement is that the transpose of the aeroelastic interface operator can be used to obtain the aerodynamic forces on structural nodes from aerodynamic forces on aerodynamic nodes:

{ f s (t)} = [I] T { f a (t)} (11)
The interface matrix can be computed using different strategies.

The approach here used consists of an interpolation scheme based on a Moving Least Squares (MLS) technique and the use of Radial Basis Functions (RBF). When using RBF, the influence between the structural and aerodynamic nodes is weighted using functions that depends only by their distance:

φ(x s , x a ) = φ(||x s -x a ||) (12)
where x s is a structural point and x a is an aerodynamic point. Different kind of functions can be used, allowing to adjust the smoothness of the interpolation.

It must be noted that the aeroelastic interface binds only the aerodynamic and structural nodes that reside on the surface of the deformable geometry. Thus from the displacements of the structural nodes it is possible to compute only the displacements of the aerodynamic nodes on the moving walls. Furthermore, the aeroelastic interface is computed only once and can be used for all the subsequent FSI simulations like trim, aerodynamic damping and flutter analysis.

The aeroelastic interface matrix is stored and reused each time the aerodynamic loads and the structural displacements have to be exchanged between the aerodynamic and structural meshes. Storing this matrix on the GPU memory would lead to an unacceptable memory overhead. Thus, the interface matrix is stored in the system RAM and computations are performed by the CPU. Although the GPU execution is temporarily stopped to wait the CPU in this phase, the procedure is still computationally efficient. In fact the data exchanged between the CPU and the GPU is limited to the aerodynamic loads and structural displacements of wall nodes. Furthermore the computational time required by this procedure is basically negligible with respect to the time required by the GPU for the convergence of the aerodynamic solution over the updated aerodynamic mesh geometry, thanks also to the modal representation of the structural model that reduces the total number of structural d.o.f. This is due to the fact that few thousands aerodynamic iterations are needed to reach aerodynamic residuals convergence between two geometry updates. This is true both for steady trim analysis and unsteady analysis using the DTS technique.

Internal nodes displacements

Once the aerodynamic wall nodes displacements are obtained, we need to compute the displacements of the aerodynamic mesh internal nodes of the fluid domain. This problem can be addressed using different strategies such as laplacian smoothing. This phase during a typical FSI simulation could easily need to be repeated hundreds or thousands of times, thus a computational efficient algorithm is required. It must be also kept in mind the reduced amount of memory available on GPUs: an explicit algorithm is preferred over a strategy that involves the solution of a system of equations. In particular an Inverse Distance Weighting (IDW) algorithm is adopted for this purpose.

Basically, the displacements u k a of an internal aerodynamic node k are computed weighting the displacements of wall nodes u i a through a function of the inverse of the distance between node k and wall node i:

u k a = N w all i=0 W k-i u i a N w all i=0 W k-i (13)
W k-i is the weighting function and can be expressed as follows:

W k-i = 1 |x k -x i | p (14)
The exponent n can be changed to adjust the smoothness of the results. Usually p = 2 or p = 3 provides good results. IDW weights are constants during the FSI simulation, so they could be stored in a matrix that can be accessed every time the aerodynamic mesh has to be updated. However, storing this matrix would require a large memory overhead. Since the limited amount of GPU memory, IDW weights are recomputed every time they are needed by GPU cores instead of being stored. This would lead to large computational overhead in a typical CPU architecture composed by few cores. However this is not a problem in a typical GPU architecture where several floating point computations can be performed while waiting for data to be recovered from GPU memory. Furthermore, the computational time spent to update the mesh with this algorithm is basically negligible with respect to the time required by the explicit aerodynamic iterations to reach convergence over the updated mesh. This is true both for trim simulations and FSI simulation where the mesh deformation is enforced (aerodynamic damping and flutter analysis).

After this phase, the aerodynamic mesh is changed. The connectivity between elements is preserved, however mesh metrics have to be recomputed from the new positions of the aerodynamic points. This can be efficiently done by the GPU since recomputing the metrics of millions of cells can be concurrently done by the thousands of available GPU cores.

Trim analysis

The aim of this section is to briefly describe the strategy used to perform trim analyses. The FSI algorithms and strategies previously introduced are combined together to perform this numerical simulation. The main differences between a trim analysis and an aerodynamic damping or flutter analysis are basically two. First of all, a trim analysis is basically a steady analysis while aerodynamic damping and flutter analyses require unsteady solutions. Furthermore, in aerodynamic damping and flutter analyses the structural displacements given by modal shapes are used to enforce the deformation on the structure. Conversely, in trim analysis, aerodynamic loads are responsible for the structural deformation, thus the steady solution of system 8 is required. Moreover, trim analysis is performed using an iterative procedure were aerodynamic loads are used to deform the structure. The new obtained shape is then used to recompute the aerodynamic solution. The procedure is repeated until modal and aerodynamic residuals reach convergence.

Before starting the trim analysis two pre-processing steps are required. In particular the modal representation of the structural model of the deformable geometry is needed. This is done only once before any FSI simulation, i.e. trim, aerodynamic damping or flutter analysis. Modal shapes alongside with generalized masses, stiffness and frequencies can be obtained using several available commercial and free FEM solvers like NASTRAN or Code_Aster [START_REF]Code aster website[END_REF]. Different kind of FEM elements can be used, e.g. beam, shell/plate, solid elements, in a trade-off between accuracy and computational efficiency. It is worth to remind that since we are studying rotating blades, centrifugal effects must be taken into account in the modal analysis. FEM solvers can easily account for centrifugal effects as stiffness contributions. Another fundamental step is performing the computation of the aeroelastic interface matrix. When dealing with several thousands of aerodynamic and structural wall nodes, the assembly of this matrix can be a computationally intensive job, requiring minutes on a modern multi-core CPU. However, it is needed to compute the aeroelastic interface matrix only once, as it can be saved on disk and re-loaded every successive FSI simulation.

After the modal representation of the structural behavior and the aeroelastic interface matrix computation, it is possible to start the trim analysis for the entire rotor characteristic curve. First of all the pure aerodynamic solution at choking condition is computed, providing a guess solution to start the entire trim analysis. This is done using the explicit timestepping procedure with the help of convergence acceleration techniques since we are searching a steady solution. At this point, aerodynamic loads over the aerodynamic mesh wall boundary points { f a } are obtained. Using relation 11, thus using the aeroelastic interface matrix, aerodynamic loads are projected over the structural mesh wall points of the deformable blade in order to obtain { f s }. Since the structural model is represented by modal properties, generalized aerodynamic forces have to be computed. Relation 9 allows the solver to obtain forces in modal coordinates {Q} from forces in structural nodal coordinates { f s }. Now it is possible to compute the response, in term of modal displacements {q}, of the structure due to aerodynamic loads. Basically we need the steady version of the forced system 8:

[K] {q} = {Q} (15)
Since usually just few modes are sufficient for an accurate representation of the structural behavior, the solution of system 15 requires very small computational times. Once modal displacements {q} are obtained, in the next step the structural nodal displacements {u s } are computed using relation 5. Using again the aeroelastic interface matrix and relation 10, aerodynamic mesh wall nodes displacements {u a } are recovered from structural mesh wall nodes displacements {u s }. Now it is possible to compute the aerodynamic mesh internal nodes displacements using the IDW strategy and recompute the aerodynamic mesh metrics. With the updated aerodynamic mesh a new steady solution is computed. Convergence is reached when aerodynamic residuals, generalized aerodynamic loads and modal displacements residuals are under a user-specified tolerance.

For the Rotor 67 test case about 500 explicit aerodynamic iterations are performed between two mesh updates and about 15-20 mesh updates are required for the aeroelastic convergence of one characteristic curve point. It must be underlined that thanks to the efficient IDW algorithm and the reduced amount of data transferred between the CPU and the GPU at each mesh update, the computational time required by the mesh update is basically negligible with respect to the time required by explicit aerodynamic iterations.

Hardware and software aspects

Since the solver is written in OpenCL, its source code is composed by two sets of files: one for the "host" and one for the "device". The host is basically the CPU that starts the solver and enqueues work for the device. The host code is written in standard C/C++ language. The device is the GPU or the CPU that is actually used to perform the CFD/FSI computations. The device code is composed by functions called "kernels" written in OpenCL C, a language similar to the standard C99 but with some restrictions. Functions and types of the OpenCL API are used in the host source code in order to organize the work (kernels executions) that is sent to the device and to exchange "buffers" (basically memory arrays) between the host and the device. Most CPU and GPU vendors, e.g Intel, AMD, NVIDIA, have their own OpenCL-compliant implementations, allowing the solver to be executed natively on a wide range of devices [START_REF] Mangani | Comparison of acceleration techniques on cfd open-source software for aerospace applications[END_REF][START_REF] Romanelli | Implementation of explicit density-based unstructured cfd solver for turbomachinery applications on graphical processing units[END_REF]. During kernel executions, the whole work is split in chunks, called "work-items", that are executed in parallel by GPU or CPU cores. Since nowadays GPUs have thousands of cores, their architectures are well suited for the data parallelism where the same operation have to be performed on a large dataset. This is basically what is needed for an explicit solver. It must be noted however, that if different "work-items" have to perform different computations (branch divergence), a performance loss is inevitable with GPUs. Thus, device code must be opportunely optimized for the GPU architecture. Another typical aspect of GPUs is the memory optimized for sequential access. Again, kernels must be optimized to reduce performance loss given by a typical hybrid unstructured mesh addressing. CPUs, nonetheless, are less susceptible to this kind of problems thanks to their MIMD (Multiple Instructions Multiple Data) architectures and bigger caches. Furthermore, GPUs main limitation is the reduced amount of "global memory", basically the amount of memory available to store buffers. A continuous data flow between the CPU and the GPU during the simulation would not be a good strategy since the PCI-Express bus that connects them would be a bottleneck. With these restrictions in mind, the explicit formulation is the better choice since the matrix storage is avoided. Despite these types of GPU problems and limitations, a typical 400 USD gaming GPUs (e.g AMD 290X) can provide over 5 TFLOPS of single precision performance, an order of magnitude more than a CPU of the same price level (e.g. Intel i7 3930K). Since the main advantages in using GPUs to perform numerical computations can be achieved with single precision, all the equations in the solver are opportunely kept dimensionless to reduce numerical problems. Finally it is important to notice that with the ability to perform computations in single precision floating point format and having low memory requirements, the solver is compatible with both HPC (High Performance Computing) GPUs and cheap gaming GPUs. HPC GPUs typically are an order of magnitude more expensive than gaming GPUs but exhibit more memory, ECC memory feature and higher double precision performance. However, HPC GPUs usually have about the same single precision peak computational power as gaming GPUs. All the simulations performed in this article were carried out using cheap gaming GPUs and single precision. Results were also reproduced with double precision and using an ECC memory-compliant GPU (Tesla C1060) but appreciable differences were not found [START_REF] Romanelli | Implementation of explicit density-based unstructured cfd solver for turbomachinery applications on graphical processing units[END_REF].

The solver uses the OpenFOAM framework to carry out pre-processing and post-processing phases. In particular mesh loading, metrics computations, domain decomposition for multi-GPU executions, AMI (Arbitrary Mesh Interface) weights and addressing computations, and monitoring relevant performance metrics (i.e. mass-flow, pressure ratio) are computed with the help of the OpenFOAM API. The device code that performs the actual CFD/FSI computation is however written from scratch with OpenCL API and OpenCL C. Thanks to the OpenFOAM framework and an the opportunely tuned GPU code, the solver is compatible and computationally efficient with hybrid unstructured meshes [START_REF] Romanelli | Implementation of explicit density-based unstructured cfd solver for turbomachinery applications on graphical processing units[END_REF].

RESULTS AND DISCUSSION

NASA's Rotor 67 trim

Here we describe the results of the trim analysis of the well known NASA's Rotor 67 case. Rotor 67 represents a typical fan blade geometry, here modeled as an isolated rotor, for which two sets of experimental data are available (one taken in 1984 and one lately in 2004). For this transonic axial fan test case, a number of CFD performance predictions have been published [START_REF] Arnone | Viscous analysis of three-dimensional rotor flow using a multigrid method[END_REF][START_REF] Grosvenor | Rans prediction of transonic compressive rotor performance near stall[END_REF]. The fan angular velocity adopted for the tests is 16043 rpm.

For what concerns the aerodynamic mesh, the details for the hub and the blade are shown in figure 1. The adopted unstructured grid is composed by about 1.1 • 10 6 hexahedral cells and wall regions are discretized with an average y + 0.5. Respectively, total and average pressure boundary conditions are imposed at inlet and outlet. The adopted turbulence model for aerodynamic and aeroelastic simulations over the entire characteristic curves is Spalart-Allmaras. However, a comparison between Spalart-Allmaras and Mixing Length turbulence models is also provided to check the mass flow convergence at choking conditions. The structural mesh is instead composed by about 8 • 10 4 tetrahedrons (solid elements). Only the blade was modeled with FEM elements while the hub is supposed to be non deformable. The modal analysis was performed with the Code_Aster FEM solver. Just 3 modes are sufficient for an accurate representation for trim purposes. However, 4 modes are used in this simulation to check modal convergence. Using 5 or more modes leads basically to negligible differences in the results. Table 1 shows modes frequencies and describes their shapes. With the aim of the characteristic curves computation, let us first discuss the choking point. The solver needs approximately 50 • 10 3 explicit iteration to reach pure aerodynamic steady solution for this particular point. Using an AMD 290X GPU the required time is about 25 min. Figures 2 show the residual and mass flow values during convergence at choking point using Mixing Length and Spalart-Allmaras turbulence models. From this solution, the trim analysis is started and another 10min of wall time are required by the same GPU to reach aeroelastic convergence. Experimental value of the mass flow is 1.589 kg/s. The pure aerodynamic solution and the aeroelastic solution converges to a mass flow of 1.570 kg/s with an error of 1.76%. In figures 5 and 6 the efficiency and pressure ratio characteristic curves as a function of non-dimensional mass flow are compared with the two sets of experimental results. It is possible to see that both the aerodynamic and aeroelastic solutions are in good agreement with the experimental points. It must be stressed however the measurement uncertainty of the 1989 data set. It is possible to see that the trim solution basically provides the same results as the aerodynamic solution for the efficiency curve, and a slightly higher curve for the pressure ratio, especially near the efficiency peak region. Generally, however, the differences between the aerodynamic and the aeroelastic solution are small, due to the high stiffness of the blade. Figure 7 shows the displacements ({q}) of the four considered modes during convergence over the trimmed characteristic curve. It is possible to see that the curves are made by multiple steps, representing the convergence over different points of the characteristic curves, starting from the choking point and moving to the stall region. From modal displacements and modal stiffness it is possible to compute the elastic energy contribution of each mode. Basically the main contribution to the total elastic energy is given by the first mode, while the second, the third and the fourth mode contributes are at least two order of magnitude lower. Using just 3 modes basically provides the same results. In particular, we can see that mode 4 displacements are basically 0 over the entire characteristic simulation. Using 5 or more modes with higher stiffness leads to further negligible contributions to the results. To give a more complete view to the trim analysis using modal representation of the structural behavior, figure 8 shows generalized aerodynamic forces ({Q}) of the four considered modes.

Figures 9 and10 shows the detail of efficiency and pressure ratio characteristic curves for the Rotor 67 near the peak efficiency point. It is possibile to see that when using different number of modes the curves are slightly different. However, using 3 or 4 modes the results are basically identical.

Figures 11 and12 show the relative Mach number contours predicted at the peak-efficiency point. A comparison with experimental data is also provided at 10% and 30% span from the shroud wall. It is possible to see that both the pure aerodynamic and the aeroelastic solutions are in agreement with the measurements, in particular for what concerns the shock position in figure 11. The differences in shocks locations between the aerodynamic and trim solutions are very small, basically negligible.

In order to highlight the importance of a trim analysis during the design and optimization phases of a rotor blade let us try to see the effects given by changing the blade stiffness. Since the limited amount of space available we will not perform a detailed quantitative analysis but just a qualitative study of the changing in the blade tip displacements obtained through a stiffness reduction of the blade. As an example we can reduce the first mode (flexural shape) stiffness of about ten times and see the blade shape at choking conditions. Figure 13 shows the detail of the tip displacements. It is easy to see the differences from figure 4: the blade tip is now characterized by bigger tip displacements. It is straightforward that one of the possible aims of a rotor trim analysis would be an accurate study of the interaction between the blade tip and the shroud wall.

GPU acceleration

Here we briefly present the GPU speedup results for the Rotor 67 test case. For more detailed benchmarks regarding the GPU solver, the reader is referred to [START_REF] Romanelli | Implementation of explicit density-based unstructured cfd solver for turbomachinery applications on graphical processing units[END_REF][START_REF] Mangani | Comparison of acceleration techniques on cfd open-source software for aerospace applications[END_REF]. All the convergence acceleration techniques (MG, RS, LTS) were active during these benchmarks in order to honestly represent a true simulation. The adopted cost metrics is the elapsed time per iteration per cell and the comparison involves CPUs in single thread execution, CPUs with multi-threading and GPUs. Thanks to OpenCL, the solver can be natively executed by several different CPUs and GPUs architectures without modifications. Benchmarks were carried out with single precision to fully exploit peak computational power provided by the adopted devices. Furthermore, when available, the compiler autovectorization is used to exploit SSE/AVX capabilities of the modern multi-core CPUs. Table 2 shows the CPUs used for the comparison, basically an entry level CPU (AMD Phenom II X4 840) and an high end CPU (Intel i7 3930k). Table 3 shows instead the GPUs used for the comparison, raging from an entry level gaming GPU (NVIDIA GTX 650) to an high end gaming GPU (AMD 290X). It is worth remembering that gaming GPUs are here preferred thanks to their high single precision performances provided at relatively low cost with respect to typical HPC GPUs like NVIDIA Tesla or AMD FirePro series.

Tables 2 and3 also show prices and theoretical singleprecision peak performance of the adopted devices for a more honest comparison. For the AMD CPU, however, data regarding the theoretical peak performances where not found on the official website (although it can be estimated to be about 30-50 GFLOPS).

Table 4 shows the results for the time per iteration per cells for each available device. From these results it is possible to build table 5. In this table, at each row, CPUs in single and multi-thread execution timings are taken as reference for the speedup value computation. As an example, on the last row, we can see that the AMD 290X is 13.06 times faster than the Intel i7 3930k CPU using all the 6 cores. Considering that this GPU costs about 400 USD and the CPU costs about 500 USD, it is clear the advantage to perform a simulation using the GPU. Even on a possible low-end workstation composed by the AMD Phenom II X4 840 CPU and the NVIDIA GTX 650 GPU we can see that using the GPU allows the solver to complete simulations about 16 times faster. Furthermore it is possible to see that the cheapest GPU (GTX 650) is about 3.5 times faster than the most expensive CPU (i7 3930k) used for these test, but costs only 1/5 of the price. Alongside with the price, it is possible to consider the TDP (Thermal Design Power), that gives an idea of the power consumption. For instance, the Intel CPU and the NVIDIA GTX 660 GPU have the same TDP but the GPU is about 6 times faster than the CPU. This also underlines the power consumption advantages of using GPUs for numerical simulations when the particular problem allows an efficient execution on these kinds of devices, a non secondary aspect in clusters and supercomputers.

CONCLUSIONS

In this paper we have described the architecture and implementation of a GPU-accelerated explicit density-based aeroelastic RANS solver for turbomachinery applications. The solver is based on the OpenFOAM framework but the code that actually performs CFD/FSI computations is written from scratch with OpenCL API and OpenCL C in order to be compatible with a wide range of multi-core CPUs and GPUs. The solver was validated using a well known turbomachinery case, the NASA's Rotor 67 axial fan. For this particular case two experimental data sets are available. A comparison between the pure aerodynamic solution and the trim solution for the entire characteristic curve was provided in order to asses the possible advantages given by simulating the interaction between the flow and the blade structure. However, as the results suggest, for this particular case the differences between the aerodynamic and aeroelastic solution are very small. This is basically due to the fact that with respect to other typical aeronautical cases (e.g. helicopter blades and planes) these fan blades are very stiff. Thus, blades displacements are very small in all the considered operating conditions and their effects on the flow are basically negligible.

Alongside with the aeroelastic investigation, a comparison between the time per iteration per cell given by different CPUs and GPUs was provided, highlighting the advantages given by GPU executions in term of wall time and costs.

Figure 1 .

 1 Figure 1. Detail of the computational mesh for Rotor 67 test case.

Figure 2 .

 2 Figure 2. Residuals and mass-flow convergence history for Rotor 67 test case. Comparison between Mixing Legnth (ML), Spalart-Allmaras (SA) turbulence models and experimental data (EXP).

Figure 3

 3 Figure 3 and 4 shows the displacement field and a detail of the blade edge of the deformed configuration obtained with the trim solution at the choking point. It is possible to see from the magnitude (colors from blue to red) that the deformation of the blade is mainly flexural, basically the shape of the first mode, with smaller influence from higher modes.

Figure 3 .

 3 Figure 3. Blade displacements field at choking conditions, colors from blue to red represent displacements magnitude.

Figure 4 .

 4 Figure 4. Detail of blade deformation. Grey surface: original non deformed blade. Surface with edges: deformed blade at choking.

Figure 5 .

 5 Figure 5. Efficiency as a function of mass flow for the Rotor 67 test case.

Figure 6 .

 6 Figure 6. Pressure ratio as a function of mass flow for the Rotor 67 test case.

Figure 7 .

 7 Figure 7. Modal displacements of the considered modes for the Rotor 67 trim analysis

Figure 8 .

 8 Figure 8. Generalized aerodynamic forces of the considered modes for the Rotor 67 trim analysis.

Figure 10 .

 10 Figure 10. Pressure ratio as a function of mass flow for the Rotor 67 test case, detail near peak efficiency using different number of modes.

Figure 11 .

 11 Figure 11. Comparison of relative Mach number contours with experiments at peak efficiency, 10% span from shroud. From left to right: numerical trim solution, experimental data, numerical aerodynamic solution.

Figure 12 .

 12 Figure 12. Comparison of relative Mach number contours predicted with experiments at peak efficiency, 30% span from shroud. From left to right: numerical trim solution, experimental data, numerical aerodynamic solution.

Figure 13 .

 13 Figure 13. Mesh deformation at choking point using 1/10 the stiffness of the first mode.

Table 1 .

 1 Rotor 67 modes frequencies and description.

	Mode Frequency (Hz) Description
	1	760.3	1st flexural
	2	2174	2nd flexural
	3	3146	1st torsional
	4	4920	Mixed

Table 2 .

 2 CPU used for the benchmark.

	Vendor	Model	Cores SP Performance Frequency	Price	TDP	OpenCL runtime
	AMD Phenom II X4 840	4	N/A	3.2GHz	∼ 100 USD 95 W	AMD APP SDK
	Intel	i7 3930k	6	∼ 170 GFLOPS	3.2GHz	∼ 500 USD 130 W Intel OpenCL Runtime

Table 3 .

 3 GPU used for the benchmark.

	Vendor	Model	Global Memory Cores SP Performance Frequency	Price	TDP	OpenCL runtime
	NVIDIA GTX 650	1 GB	384	∼ 800 GFLOPS	928 MHz ∼ 100 USD 110 W	CUDA SDK
	NVIDIA GTX 660	2 GB	960	∼ 1800 GFLOPS 980 MHz ∼ 200 USD 130 W	CUDA SDK
	AMD	290X	4 GB	2816 ∼ 5600 GFLOPS 1040 MHz ∼ 400 USD 290 W AMD APP SDK

Table 4 .

 4 Time per iteration per cell.

	Phenom 1-core Phenom 4-core i7 1-core i7 6-core GTX 650 GTX 660	290X
	5.99e-6	1.74e-6	1.82e-6	3.67e-7	1.04e-7	6.37e-8	2.81e-8

Table 5 .

 5 Relative speedups.

	Phenom 1-thread Phenom 4-threads i7 1-thread i7 6-threads GTX 650 GTX 660	290X
	1.00x	3.44x	3.29x	16.32x	57.6x	94.04x	213.17x
	0.29x	1.00x	0.96x	4.74x	16.26x	27.32x	61.92x
	0.30x	1.04x	1.00x	4.96x	17.5x	28.57x	64.77x
	0.06x	0.21x	0.20x	1.00x	3.53x	5.76x	13.06x


Yiwei Wang

Xiaocui Wu

Chenguang Huang
email: huangcg@imech.ac.cn

S Y Mpo

O Tatin

Ventilated partial cavitating flow around a blunt body near the free surface

Keywords:

Ventilated cavitation occurs as an extremely complicated problem if the free surface is close to the cavity boundary around the high-speed underwater vehicle. The cavitating flow around a blunt axisymmetric body very near the free surface is investigated in the present paper. The typical experiment is carried out in a launching system based on an SHPB device, and the numerical scheme is established on the basis of large eddy simulation (LES) and volume of fraction (VOF) methods. Unsteady behaviors such as air injection and shedding of the cloud cavity were obtained, and good agreement was achieved between the numerical and experimental results. The distinctions of evolution features between the cavities on the up and down sides are presented and analyzed. Firstly, strong injection of non-condensable air occurs and changes the media property inside the cavity, which makes the cavity much larger and more stable. Furthermore, because of the small distance between the vehicle and the free surface, the re-entry jet generated in the upper part is very thin and cannot cut off the main cavity clearly, which induces that the upper part of the cavity is approximately unchanged after the growth stage and no shedding phenomenon occurs inside it.

INTRODUCTION

The interaction between the cavitation region and the free surface is an important issue for the fast cruise of surface vehicles. Ventilated cavitation occurs as an extremely complicated problem if the free surface is close to the vehicle or the cavity boundary as Faltinsen mentioned [START_REF] Faltinsen | The effect of gravity and cavitation on a hydrofoil near the free surface[END_REF]. The gas above the free surface can inject into the cavity, which may change the dynamic behaviors of the cavity significantly. The effect of wave elevation is also coupled with the distribution of cavitation region.

Relevant research works on the interactions between the free surface and the cavitating flow is very limited in literature. Most of these works were focused on the stable cavity, and the interfaces of the cavity and the free surface are separated between each other. Theoretical and numerical methods were established, and the influences of the free surface on the cavitating flow around hydrofoils and axisymmetric bodies are obtained [START_REF] Faltinsen | The effect of gravity and cavitation on a hydrofoil near the free surface[END_REF][START_REF] Bal | A bem for the prediction of free surface effects on cavitating hydrofoils[END_REF][START_REF]High-speed submerged and surface piercing cavitating hydrofoils, including tandem case[END_REF][START_REF] Amromin | Analysis of body supercavitation in shallow water[END_REF]. On the other side, neglecting the influence of the free surface, the dynamic characteristics of the cloud cavitating flow are often linked with the motion of vortexes. Based on the experimental and numerical approaches, the evolution of cavities including growth, re-entry, shedding and collapse were gained [START_REF] Stutz | X-ray measurements within unsteady cavitation[END_REF][START_REF] Ji | Large eddy simulation and theoretical investigations of the transient cavitating vortical flow structure around a naca66 hydrofoil[END_REF][START_REF] Roohi | Numerical simulation of cavitation around a two-dimensional hydrofoil using vof method and les turbulence model[END_REF][START_REF] Yu | Study of characteristics of cloud cavity around axisymmetric projectile by large eddy simulation[END_REF].

The cavitating flow around a blunt axisymmetric body very near the free surface is investigated in the present paper. The typical experiment is carried out in a launching system based on an SHPB device, and the numerical scheme is established on the basis of large eddy simulation (LES) and volume of fraction (VOF) methods. The distinctions of evo-lution features between the cavities on the up and down sides are presented. The effects of the free surface and the gas injection on the cavity shape are analyzed.

EXPERIMENT SETUP

The experiment was performed by using a launching system (as shown in Figure 1) on the basis of the SHPB technology [START_REF] Yu | Study of characteristics of cloud cavity around axisymmetric projectile by large eddy simulation[END_REF][START_REF] Wei | A scaled underwater launch system accomplished by stress wave propagation technique[END_REF], which can transiently accelerate the projectile with slight disturbance on the water. The projectile is driven by a strong stress wave. The acceleration process is very short which is about 200 s for 20m/s. Therefore the acceleration cannot be captured by using the high-speed camera. Moreover, because solid steel projectile is used which is rather heavy, the deceleration is also very small in the observation period. So the speed of projectile is considered approximately uniform for numerical simulation.The projectile used in this study is a blunt axisymmetric body. The total length is 150 mm, and the diameter is 37 mm. Photographs of typical cavitation can be obtained using a high-speed camera with 12000 fps. For the typical condition investigated in the present paper, the distance between the free surface and the up side of the projectile is 5mm, and the analysis of obtained images indicates that the speed is approximately uniform at 19.1 m/s. The cavitation number calculated as 0.537.

The cavity is nonaxisymmetrically affected by the free surface. Thus, the length on the upper and lower sides of the cavity is measured, as shown in Figure 2. The precision of the length and thickness is approximately a pixel of the image, which stands for about 0.65 mm.

NUMERICAL METHODS

Governing equations

To simulate the motions of liquid water and vapour and air, including the phase change, mixture/multiphase flow equations are adopted. Continuity and momentum equations for the mixture are established as follows.

∂ ρ ∂t + ∂( ρu j ) ∂ x j = 0 (1)
∂( ρu i ) ∂t + ∂( ρu i u j ) ∂ x j = - ∂p ∂ x j + ∂ ∂ x j ( µ ∂u i ∂ x j ) (2)
where u i is the velocity component in the i direction, and p is the pressure. Laminar viscosity µ is defined as a volumeweighted average of the three components as

µ= (1-α v -α a ) µ l + α v µ v + α a µ a (3)
where α v and α a are the vapour and air volume fractions, respectively. Mixture density ρ is defined as

ρ= (1-α v -α a ) ρ l + α v ρ v + α a ρ a (4)
The volume fractions of vapour and air are governed by the following mass transfer equations:

∂( ρ v α v ) ∂t + ∂( ρ v α v u j ) ∂ x j = ṁ+ -ṁ- (5) ∂( ρ a α a ) ∂t + ∂( ρ a α a u j ) ∂ x j = 0 (6)
where source terms ṁ+ and ṁ-represent the mass transfer rate of evaporation and condensation, which are derived from the bubble dynamics equations of generalised RayleighPlesset equation by Zwart et al. [START_REF] Zwart | A twophase flow model for predicting cavitation dynamics[END_REF] as follows.

ṁ+ = F va p 3a nuc (1 -α v ) ρ v R B √ 2 3 max(p v -p, 0) ρ l (7) ṁ-= F cond 3α v ρ v R B √ 2 3 max(p -p v , 0) ρ l ( 8
)
where generalised bubble radius R B is set at 10 -6 m, nucleation site volume fraction a nuc is set at 5 × 10 -4 , the evaporation coefficient is set at 50 and the condensation coefficient is set at 0.01.

Simulation procedure

The governing equations are solved by a LES approach based on Smagorinsky-Lilly model. Unsteady numerical simulations are performed on the basis of finite volume method with coupled scheme by using the commercial CFD software ANSYS-FLUENT. The equations are discretised by a secondorder implicit scheme in time and a bounded central differencing scheme in space. The pressure staggering option was selected for pressure interpolation with the modified highresolution interface-capturing scheme used for the volume fraction. The unsteady cavitating flow simulations were started from a uniform flow field, and the time step was set at 10 -5 s. The computational domain is discretised with a blockstructured grid, which is refined around the model and near the free surface, and only half of the model is considered (as shown in Figure 3). A semi-infinite projectile model was used, and the effect of the on the shoulder cavity is neglected. The cell number is approximately 4 million with good orthogonality. Independence analysis and validation of similar grid have been performed on a simulation of the typical evolution of cloud cavitating flow [START_REF] Yu | Study of characteristics of cloud cavity around axisymmetric projectile by large eddy simulation[END_REF]. The model was fixed, with the free surface moving towards the model. The inlet velocity is set as 19.1 m/s, with no turbulent perturbations for the velocity inlet boundary condition.

RESULTS AND DISCUSSION

Overall cavities evolutions

From the experimental and numerical results, the quasi-periodic development of cavity shape is obtained. Photographs on typical moments of different stages are as shown in Fig. 4, and the volume fraction of liquid water is 0.9 on the isosurface. The contours of water volume fraction on the symmetric plane and wall surface are also given, in which the red The variations of cavity lengths on the up and down sides are as shown in Figure 5. The length development and variation of the obtained numerical results agree fairly with the experiments, including the breaking phenomenon of the lower side cavity and the relationship between cavities on upper and lower sides. These characteristics can affect the nature of the bubble within the media, the impact of the return jets and other characteristics of the unsteady factors.

Refer to the evolutions aforementioned, we can see the cavity on the upper side stays approximately stable, and its length varies in a very small range after the growth stage. The length variation on the lower side is similar with the case under the same cavitation number without the free surface but has much larger amplitude. These characteristics are relevant to the particularities of the current case, for example the free surface are close to the cavity boundary and the air above the free surface has injected into the cavity. These characteristics can affect the media property inside the cavity, and have impact on the re-entry jet and other instability factors, which will be analyzed and discussed in detail separately in the following contents.

Air entrainment into the cavity

The thickness of the cavity increases gradually in the growth stage, and regions of free air space and the cavity will connect after the cavity boundary intersects with the free surface. Because the pressure above the free surface is about 1atm, which is much higher than the pressure inside the cavity, highspeed air entrainment forms in the connecting channel. At the beginning of the air entrainment, high-speed air flow impacts on the leading edge of the cavity, which raises the upper cavity boundary and thickens the cavity (as shown in Figure 6). At the same time, air also flows to the low pressure region inside the cavity on the lower side and fills it (as shown in Figure 7). The region influenced by the air entrainment is foam-like in the experiment photograph (as shown in Figure 4-1a), the shape of which is similar with the region with large negative Y velocity in numerical results (as shown in Figure 7). Thereafter, the water layer between the cavity and the outside space is affected by gravity, and moves downward to the wall. The reattachment occurs at the time of about 4ms, which ends the air entrainment finally. As shown in Figure 8, the water layer was curled under the joint effects of gravity and air entrainment, and the air entrainment is compressed just inside a very thin gap. The circular velocity inside the cavity also decreases remarkably at this moment (as shown in Figure 9). The whole air entrainment process lasts about 2ms in the cavity growth stage, carrying large amounts of air to fill the cavity. On the one hand, it increases the inner pressure significantly, so the inner pressure is actually remarkably higher than the saturation pressure in the natural cavitation bubbles, which can be considered as that the cavitation number is decreased. On the other hand, the ventilated air cannot condense, which will weaken the condensation and collapse processes on the cavity trailing edge remarkably, and will also increase the cavity length and stability.

Influence on re-entry jets and cavity shedding

Strong adverse pressure gradient forms in the cavity growth stage in the area between the high pressure region at the cavity closure and the low pressure region inside the cavity. Flow detachment will occur, and re-entry jets is generated as a water layer flowing upstream near the wall inside the cavity. The cavity profile and the re-entry jet velocity on the symmetric plane are as shown in Figure 10. It demonstrates that velocity of the re-entry jet on the upper side is high but the water layer is very thin. In contrast, re-entry jet velocity on the lower sides is slower, but the water layer is much thicker (as shown in Figure 10). And the middle and lower parts of the cavity are significantly disturbed by the re-entry jet (as shown in Figure 11). With the further development of the re-entry jet, the cavity on the upper side is cut off at first. However, because re-entry jet is very thin, there is only a small gap between the newly incepted and the remaining cavities (as shown in Figure 12). The gap will be filled soon by the growth of the new cavity. In contrast, the re-entry jet in the lower part of the cavity is thick and cut off the cavity by intersecting with the outside cavity boundary. It is also seen in the experiment photograph (as shown in Figure 4-1b) that the upper part of the cavity is smooth while the middle and lower parts have many small holes on the surface.

When the re-entry jet reaches the leading edge of cavity on the lower side, the remaining bubbles only retain a small space of the origin cavity. A region with high content of liquid water can separate the new cavity form the remaining part, inducing obvious shedding phenomenon. However, only small fluctuation occurs on the cavity boundary on the upper side, and the cavity approximately remains in a constant area.

Aforementioned results indicate that, if the vehicle moves closely to the free surface, the water layer flowing around the upper side is thin and the restriction of water layer to the flow is weak. Thus, the re-entry jet generated does not also have enough strength and thickness to cut off the main cavity clearly, so notable fluctuation of cavity shape cannot be induced on the upper side. However, the influence of the free surface on the cavity on the lower side is relatively small, so powerful re-entry jet can be generated, and unsteady evolutions including breaking and shedding are still exist.

CONCLUSIONS

An experiment on the ventilated cloud cavitating flow around the axisymmetric projectile near the free surface has been presented in this paper. A relevant numerical simulation was performed based on VOF and LES methods. Unsteady behaviors such as air entrainment and shedding of the cloud cavity were obtained, and good agreement was achieved between the numerical and experimental results.

Strong air entrainment occurs when the cavity grows and contacts the free surface. The ventilated non-condensable air makes the cavity much larger and more stable than the cloud cavity under similar cavitation number without the influence of the free surface.

The upper part of the cavity is approximately unchanged after the growth stage and no shedding phenomenon occurs inside it. This is because the restriction of water layer to the flow flowing around the upper side is weak, and the re-entry jet generated in the upper part is very thin and cannot cut off the main cavity clearly.

However, there are still some differences on the phenomena after cavity shedding, and the current numerical simulation method should be further validated by comparing with classic experimental results of ventilation cavitation and also be improved in the future. Through other experimental and numerical results with various submergence, it is found that this specific phenomenon discussed in the present paper just may occur in a very small submergence range. The ventila-tion dose not happen if the submergence is large, while the stable ventilated supercavity appears if the submergence is smaller. The detailed results with other submergence may be provided in other subsequent papers.

Moreover this present work can be extended to other shape vehicles including the hydrofoils, and the cavity may be controlled to remain stable by adjusting the distance between the vehicle and the free surface.

Figure 1 .Figure 2 .

 12 Figure 1. Underwater launch system

Figure 3 .

 3 Figure 3. Mesh near the head of the projectile

Figure 4 .

 4 Figure 4. Time evolution of cavitation patterns (the breaking location is marked by the arrows in the bottom views)

Figure 5 .Figure 6 .

 56 Figure 5. The comparison of the variation of cavity length between numerical and experimental results

Figure 7 .Figure 8 .

 78 Figure 7. Distribution of the air entrainment at the beginning (the flow time is equal to 2ms, translucent surfaces represent the cavity boundary and the free surface, line represents the volume fraction of air, and the color represents the velocity in the Y direction)

Figure 9 .Figure 10 .

 910 Figure 9. Distribution of the air entrainment at the end (the flow time is equal to 4ms, translucent surfaces represent the cavity boundary and the free surface, line represents the volume fraction of air, and the color represents the velocity in the Y direction)

Figure 11 .Figure 12 .

 1112 Figure 11. Shape of the re-entry jet (the flow time is equal to 7ms, translucent surfaces represent the cavity boundary and the free surface, line represents the volume fraction of air, and the color represents the velocity in the X direction

Figure 13 .

 13 Figure 13. Cavity breaking on both sides on the symmetric plane (the flow time is equal to 13ms, and the line represents the interfaces where the volume fraction of the liquid water is equal to 0.5, while the color represents the velocity in the X direction)

ACKNOWLEDGMENTS

This research was sponsored by National Natural Science Foundation of China (Grant numbers 11332011 & 11202215) and the Youth Innovation Promotion Association of CAS (Grant number 2015015).