Numerical Modelling of the 3D Unsteady Flow of an Inlet Particle Separator for Turboshaft Engines

: Helicopter and turboprop engines are susceptible to the ingestion of debris and other foreign objects, especially during take-off, landing, and hover. To avoid deleterious effects, ﬁlters such as Inlet Particle Separators (IPS) can be installed. However, the performance and limitations of these systems have to be investigated before the actual equipment can be installed in the aircraft powerplant. In this paper, we propose different numerical methods with increasing resolution in order to provide an aerodynamic characterization of the IPS, i.e., from a simple semi-empirical model to 3D large eddy simulation. We validate these numerical tools that could aid IPS design using experimental data in terms of global parameters such as separation efﬁciency and pressure losses. For each of those tools, we underline weaknesses and potential beneﬁts in industry practices. Unsteady ﬂow analysis reveals that detached eddy simulation is the trade-off choice that allows designers to most effectively plan experimental campaigns and mitigate risks.


Introduction
The ingestion of dust, sand, and other debris is one of the most threatening issues for helicopter engines when operating in extreme conditions.Particles jeopardize the integrity of the internal components, potentially inducing blade erosion and combustor wall glazing [1,2].Inlet Particle Separators (IPS) are one method to avoid particle contamination in helicopter and turboprop engines, and they have become popular thanks to their compactness, low pressure loss, and weight savings [3] when compared to other filtering systems such as vortex tube separators and inlet barrier filters [1,3].The working principle of an IPS is quite simple; the unfiltered airflow, composed of a gaseous phase and a particulate phase, is directed through a bifurcating channel.The rapid change of curvature in the geometry forces the particle trajectories to change and concentrate only into a small part of the gas flow, while the remaining debris is separated by impingement onto the lateral surfaces and redirected towards a bypass channel.As a result, the air flowing into the main channel is cleaner, and the process can be enhanced either by increasing the geometry curvature or by introducing additional obstacles to the flow path, such as flaps.The described features are visualized in Figure 1 (left), which provides a schematic representation of the geometry inside an IPS built at the von Karman Institute for Fluid Dynamics (VKI).This model is designed for the latest generation of turboprop engines, in which the presence of additional barriers (the two flaps) fosters the separation process [4].
Despite the benefits of these filters, IPS designers need to tackle several issues, such as additional weight and drag, increased engine power requirements, the constant need for inspection, and supplementary costs for logistics and installation.Moreover, the IPS geometry has to be carefully defined to avoid excessive pressure drop in the flow path from the inlet to the outlet channel, which is turbulent, highly 3D, and unsteady.
In response to the industrial need for IPS characterization and validation, researchers have devoted a great deal of attention to assessing separation efficiency through sand injection tests [5,6] and to the aerodynamic characterization of the IPS using CFD tools [7][8][9].Barone et al. [3,5,6] performed experimental characterization of a 2D IPS by injecting sand of various sizes.The procedure consisted of weighing the filtered sand and comparing it with the mass of the total number of particles entering the intake.However, the aforementioned studies are limited to a bidimensional test section in a vertically-designed wind tunnel.Connolly et al. [8] recently published the first 3D simulation of IPS behaviour, including both particle tracking and detached eddy simulation (DES) approaches.The simplified geometry from the previous studies of Barone et al. [3,5,6] was extruded to account for three-dimensional flow investigations.The results clarify many flow features that were previously uncovered.Nonetheless, the study focused on a simplified bifurcating channel with neither flaps nor complex 3D surfaces.
To the best of the authors' knowledge, most studies available in the literature devote scarce attention to the filters actually installed in helicopter and turboprop engines, which are more complex in terms of geometry and unsteadiness.Moreover, there is no systematic exploration of the analytical and CFD approaches that can be employed in the preliminary IPS design phase to support and complete experimental campaigns.
Therefore, this paper aims to provide a set of numerical models for the aerodynamic characterization of an industrial IPS (Figure 1, right); these models have increased level of representativity for the flow and computational effort.The semi-empirical and steadystate models allow designers to carry out a preliminary assessment of the flow features at the boundaries of the domain (bypass, outlet, and inlet channel), while the unsteady simulations deepen the flow separation and other important phenomena.The different approaches are validated with experiments in terms of global parameters (pressure losses and separation efficiency) to ensure their potential for use in parallel with experimental campaigns.The final goal is to define the trade-off method for the IPS design phase, i.e., a high-fidelity approach which helps to reduce the time of experimental campaigns while at the same time being characterized by reduced computational effort to meet the most stringent industrial timescales.The remainder of this paper is structured as follows.Section 2 describes the experimental apparatus used for determining the separation efficiency and pressure losses of the IPS system.The numerical setup is presented in Section 3, first discussing the CFD models and then the semi-empirical 1D model.The results presented in Section 4 provide insights into the validation of the numerical models and the unsteady flow features inside the IPS.Finally, Section 5 concludes the paper with a discussion on the numerical modeling of the system, underlining weaknesses and potential benefits of each model in industrial practice.

Experimental Layout
The IPS designed for testing (shown in Figure 2, left) is installed in the C3 facility at the VKI, which is depicted in Figure 3; this filter is the scaled test section of an IPS manufactured by Safran Helicopter Engines in the framework of the EU-supported IPANEMA programme (Inlet PArticle Separator Numerical ExperiMental Assessment).The rig is a blow-down facility fed by compressed air at 40 bar in a 72 m 3 reservoir.This allows typical test runs of 5-10 min depending on the mass flow required in the rig.The inlet vessel has the main function of stabilizing the flow coming from the pressure line system.The filtered air and the air flowing in the bypass are discharged into two downstream vessels, where sand collection is easier and more efficient.Each of the vessels is equipped with a frame mounting a particle filter.The reason for using vessels at the exit of both streams is to reduce the velocity, lowering the inlet pressure required to achieve the target mass flow and allowing the particles to settle on the filters as the flow speed is reduced.The locations of the measurement sections are indicated in Figure 2    The IPS is tested over a range of flow conditions as close as possible to the operating conditions of a helicopter or turboprop engine in harsh environments and/or off-design situations.In this facility, these conditions are realized mainly through variation of the inlet total pressure p 0 I N and outlet Reynolds number Re OUT , defined as follows: where ṁOUT is the mass flow rate, µ is the dynamic viscosity, and P OUT is the wetted perimeter, all of which are assessed in the cross-section at the outlet of the IPS (Figure 2, right).In order to compare different test cases, the reduced mass flow at the outlet can be defined as where T 0 OUT and p 0 OUT are the total temperature and pressure at the outlet, respectively.The mass flow rates in the channels are measured with Venturi nozzles, while Re OUT is manually regulated using a butterfly valve located after the outlet mass flow meter.There is another control valve installed downstream from the bypass mass flow meter, which is operated automatically and sets the bypass ratio (BPR) at each test (visible in Figure 3).The BPR is provided by the ratio of the mass flow rate measured in the bypass channel ( ṁBY ) and the total mass flow rate injected at the inlet ( ṁIN ): Moreover, the total quantities are measured in the upstream and downstream vessels during each test (locations depicted in Figure 2, right), allowing the total pressure losses ζ OUT/BY to be calculated: where the subscripts for the total pressure p 0 indicate either the outlet (OUT) or bypass (BY) exit cross-section, while ρ I N and v I N respectively indicate the density and velocity measured at the inlet cross-section.To validate the numerical models, the relative total pressure difference PL OUT/BY is taken into account as follows: Regarding the particles, within the scope of this paper only one size range is employed; more details will be presented in a future work.We selected soda lime glass microspheres of 20-30 µm and density 2500 kg/m 3 ; due to their reduced size, they represent a critical test case.This fine sand is characterized by a Stokes number in the range of 2-5, computed as the ratio between the particle time and flow characteristic time.The experimental test matrix is summarized in Table 1.The separation efficiency η is calculated as follows: Weight of sand recovered from the bypass channel Total weight of sand injected from the inlet channel .The uncertainty of the derived quantities (e.g., mass flow rate, pressure loss, separation efficiency, Mach and Reynolds number, etc.) are calculated using Taylor's formula [10], while the absolute uncertainties of the measurement devices are defined through calibration.The uncertainty of the separation efficiency is strongly dependent on the uncollected mass of sand after each test, which is mainly caused by leakages and manufacturing imperfections.

Numerical Setup
For the success of an experimental campaign, it is necessary to forecast potential issues during the IPS testing.Analytical models can provide rapid assessment of the IPS flow conditions and allow for leaner design procedures, as the designer does not need to resort to CFD tools to verify the effect of small variations in the initial conditions or in the IPS geometry.In this paper, we propose a semi-empirical 1D model (Figure 1, right), which is combined with CFD data to enhance its accuracy.However, as the method does not resolve the physics inside the IPS, it should be used in conjunction with higher-resolution methods.For investigations on industrially relevant problems, Reynolds-averaged Navier-Stokes (RANS) simulations are most frequently used.However, their main limitation is that all quantities are averaged; therefore, the obtained information is limited to averaged fields, i.e., data on secondary statistics cannot be straightforwardly obtained via RANS simulations.On the other hand, unsteady RANS (URANS) can offer time-accurate solutions; however, these remain limited to the unsteady fluctuations in the mean quantities.Consequently, RANS statistical approaches might not be suitable for the description of the IPS system, which is characterized by 3D, turbulent, and unsteady separation regions.Time-resolved simulations such as DES and large eddy simulation (LES) are capable of resolving those areas; however, they are significantly more expensive in terms of computational cost, especially for high-resolution grids.Due to the high Reynolds numbers (Table 1) and large size of the IPS domain, direct numerical simulations are computationally prohibitive, and have not been performed.The commercial software programs Fluent 2021R2 and CFD-POST 2021R2 are used for simulations and postprocessing respectively, while the 1D model is implemented in Python 3.9.7.

CFD Models
As a first step, we generated several tetrahedral meshes, with the scope gradually increasing the resolution in the streamwise, transverse, and spanwise (x, y, and z) directions.The grids are refined in the areas with high curvature without refining the boundary layer region with prismatic layers, as this paper focuses on the core flow regions.The fact that the regions close to the walls are not resolved certainly implies some uncertainty in the identification of the exact separation points, and small reattachment bubbles and vortexes close to the walls are not well captured.Nevertheless, this modeling proves to be effective for representing the flow features inside the IPS (Figures 4 and 5).Details on the generated meshes are available in Table 2; the values of y + refer to k-ε RANS calculations at Re OUT = 580,000, p 0 I N = 124,200 Pa, BPR = 21%, and G OUT = 14.4 kg s K 0.5 bar .In the two extreme cases, i.e., the coarsest and finest grids of around 8 million (M) and 95 M cells, respectively, the mesh is overall refined by around two times in each direction.However, excessively refining the grid up to 95 M cells causes a reduction of the orthogonal quality, mainly related to the absence of inflation layers in the high curvature walls, which implies an increase of the aspect ratio in those areas.Nonetheless, the enhancement of the resolution of the boundary layer for the 95 M mesh would make the CFD solution too expensive from the computational viewpoint.By focusing on the y + distribution in the fluid domain, y + is shifted towards lower values for more refined meshes.Regarding the coarsest mesh (8 M), although the maximum y + value is quite high, only 4% of the cells have y + > 300, i.e., most of the domain is in the suggested limit for the k-ε model.Table 2 provides information on the relative total pressure losses (PL OUT ) and approximate computational time for the different grids.The latter is expressed as a function of the time needed for the RANS calculation with the 8 million cell mesh (φ), regarded as a common industrial practice; φ is around 40 min with 64 CPU cores.The difference between PL OUT predicted by the coarsest mesh and the finest mesh is less than 5%.Preliminary RANS simulations allow the initial LES resolution factor S [11] to be estimated as the ratio between the integral length scale and the grid size: where l 0 is the integral length scale, ∆ = (cell volume) 1/3 , k is the turbulent kinetic energy (TKE), and ε is the dissipation rate.Practically speaking, the resolution factor provides an indication of how many cells are surrounding an eddy of size l 0 .This can provide only an initial estimate, as the length from the RANS computation is an approximation of the integral length scale specifically adapted for homogeneous isotropic turbulence.As a rule of thumb [12], S should be greater than 5 to resolve at least 80% of the TKE.By plotting the contours of S in the symmetry plane (Figure 4), it can be observed that the mesh resolution increases with the grid size, i.e., the area in yellow (S > 5) in the separation region is wider.
The different grids can be analyzed by plotting the velocity streamlines and the normalized vorticity contours (Figure 5).In this study, the normalized vorticity ω is defined as where ω is the vorticity magnitude, H is the inlet channel height, and v I N is the inlet velocity; ω varies with Re I N , and for these preliminary RANS simulations the ratio v I N /H is approximately 680.As depicted in Figure 5, all grids predict a separation bubble near the entrance of the bypass channel located in the collector and above the downstream flap.However, the vortexes have different sizes, and the unsteadiness of these phenomena suggests the impossibility of fully predicting the flow behaviour in the IPS with steady RANS simulations.Moreover, ω has a very similar trend and value for all grids.Although Figure 5 and Table 2 demonstrate the grid convergence of the RANS computations, no accuracy conclusions can be drawn, as this analysis is limited to steady-state flow fields.Therefore, two meshes have been chosen to conduct the unsteady analysis and evaluate their differences.On the one hand, we selected the lowest resolution mesh (8 M), as it predicts the presence of vortical structures well and has an acceptable resolution factor S (Figures 4 and 5).On the other hand, the fine mesh of 70 M cells is used for high-resolution LES and URANS, as the improvement in the resolution factor S from 70 M to 95 M is small and mainly related to low-turbulence regions.According to [13], the 70 M mesh is suitable for wall-modeled LES, as its number of grid points is proportional to Re L (Reynolds calculated using the longitudinal length of the IPS L).It is important to underline that neither this grid nor the finest mesh (95 M) is able to accurately resolve the boundary layer: a wallresolved LES would require a number of grid points proportional to Re 13/7

L
. The numerical and subgrid-scale parametrization requirements for a wall-resolved simulation are beyond the relevant computational capabilities, principally due to the high Reynolds number and complex geometry.The technical details of the CFD cases are summarized in Table 3, which provides information on the boundary conditions, mesh, turbulence models, and approximate computational effort.The latter increases nonlinearly with the Reynolds number, BPR, and p 0 I N ; for this reason, a range of φ values is indicated for cases 5 and 6 (unsteady particle tracking).The realizable k-ε model is chosen for the RANS and DES (hybrid RANS-LES) simulations, and wall modelling is used in order to remain consistent with the y + values of the meshes used.To characterize the particle-laden flow, Lagrangian particle tracking is enabled and a Rosin-Rammler particle distribution is defined to be as close as possible to the experimental campaign (d min = 20 µm and d max = 30 µm).The Wall-Adapting Local Eddy Viscosity model is chosen for the LES, as it is suited to the complex wall-bounded flow in the IPS [14].This approach represents the eddy viscosity with a local formulation based on the square of the velocity gradient tensor, and accounts for the effects of both the strain and the rotation rate of the smallest resolved turbulent fluctuations.As a result, no explicit filtering is needed and only local information is required to build the eddy viscosity.The convergence criteria for the residuals are set to 10 −3 .For the unsteady simulations, we monitored quantities such as static pressure and mass flow rate at the boundaries and volume average of the turbulent kinetic energy.The temporal resolution of the unsteady simulations (∆t = 1 × 10 −6 s) was chosen to maintain a CFL condition of order unity, even though implicit time schemes are used, and to capture in time most of the fluctuations of the turbulent eddies that are resolved in space by the different grids.Explicit time integration is not available for the pressure-based solver in Fluent.The monitors allow the minimum simulation time needed to reach convergence to be defined, i.e., 5τ; τ is the through-flow time, defined as the ideal time a fluid particle takes to go past all the IPS: where v avg is the average velocity magnitude in the longitudinal direction and L is the longitudinal length of the IPS.The accumulation statistics for the unsteady characterization started after 5τ and proceeded for at least 5τ.For clarity, the CPU time indicated for the unsteady simulations in Table 3 refers to a simulation time of 5τ.Regarding the unsteady particle tracking, around 15 million particles are injected during a time window of (1/20)τ.
The particles are dispersed in dilute conditions with small volume and mass fractions (<0.1%), and consequently one-way coupling between fluid and sand can be assumed.The first six cases in Table 3 correspond to the experimental cases in Table 1.Regarding case 7, all the meshes in Table 2 were tested.Case 8 was set with the highest inlet total pressure structurally allowed in the VKI facility.Case 9 was considered because it represents a limit condition in terms of Reynolds number; this case presents more critical flow behaviour with respect to the lower Reynolds number cases (as it will be discussed in Section 4.2),and for this reason it was chosen to validate both the low-resolution (8 M) and high-resolution (70 M) grids.

Semi-Empirical 1D Model
In the IPS focus of this paper, the flow is accelerated and forced through a geometrical restriction from the inlet to the outlet, while it is discharged past the bypass channel with an initial contraction followed by a diffuser.Therefore, we can simplify this IPS into a 2D nozzle system, shown in Figure 6, and solve the equations for compressible flows (1D Euler equations) to find a first approximation of the flow behaviour in the IPS.The flow has been modeled as steady-state and inviscid, and the total temperature is assumed to be constant during the filtering process.Moreover, the semi-empirical 1D model is non-isentropic, as the algorithm computes the losses inside the IPS channel with a modelling derived from the specific test section geometry.As inputs, the model receives the BPR, outlet mass flow rate ṁOUT , inlet total pressure p 0 I N , total temperature T 0 (assumed to be constant during the process), gas properties (specific heat ratio γ, specific heat at constant pressure c p , air gas constant R), and geometrical features, i.e., the cross-section of the outlet, bypass, and inlet channels (A OUT , A I N , A BY ), the throat cross-section A BY,T , and the ratio between the orifice diameter and pipe diameter β (used in the calculation of the pressure losses).The user can set the desired convergence criteria for the iterations.The model returns the following variables: (1) the static-to-total pressure ratios p OUT /p 0 OUT and p BY /p 0 BY ; (2) the flow conditions at the boundary cross-sections (density, Mach number, and Reynolds number); and (3) the pressure losses at the bypass ζ BY and outlet and ζ OUT , which are modeled as explained herein.
It can be observed that the total pressure loss can be written in general as the sum of two contributions, namely, the localized pressure losses and the distributed pressure loss.Distributed losses are verified to have a limited contribution to the overall losses (less than 1%), and can be neglected.Local pressure drops mainly occur in three locations, as shown in Figure 6: to a first approximation, the boundary layer effects and the separation phenomena can be neglected.The restriction in Area A (in green) causes most of the losses; however, the collector bend (Area C in light blue) and the small contraction before the diffuser in the bypass channel (Area B in purple) cause minor local pressure drops as well.The localized pressure drop in the bypass exit channel (Area B in purple) is always lower than the one affecting the flow past the outlet channel, and can be calculated as follows: i.e., the pressure losses are 30% of the dynamic head, computed with the arithmetic average between inlet and outlet density and velocity.This approximation has proven to be valid only at low inlet total pressure, when the bypass pressure losses are not as high.This is the case of the VKI experimental campaign.If the pressure increases, the compressibility in that region is no longer negligible, and the losses due to boundary layer and separation become very high.The factor 0.3 in Equation ( 10) is the localized loss coefficient for elbows with a laminar Reynolds number [15], which works fine for a first approximation and is validated by CFD at low pressure.
As confirmed by CFD studies with respect to the flow passage towards the compressor (of primary interest from an industrial point of view), Area A accounts for 80% of the losses and Area C for almost 20%.The former can be simplified as a passage with an orifice.The prescription provided by the ISO [15] can be used to calculate the localized losses in this case: In Equation (11), is the compressibility factor, which depends on the ratio of the pressures between the outlet and inlet channels.
On the other hand, the losses in Area C (Figure 6) can be visualized as two 90-degree bends without considering any flow swirl or vorticity induced by the collector geometrical shape.Therefore, we can write The above formula proves to be valid at low pressure, and is almost independent from the Reynolds number.As the pressure increases, the algorithm fails to provide accurate predictions.
Finally, the pressure losses are provided by

Validation of the Numerical Models
The aforementioned estimation of the resolution factor (Equation ( 7)) can only be used as an a priori calculation, and must be validated a posteriori by focusing on the ratio between the resolved TKE and total TKE from the unsteady simulations (DES and LES). Figure 7B,C depicts this ratio for the DES and LES simulations (case 9 in Table 3).Regarding the LES, in almost all the domains more than 95% of the turbulent kinetic energy is resolved, while the DES (performed with the 8 million mesh) has larger areas where the turbulent kinetic energy is less resolved, for instance, those above the two flaps and before the entrance of the bypass channel.Those regions are mostly assigned the RANS mode of solution, as indicated in yellow by the contour in Figure 7A.This contour is based on the calculation of a user-defined function (UDF): where the two length scales of the realizable k-ε and DES model are respectively provided by l RKE = k 3/2 /ε and l DES = C DES max(∆x, ∆y, ∆z) = 0.61max(∆x, ∆y, ∆z) (maximum local grid spacing).
The models are compared with experiments in terms of pressure losses and separation efficiency in Table 4. Case 7 is not presented, as it has been already discussed (Table 2), while cases 8 and 9 are taken into account for the unsteady flow characterization.
Only the relative total pressure losses in the outlet channel (PL OUT ) are discussed, due to their greater influence on the overall engine performance with respect to those in the bypass channel.In fact, the pressure drop occurring in the mainstream channel directly affects the compressor operating point and its performance.It can be seen that PL OUT increases systematically as the bypass ratio increases, and the influence of the BPR is more relevant at high Reynolds numbers.Although the experimental data may suggest that PL OUT increases with Re OUT for a given BPR, this trend cannot be confirmed by relying only on experiments; due to the facility's architecture, the test conditions are realized by varying both Re OUT and p 0 I N at the same time, and as such the effect of these two parameters is not decoupled.Nonetheless, the semi-empirical model confirms that for given BPR and p 0 I N the pressure losses increase with growing Re OUT ; this is related to the increase in ∆p 0 ORIF and ∆p 0 COLL (Equations ( 11) and ( 12)) as well.Moreover, for a given BPR and Re OUT , the model predicts a decrease in PL OUT as p 0 I N increases.However, the theoretical model shows a larger difference with respect to experiments at higher Reynolds numbers (Re OUT = 710,000) and BPR (>13%), as additional flow phenomena develop that cannot be captured by this 1D model, such as the separation bubble at the entrance of the bypass channel (as it will be discussed in Section 4.2).The pressure losses predicted by the steady-state RANS simulations are different from the experimental results at high Reynolds numbers and BPR, and they fall outside the uncertainty range in cases 4 and 5. On the other hand, the average pressure losses predicted by DES and URANS are closer to the experiments, overcoming the weakness of RANS modeling in characterizing the development of 3D flow instabilities at high Reynolds numbers and BPR.Moreover, the 1D semi empirical model and the steady-state CFD can be used only as preliminary design tools for Re OUT = 1,280,000 (case 9), as they prove inadequate in the estimation of PL OUT (underestimation of 15% compared to LES simulation).The average pressure losses from the URANS model are similar to those from high-resolution methods (DES for cases 4, 5, and 6 and LES for case 9); concerning case 9, it can be observed that the URANS calculations performed on two different grids (8 M and 70 M) show similar results (relative difference in PL OUT of less than 1%).For simplicity, the pressure losses in the bypass channel (PL BY ) are not shown; they are one order of magnitude less than PL OUT , and show an increasing trend with increasing BPR.In the limit case BPR→0, both particles and airflow are directed towards the mainstream channel; hence, the pressure drop in the bypass channel becomes smaller and smaller (PL BY → 0).3).
Regarding the separation efficiency η (Table 4), the outcomes of the experimental campaign show that the separation efficiency is harder to evaluate at low BPR and high Reynolds numbers (higher uncertainty), as the amount of uncollected sand during tests increases; this negatively affects the validation process, especially at high Reynolds numbers.In these flow conditions, the RANS simulations with Lagrangian tracking loses accuracy and reliability; the solver does not fully converge, as particles remain trapped in the fluid domain.Therefore, RANS can only be used as an effective tool for the estimation of separation efficiency at high BPR.On the other hand, URANS coupled with unsteady particle tracking allows the separation efficiency to be predicted in better agreement with the experimental data, although it requires a longer computational time in order to let the injected particles escape from the bypass and outlet channel.It was observed that the particular shape of the collector further hinders the particles and makes the separation process more effective, as will be analyzed in a future publication.

Unsteady Flow Characterization
Although the URANS model accords with the experimental data, the model validation cannot be limited to the global flow parameters.For this reason, the purpose of this section is to analyze several of the flow features inside the IPS at high Reynolds numbers (Re OUT = 1,280,000, case 9) in order to identify the weaknesses and strengths of each CFD model.The LES performed with the 70 M grid is taken as reference case; for brevity, we focus only on the dimensionless total TKE k.The latter is computed as follows: for RANS and URANS where u x , u y , and u z are the respective velocity fluctuations in the x, y, and z directions, k avg is the average TKE, and k sgs,avg is the average subgrid-scale TKE from the WALE model.The definition of the total turbulent kinetic energy for LES and DES in Equation ( 16) assumes that the unsteadiness of the mean flow is negligible when compared to the turbulent fluctuations, as validated through URANS computations.Figure 8 depicts the total TKE k in the separation region (symmetry plane view).While the RANS models (both 8 M and 70 M) correctly identify the presence of turbulent flow structures above the two flaps, they underestimate the TKE of the separated flow, although the "improved" version of the k-ε model (realizable k-ε model) has been used here.Moreover, the location of the high-TKE region above the downstream flap changes with the mesh resolution, suggesting unsteady behaviour and indicating the inadequacy of steady-state methods.Furthermore, the URANS models show a number of differences compared to time-resolved methods such as DES and LES; these are mainly related to the very low magnitude of the captured velocity fluctuations (the resolved TKE is very small).On the other hand, the total k resolved by DES has the same magnitude as the high-resolution LES in the separation regions; the only areas where the model shows discrepancies with LES are those assigned to the RANS mode of solution (Figure 7A), i.e., close to the walls, in which unreliable gradients are predicted.16)) for different models and mesh resolutions (Case 9, Table 3).
The region of interest is the same as the one illustrated in Figure 9.
These considerations indicate that the flow in the IPS is highly turbulent and unsteady; therefore, the unsteady characterization cannot be limited to the mean flow quantities.Because of its considerably lower computational effort and relatively small difference with LES, DES has been taken into account for studying the instantaneous flow fields from the statistics accumulation; Figure 9 illustrates the development of the flow separation bubble at the entrance of the bypass channel when increasing Re OUT .This instability is detrimental for the separation process, and represents an additional threat to the safe operation of the IPS.Designers should be aware of the unsteady development of the bubble as well as of the vortexes developing above the downstream flap (suggested by the streamlines and by the high vorticity regions in Figure 9).Focusing on planes perpendicular to the streamwise direction above the downstream flap, we observed that for both low and high Reynolds numbers (cases 8 and 9), the RANS and URANS simulations underestimated the local pressure and density drop, resulting in the evolution of these profiles (pressure recovery) being more optimistic than the one predicted by the mean of LES and DES.
As illustrated in the plot in Figure 10, the mean quantities are not enough to deeply describe the separation process, as the BPR fluctuates over the acquisition window (15τ-20τ) with a frequency close to the through-flow frequency 1/τ and the targeted BPR (21%) cannot be kept constant throughout the whole separation process.The steady state simulations converge towards a mean BPR of around 21%, while the unsteady simulations have slightly lower average values due to the captured mass flow rate fluctuations (which are minimal in case of URANS).Figure 11 indicates the evolution of several separation bubbles during (approximately) one oscillation cycle.All the 3D vortexes are projected onto the symmetry plane, and their directions are schematically depicted.During a cycle, there small corner vortexes develop and dissipate around the L-strut and downstream flap as well as in the curvature regions of the collector entrance (Figure 1).The vortex above the upstream flap suggests the generation of a horseshoe vortex in which the out-of-plane rotation is almost negligible.On the other hand, the two largest recirculation bubbles above the downstream flap and at the entrance of the bypass channel mostly deploy in the direction perpendicular to the symmetry plane.In particular, the bubble above the downstream flap is visible only at t = 16.25τ;therefore, this vortex does not having the same intensity along the out-of-plane direction.Last but not least, the vortical distribution at the entrance of the bypass channel shows that several bubbles develop perpendicularly to the depicted symmetry plane, and are characterized by different intensities.All of these features can influence the performance of the IPS in terms of separation efficiency and pressure losses.

Conclusions
Before actual installation in helicopter and turboprop engines, the IPS needs to be experimentally and numerically validated.This paper has shown how experimental campaigns can be supported and completed by employing different numerical models of increasing complexity.First, 1D semi-empirical models prove to be effective only for preliminary design calculations and planning of experiments, especially when the geometry of the IPS is not known in advance and the designer needs to perform further iterations to optimize the flow boundaries (inlet, outlet, and bypass).The 1D models fail to predict global pressure losses at high Reynolds numbers (Re OUT = 1, 280, 000) and BPR ( 20%), and for this reason have to be supported by CFD.RANS simulations are beneficial in terms of computational time, as they allow losses and separation efficiency to be estimated (though not in extreme conditions) with low effort.Thus, they have been widely employed to comply with stringent industrial timescales.Nonetheless, RANS offers only a limited overview of the separation phenomena occurring in the IPS, underestimating the pressure losses and neglecting the BPR fluctuations.Moreover, RANS models such as k-ε show large discrepancies in their assessment of the turbulent kinetic energy in comparison with time-resolved simulations.A viable option for the estimation of separation efficiency and pressure losses is provided by URANS simulations, which require approximately 25% time less than DES for a given mesh, and more than three times less if unsteady particle tracking is not enabled.Although URANS has been successfully validated with experimental data, it is inaccurate for describing the unsteady behaviour of the IPS (underestimation of k), as with the chosen turbulence model, time stepping, and spatial discretization the velocity fluctuations in the separation regions of the IPS are not well captured.In these core flow regions, DES and LES provide similar average results regardless of the mesh resolution; for this reason, DES can be implemented for a computationally affordable flow analysis.The possibility of visualizing and describing features such as separation bubbles and turbulent structures is highly valuable for the designers of industrial IPS systems as well as in view of test campaigns, as it is possible to move beyond the averaging process and obtain a large amount of information within a reasonable time frame.Therefore, DES is confirmed to be an efficient tool for characterizing IPS aerodynamics, as it allows the limitations of steady-state simulations to be overcome without the heavy computational effort and memory requirements of a high-resolution LES.

Figure 1 .
Figure 1.(Left) Schematic test section of the IPS built at VKI: 1, upstream flap; 2, downstream flap with an L-strut supporting the flap in the deployment phase; 3, collector.(Right) Numerical models presented in this paper.
(right); the inlet (IN), outlet/mainstream (OUT) and bypass (BY) planes correspond to the same locations in which the boundary conditions for the CFD simulations are enforced.

Figure 2 .
Figure 2. (Left) IPS geometry manufactured with 3D printing technologies.(Right) Detail of the test geometry, instrumentation location, and measurement sections.

Figure 3 .
Figure 3. Global layout of the C3 facility (dimensions in mm).

Figure 6 .
Figure 6.(Left) IPS schematized as a convergent nozzle in parallel with a convergent-divergent nozzle.(Right) The three main loss regions in the industrial IPS.

Figure 7 .
Figure 7. (A): User-defined function UDF defining the regions in which the DES model switches to RANS (yellow).(B,C): Contours of the ratio between resolved turbulent kinetic energy and total turbulent kinetic energy (Case 9, Table3).

Figure 8 .
Figure 8. Contours of k (Equation (16)) for different models and mesh resolutions (Case 9, Table3).The region of interest is the same as the one illustrated in Figure9.

Figure 9 .
Figure 9. Instantaneous ω contours (t = 16.75τ) and velocity streamlines in the IPS regions between the two flaps from the mean of DES calculations (cases 8 and 9 Table 3) for two Re OUT .(Left) Re OUT = 710,000 (v I N /H ≈ 680).(Right) Re OUT = 1,280,000 (v I N /H ≈ 820).

Figure 10 .
Figure 10.BPR in the IPS region depicted in Figure 9 at Re OUT = 1,280,000.

Table 2 .
Summary of the grids generated for the IPS numerical modeling, with statistics of mesh quality: x = streamwise; y = transverse; z = spanwise; AR = aspect ratio; OQ = orthogonal quality.

Table 3 .
Numerical test matrix.PT = particle tracking.p 0 I N in [Pa]; G OUT in