Large Eddy Simulation of the By-pass Transition Process under Different Inlet Turbulence Conditions

The transition process of the boundary layer developing over a flat plate with elevated inlet Free Stream Turbulence Intensity (FSTI) has been studied by means of Large Eddy Simulation (LES). To this purpose, four cases with different inflow disturbances have been tested varying the magnitude and the length scale of turbulence. LES has been performed by using the finite-volume ANSYS Fluent code. The computational domain, which was constituted by a rectangular domain with a zero thickness plate, was based on an ERCOFTAC test case in order to provide a validation with a well-known set of data by comparing the boundary layer integral parameters and mean and fluctuating streamwise velocity profiles. The four cases were discussed within the paper by looking at classical statistical properties as well as advanced post-processing tools. It was shown that the decrease in the free stream turbulence level postpones the transition location, whereas the variation of the integral length scale has a very low influence on the distribution of the time-mean flow properties. Proper Orthogonal Decomposition (POD) has been applied to the instantaneous LES flow fields in order to provide a statistical representation of the structures responsible for transition and their response to free-stream turbulence intensity and length scale. The presence of vortical filaments parallel to the wall, typically referred as boundary layer streaks, is clearly identified; their characteristic dimensions and how they change as a function of FSTI properties were analyzed within the paper.


Introduction
The capability of CFD solvers in predicting real flows encountered in engineering applications strongly depends on the assumptions made for the closure of the Reynolds shear stress, requested into RANS based schemes. In this context, transitional boundary layers are well known to be difficult to be predicted, due to the variety of coherent structures that are generated during transition (e.g., [1][2][3]).
Thanks to the advancement of hi-fidelity simulations and advanced optical measurement techniques, it is nowadays well known that the by-pass transition process is driven by the breakup of streaky structures [1], while the separated flow transition mechanism is mainly driven by the formation of Kelvin-Helmholtz (KH) rolls dominating transition (e.g., [4]). From a mathematical point of view, the former represents a non-modal optimal perturbation to the Navier Stokes operator (see Jacob and Durbin [5], Schmid and Henningson [1] for example), while the latter is well represented by the dominating eigenvalue of the same operator in the limiting case of infinite Reynolds number, also known in the literature as a Rayleigh problem. In the case of by-pass transition, the elevated free-stream oscillations may (or may not) penetrate into the boundary layer through a "sheltering" process, as described first in Jacob and Durbin [2] and then also in Zaki and Durbin [6] for variable pressure gradient conditions. Basically, high frequency oscillations cannot penetrate into the boundary layer that consequently acts like a filter admitting only low frequency oscillations. The low frequency oscillations penetrating into the boundary layer are then amplified, assuming the shape of elongated low and high speed filaments, typically referred as streaks (e.g., [7]). The optimal perturbation theory of Luchini [8] provides self-similar distribution of the velocity root-mean square along the wall-normal direction, and it has been proved to work also in the case of adverse pressure gradient [9,10]. The characteristic dimensions of streaky structures like the spanwise distance and streamwise wavelength represent important information that could be used to set boundary conditions to RANS based CFD solver [11]. The works of Mans et al. [12] and that of Schlatter et al. [13] provide a statistical representation of the main geometrical parameters of a streaky pattern preceding breakup. Similarly, in the work of Lengani et al. [10,14], the streak spacing in a realistic low pressure turbine application with and without incoming wakes is provided by means of Proper Orthogonal Decomposition (POD) and auto-correlation analysis. In the present work, LES has been carried out to explore the effects due to the systematic variation of the free-stream turbulence intensity magnitude and length scale on the transition process of a flat plate boundary layer. The reference condition has been set to match the ERCOFTAC data reported in Roach and Brierley [15], thus providing validation of the accuracy of the present work, as provided by Voke and Yang in several works [16][17][18] in order to validate their LES results.
Successively, both the turbulence intensity magnitude and the length scale have been varied independently in order to clearly highlight the effects of these parameters on the statistical properties of the transition process, and on the coherent structures involved in the process itself. Particularly, the capability of POD in recognizing the dominant source of oscillation into a set of snapshot has been here invoked to provide a statistical representation of the most energetic oscillations events affecting transition (see [10,19] for example). The auto-correlation of instantaneous snapshots has been computed to identify the effects of the parameter variation on the streak spacing, following the work of Matsubara and Alfredsson [20].

Simulation Approach and Data Reduction
Large eddy simulations have been performed to analyze the influence of turbulent free-stream on the development of zero pressure gradient transitional boundary layers. Both the computational domain and the turbulent flow condition at the leading edge of the flat plate have been chosen with reference to the experimental test case ERCOFTAC T3A [15], as also adopted in Rai and Moin [21] and Jacobs and Durbin [5].
The computational domain is a rectangular box whose dimensions in the streamwise, wall-normal and spanwise directions are L x =1666 mm, L y =87 mm and L z =65 mm, respectively. The domain inlet is placed at 250 mm upstream of the flat plate leading edge, similarly to Rai and Moin simulations [21], while the outlet position and the other two dimensions of the domain have been chosen according to Jacobs and Durbin [5]. In particular, the lateral size of the domain is sufficiently large to avoid influencing the average streak spacing.
The computational grid used in the present work consists of 6.02 million nodes uniformly spaced in the streamwise and spanwise directions and clustered near the lower boundary in the y-direction. More specifically, the grid points in each direction are 806, 89 and 84 giving a resolution in wall units of x + =36, y + varying from 0.95 at the wall to 37 at the upper boundary and z + =14.
The resulting Courant number in all simulations was maintained lower than 0.25, as also indicated in Ref. [18]. All the simulations described here have been performed adopting the same kind of boundary conditions and the same numerical set up. Upstream of the flat plate, where free-stream turbulence develops, both the lower and upper boundaries are treated as symmetric. Then, a no-slip condition is imposed on the plate and periodicity is assumed in the spanwise direction. Ambient pressure is specified at the remainder part of the upper boundary and at the outlet. On this last surface, it is prescribed as average value. At the inlet, a uniform velocity profile is specified (u=5.4 m/s, v=w=0 m/s .

Fig. 1 Computational domain
On this mean profile, fluctuating velocity components, which are computed by means of a spectral synthesizer based on a random flow generation technique [22], are superimposed. These perturbation velocities are linked to the spectral synthesizer turbulent intensity (SSTI) and length scale (λ 0 prescribed at the inlet. Table 1 summarizes these values along with the calculated FSTI and length scale at the leading edge of the flat plate for the four cases analyzed in the present work. All the large eddy simulations have been performed by means of the finite-volume code Fluent from ANSYS [23]. More specifically, the SIMPLEC pressure-based solver has been used. As regards spatial discretization, the bounded central differencing scheme [23] has been adopted in order to avoid the appearance of unphysical oscillations in the flow field, even if the unbounded central one [23] would be an ideal choice thanks to its low numerical diffusion. Moreover, it has been adopted a second order implicit transient formulation for the advancement in time and the WALE subgrid scale model. Finally, it is worth noting that the time step for each calculation has been set to 0.1 ms. The flow field has been allowed to develop for seventeen throughflow time before starting data collection for approximately 1 s of physical flow time at 1 kHz sampling rate. These latter data have been further analyzed by means of POD.

Proper Orthogonal Decomposition
The Proper Orthogonal Decomposition (POD) is nowadays a well-established mathematical procedure which has been largely used for the detection of coherent structures embedded within the flow (since Lumley [24]). The POD is performed in the present work following the mathematical procedure described in Sirovich [25]. According to this procedure, given N LES snapshots, the POD can be computed looking for the eigenvalues  and the eigenvectors X of a cross-correlation matrix C[NxN]. The i, j element of the matrix C is defined as the volume integral of the inner product of the vector fields at times i and j. This operation can be performed in matrix form collecting the LES data in a proper matrix V N . The rows of V N are built in order to contain the velocity information in the spatial domain for each snapshot; the columns of V N contain the temporal information, i.e., the N=1000 snapshots in the present case. The matrix C is then computed as where the matrix M is diagonal and its entries are given by the local cell volume (see [26]). This operation is particularly convenient for LES data since the rows of V N have dimension O(10 7 ), while its columns have dimension O(10 3 ), consequently the matrix C has then square dimension O(10 3 ).
The POD modes are finally obtained projecting the original data on the computed eigenvectors with a scalar product involving the weighting matrix M. More details on this procedure may be found in a large amount of works in literature (e.g. [27][28][29]).The POD provides a triplet of information: the eigenvalues , the eigenvectors X, and the POD modes . The eigenvalue of the i th mode ( (i) ) represents the energy contribution of the mode to the total kinetic energy of velocity fluctuations in the whole measuring plane. The eigenvector of the i th mode ( (i) ) retains the temporal information related to each mode. The POD modes constitute an orthonormal basis that provides the spatial information identifying coherent structures in the flow. Hence, the method extracts and separates spatial from temporal information determining the most energetic structures and ordering them by their energy content. Since the energy rank is based on the total kinetic energy, the amplitudes of the POD modes of the streamwise ( u ) and of the wall-normal ( v ) velocities are not equally distributed: i.e., a POD mode may show large fluctuations of the streamwise velocity ( u ) and null fluctuations of the wall-normal velocity ( v ).

Auto-correlation
In order to characterize the streaky structures within the boundary layer, the auto-correlation of the instantaneous perturbation velocity fields is computed similarly to what has been done in Ref. [10]: where  u is the standard deviation of the velocity component u, and ∆z is the distance between two observations along the spanwise direction. The notation A B  represents an ensemble average over the 1000 LES snapshots (i.e., using the same ensemble of data where POD has been computed). This operation is performed after applying a POD filter that removes the most energetic POD modes that are linked to large scale structures as discussed in the result section.

Results and Discussion
The instantaneous velocity distribution on the wall-normal plane at the center of the computational domain is reported in Fig. 2  The distribution of the shape factor H 12 ( Fig. 5(a)) provides a clear picture of the boundary layer state for the different flow conditions. Because of the high turbulence level, the shape factor is above 2.5 for a very short length and starts to decrease just downstream of x=0.1 m for all conditions. Cases 1 and 2 show the fastest decrease of the shape factor, coherently with the average results of Figs. 3 and 4. Particularly, for these two cases in the position of maximum RMS within the boundary layer region, the shape factor is below 1.8 indicating that the boundary layer state is rapidly approaching the fully turbulent condition. This condition is slightly delayed for Case 3. The further reduction of the inlet FSTI leads to a further delay of the transition process. The shape factor is above 2.4 up to x=0.4 m even though the RMS of the streamwise velocity fluctuations is above 7% within the boundary layer. This case, as it will be shown in the following, will provide the clearest visualization of the streaky structures that populate the boundary layer before transition. At the end of the measurement domain in analysis, Cases 1-3 are in a fully turbulent condition (H 12 =1.45 , while Case 4 is still not. In Fig. 5(b), the shape factor distribution is plotted against the logarithmic mapping of the local momentum thickness Reynolds number. This plot is provided in order to discuss a direct comparison of the present data against the ERCOFTAC test case of Roach and Brierley [15], since they provided the same plot in their paper.   for Case 2 at different streamwise locations for LES data (continuous lines) and data from ERCOFTAC test case (dots) [15] Therefore, since a synthetic turbulence generator has been used, different from the two turbulence generating grids considered by Roach [15], the present data can be considered successfully validated. In Fig. 6, the velocity profiles at Re x =32 400 show perfect agreement and a laminar like distribution. Also the rms u is well captured. Good agreement is also shown at Re x =273 500, where the boundary layer is in a fully turbulent condition. The largest discrepancies are in the transitional region (grey lines), where LES results show a slightly more advanced transition process (the peak in rms u is higher than experimental data).

POD analysis
Proper orthogonal decomposition has been applied to the instantaneous LES flow fields in order to provide a statistical representation of the streaky structures that lead to transition and their response to free-stream turbulence intensity and length scale variation. Fig. 7 shows the POD eigenvalues normalized with the total kinetic energy for the different cases. Case 4 presents higher values of normalized energy for the first 10 modes with respect to other cases, while for all cases the energy of the modes above the 10 th becomes rapidly lower than 1%. Fig. 8 provides a 3D view of the first, third and seventh POD modes of the streamwise velocity component for the four different cases in analysis. Modes are shown from top to bottom of each case in decreasing order; the first plot on top is the first mode, the most energetic one; the second is the third mode and on bottom the seventh mode (where for all cases the energy is the same). Iso-surfaces at low or high values of the POD mode are shown to identify the most energetic structures. The origin of the plate is marked on the left bottom corner of each plot displaying the 3D axes.
The structures of the POD modes under analysis are very similar for the first three cases. The first mode highlights large scale structures that are elongated in the streamwise direction. These structures are clearly larger than the typical streaks dimension that has been presented in Fig. 2. The high-energy POD modes (the first one and the third one as well) may statistically represent "agglomeration" of these elongated structures as also discussed in Lengani and Simoni [19].
Furthermore, the modes of the first three cases show the largest values in the front half part of the computational domain since in the aft part the boundary layer is turbulent and large scale structures broke down. In this part of the plate, finer scale structures, with a significantly smaller elongation in the streamwise direction, may be observed for Cases 1-3.
Case 4 shows again the most notable differences with respect to the other cases, since at this lowest turbulence level the largest scale structures (Modes 1 and 3 of Fig.  8(d)) are localized in the aft part of the computational domain. Mode 3 shows very long ordered structures that are elongated in the streamwise direction in the middle of the computational domain and that maintain their coherence up to the exit of the flat plate.
The higher order modes (the bottom plot of each subfigure) provide a very similar representation of smaller scale that populate the first half of the LES domain for Cases 1-3 and are localized further downstream for Case 4. Due to the delayed transition, Case 4 exhibits streamwise oriented structures also on the 7 th mode.

Instantaneous flow field
Sequences of the LES snapshots that have been filtered by means of POD have been analyzed by looking at their instantaneous distribution and performing auto-correlation along the z-direction. Figs. 9 and 10 show the perturbation streamwise velocity (u′ according to the Reynolds decomposition) for Cases 2 and 4, respectively.
These velocity contours are obtained within the boundary layer region at a fixed y coordinate, where the maximum RMS is observed. Cases 2 and 4 only are considered in the following discussion since the differences between the first three cases are limited. Boundary layer streaks can be clearly observed in both series of images.
They may be identified as small structures that are elongated in the streamwise direction with a perturbation velocity magnitude that is about 20% of the freestream velocity; the blue stripes represent structures that move slower than the mean flow, while the red ones move faster. When the FSTI is high (Fig. 9), the streaks populate the boundary layer in the range 0.18 m<x<0.4 m.
Interestingly, low speed filaments are more present into the flow. Downstream of x=0.4 m most of the streaks undergo break down thus providing the largest turbulence generation. In fact for this case the maximum RMS was observed at x=0.45 m. Downstream of this position the perturbation velocity assumes a random pattern with nuclei at positive and negative velocity that loose a preferred elongation direction, while they tends to assume a circular shape.
The propagation of the streaks within the boundary layer for the low turbulence case (Fig. 10) is different. Boundary layer streaks maintain their coherence in the range 0.2 m<x<0.6 m for a much longer distance than in the high turbulence case. Also in this case, the low speed streaks are more clearly recognizable and more frequent than high speed ones. The region of breakdown of the streaky structures is downstream of x=0.6 m, leading to the maximum RMS at x=0.65 m, coherently to what observed in Fig. 4. Furthermore, the magnitude of the perturbation velocity u′ seems slightly lower than that observed for the high FSTI case.

Auto-correlation analysis
The auto-correlation analysis (e.g., [20]) has been performed on the perturbation velocity field after applying the POD-filter previously discussed. According to Matsubara and Alfredsson [20], the locus of the minimum of the auto-correlation is strictly connected to the distance between a high speed streak and the neighbouring low speed streak. This is used as characteristic dimension of the streaky structures. It has been also shown that this dimension often scales between different conditions with the boundary layer integral parameters. In the present case, the auto-correlation analysis has been performed along the z-direction for different axial position at fixed y coordinate averaging the results among the 1000 LES snapshots for each condition.
The  Particularly, once the auto-correlation is scaled with the local value of the displacement thickness, the streaks spacing (minimum of the red lines) is about 3. After breakdown to turbulence, the auto-correlation in the spanwise direction changes and provides a minimum that is below 2 for each condition. Similar results have been also found for the other two conditions. Therefore, the typical dimension of the boundary layer streaks may be scaled with the displacement thickness and be self-similar among different inflow conditions (i.e. varying the turbulence or its length scale).

Conclusions
LES data of the boundary layer developing over a flat plate have been discussed in the paper. This study concerned four different cases where the inlet free-stream turbulence intensity and length scale have been varied. The numerical results have been validated against a well-know ERCOFTAC test case by comparing the numerical and experimental distributions of the boundary layer shape factor and mean and fluctuating streamwise velocity profiles. The LES data have been further analyzed by means of time mean properties and providing advanced statistical analysis by means of POD and auto-correlations. The increase of the free stream turbulence intensity leads to a faster transition. Particularly, it has been shown that the transition is anticipated for the high turbulence cases by a prompter breakdown of the boundary layer streaks. Such structures, that are initially ordered and elongated along the streamwise direction, may be represented by the POD modes that identify the regions where the largest coherent velocity fluctuations occur. The POD modes identify even larger scale structures when the turbulence is decreased. A similar description of the flow physics is provided by sequences of the instantaneous velocity distribution. The peak RMS of velocity fluctuations corresponds to region downstream of the streak breakdown, while, apparently the length scale of such structures is different among the cases in analysis. However, the scaling of the characteristic dimension of the streaks with the boundary layer displacement thickness shows a self-similarity property among the different inflow condition.