Validation of large eddy simulation of flow behind a circular cylinder

The flow in the wake behind a circular cylinder in a cross-flow at Reynolds number of 4815 was studied both experimentally and via mathematical modeling. The mathematical model was performed as a Large Eddy Simulation (LES), while the experiments were carried out using the time-resolved variant of the Particle Image Velocimetry (PIV) method. Both the simulation and experiment took into account the dynamical aspects of the studied phenomenon, which enabled a detailed validation of the mathematical model. The overall statistical properties of the simulated flow were validated via comparing the time-averaged measured and computed velocity and vorticity fields. To validate the dynamical behavior, the velocity spectra were examined first. Next, the Proper Orthogonal Decomposition (POD) of the spatio-temporal velocity data was performed on both the experimental and numerical data and a comparison of the obtained energetic modes was carried out. All the performed validations have shown a satisfactory agreement between the simulation and the experiment.


Introduction
Two strong trends have been observed in fluid dynamics during the past decade. First, the increasing available computing power have enabled applications of high-fidelity modeling approaches, such as the large eddy simulation (LES), even to industrial-scale problems [1][2][3]. Second, improving measurement techniques, e.g. the particle image velocimetry (PIV), provide more and more detailed data on flow dynamics. Nevertheless, it seems that even the dynamic LE simulations are commonly validated against static data, see e.g. [2,[4][5][6]. Furthermore, if any analysis of the system dynamics is performed, it is usually limited to the comparison of turbulence spectra as in [7,8]. Such a situation follows from the fact that most existing verification and validation methods for fluid flow simulations are derived for the computationally cheaper and better established Reynolds-averaged Navier-Stokes equations (RANS) [9].
In the present paper, we aim to promote a discussion on LES data validation. Especially with respect to the flow dynamics and validation of the simulated turbulent structures. A geometrically simple set-up of a wide and thin circular cylinder in a cross-flow is used in both experiment and LE simulation. The experiment is designed in order to provide detailed data on flow dynamics and is similar as published in [10,11]. Furthermore, the simulation is prepared with a level of detail corresponding to the industrial LES, i.e. the whole problem geometry (experimental test section) is not modeled using the computationally expensive LES. Instead, we employed the k − ω shear stress transport (SST) detached eddy simulation (DES) model based on the work [12], which is a hybrid LES-RANS model.
The selected modeling approach provides high-resolution data in the zone of interest corresponding to the measured region. Both experimental and numerical data are, similarly to [7,8] analyzed with respect to (i) time-averaged quantities, and (ii) spectra of velocity components at given locations. Furthermore, the measured and modeled flow fields are studied via the means of proper orthogonal decomposition (POD) [13][14][15][16], where we focus on the most energetic coherent flow structures. Note that POD analysis of both experimental [10,11,17] and simulation [18][19][20] results is present in the literature. However, the possibility of comparing measured and computed POD results in order to validate flow simulation seems not to be well studied. The studied problem comprises the cross-flow around a circular cylinder. The experimental data were measured at a blow-down facility at the Institute of Thermomechanics of the Czech Academy of Sciences with a closed test section of 250 × 250 mm 2 . The cylinder model of diameter d = 15 mm and height h = 250 mm (cross-section blockage of about 6 %) was placed at the test section center and oriented transversely with respect to the main flow direction. The situation is schematically depicted in Figure 1a.

Problem geometry and experimental setup
The test section inlet nozzle provides a flow field with a low turbulence of intensity I ≈ 0.2 %, a good regularity with departures below 1 %, and mean velocity U in = 5ms −1 resulting in the Reynolds number Re = U in d/ν = 4815, where ν is the fluid (air) kinematic viscosity. The inlet nozzle provides a well specified flow at the contact with the cylinder model. However, it causes a geometrical discrepancy between the experimental and simulation setup. In particular, the nozzle itself cannot be modeled as the conditions at the nozzle inlet are not known. At the same time, the region in front of the cylinder needs to be included in the computational domain. To extend the computational domain in front of the cylinder while not having to resolve the nozzle inlet, the test section was prolonged in the upstream direction. However, the test section walls in front of the cylinder were considered as impermeable but not posing any tangential resistance to the flow. The situation is depicted in Figure 1a, where the special computational boundary (virtual walls) is shown in red.
The experimental velocity fields were obtained using the Particle Image Velocimetry (PIV) method. The measurement apparatus consists of a laser and CMOS camera by the Dantec company. The laser is the New Wave Pegasus, Nd:YLF double head with wavelength of 527 nm, maximal frequency 10 kHz, shot energy of 10 mJ (for 1 kHz) and the corresponding power of 10 W per one head. The camera Phantom V611 has the resolution of 1 280 x 800 pixels, is able to acquire double snaps with frequency up to 3000 Hz (at full resolution), and it uses an internal memory 8 GB. The data were acquired and post-processed in the DynamicStudio software [21]. To visualize the flow, the SAFEX particles, oil droplets 1 μm in diameter, have been used.
The PIV measurements were performed in the plane perpendicular to the cylinder axis and parallel to the main flow direction, i.e. in the x − y plane as indicated in Figure 1. The evaluated velocity field comprised 159 × 99 vectors and was acquired at the sampling frequency of 2 kHz. One measurement consisted of 4000 double-snapshots representing 2 s of the real-time flow. More details on the experiment may be found in [10,11].

Flow governing equations
The flow in the problem geometry shown in Figure 1 is assumed to be governed by Navier-Stokes equations for an incompressible Newtonian fluid of the standard form where u is the fluid velocity,p = p/ρ is the kinematic pressure, where p is the pressure and ρ the fluid density. Finally, T is the viscous stress tensor. In general, the turbulent flow was solved within the large eddy simulation framework in which u andp are replaced by their filtered counterparts u and p and the stress tensor T is decomposed as where σ is the filtered, or resolved scale, strain rate tensor and τ is the unknown subgrid scale (SGS) stress tensor that represents the effects of SGS motion on the resolved LES fields [1]. Furthermore, in order to correctly and efficiently account for the computational domain being relatively large compared to the region of interest containing the plane of measurement, we use a hybrid LES-RANS turbulence model based on k − ω SST approach introduced by Strelets [12]. This particular approach known as the detached eddy simulation (DES) turbulence model automatically adapts between LES and RANS computation in dependence on a DES length scalel, which is determined by the local cell size Δ and turbulent kinetic energy (TKE, k) as described in [12,22]. The flow governing equations were solved utilizing the open-source CFD library OpenFOAM [23], which is based on the finite volume method.

Computational mesh and boundary conditions
The computational mesh was constructed such that the region of interest (RoI: as well as the domain parts directly affecting the flow in RoI comprise a structured hex-dominant grid with a resolution sufficient to fully simulate all the relevant flow scales. On the other hand, the main purpose of the computational domain outside of RoI is to provide physically correct boundary conditions at the RoI boundaries. Thus, the mesh in these regions is coarsened and the flow there is simulated using the k − ω SST unsteady RANS model. The overall mesh structure is shown in Figure 2, where all the distances are scaled by the cylinder diameter d. The center of the coordinate system is located in the cylinder center with the z axis aligned with the cylinder principal axis, and the x and y axes oriented in the streamwise and spanwise directions, respectively.
Note the high mesh resolution directly in front of the cylinder and the grid refinement towards the test section walls, see the blue regions in Figure 2. The high-resolution mesh in front of the obstacle is necessary to correctly resolve the pressure field, which affects the flow separation. The regions close to the test section walls needed to be fully resolved in order to account for the wall effects that, for the studied flow conditions of Re = 4815, affect the wake behind the cylinder already in the region of interest. Table 1: Applied boundary conditions. The virtual walls are shown in red in Figure 1.
The used flow boundary conditions were mostly standard and are summarized in Table ??. However, it was necessary to resolve the geometric discrepancy between the experiment and simulation, shown in Figure 1. The flow field at the inlet is well defined with U in = 5ms −1 and I = 0.2 %. Still, in the simulation, the inlet needed to be moved upstream, see the red (virtual) walls in Figure 1a. To keep the flow uniformity at the contact with the cylinder, the slip boundary condition, where t w is the wall stress vector, is applied at the virtual walls. Finally, note that although the wall functions are used for the turbulence variables at the test section walls, the mesh resolution is such that they are effectively turned-off for most of the test-section and affect the flow solution only downstream from the region of interest.

Proper Orthogonal Decomposition
The proper orthogonal decomposition (POD) is a method suitable to (i) identify high-energy coherent structures in spatio-temporal data, and (ii) obtain a low-dimensional approximation of a high-dimensional process [24][25][26]. As such, it is often applied in the fields of turbulence analysis [10,27] and modeling [28,29] or in model order reduction [26,30]. The main advantages of POD are its robustness and the fact that it is a purely algebraic matrix-based method.
Both the experiment and simulation generate data through sampling of the studied system at given times, where S is a set of system snapshots and y j is the vector of data at time t j . If all the system snapshots have the same dimension, i.e. y j ∈ R m , j = 1,...,n, the set S may be rearranged into the matrix of snapshots where we assumed that the number of system snapshots n is lower than the data dimension m, which is commonly true for both experimental and numerical data. The standard POD algorithms are based on the singular value decomposition (SVD) of the matrix Y .
where E m and E n are identity matrices in R m×m and R n×n , respectively, and the first n columns of the matrix Ψ form an orthonormal basis of the column space of the matrix Y , i.e.
In (7), each column ψ j , j = 1,...,n of the matrix Ψ n represents a temporally coherent spatial structure, topos. Furthermore, at any time t j , the corresponding data vector y j may be reconstructed as a linear combination of toposes, utilizing the coefficients η j . Thus, each row j of the matrix H ∈ R n×n represents the temporal evolution of the contribution of the topos ψ j to the reconstruction of the vector y j , chronos.
Finally, let us note that given the SVD properties, the toposes, or POD modes, are sorted in descending order by the amount of variance in the original data they represent. In particular, if the data y j , j = 1,...,n correspond to the flow velocity field fluctuations u � = u − mean u, the POD modes are ordered by the amount of total kinetic energy they hold. The importance of each topos ψ j is proportional to the corresponding singular value σ j .

Computational remark:
The standard approach to POD analysis is to first compute the full matrix Ψ n and to subsequently truncate as Ψ � = [ψ 1 ,..., ψ � ] ∈ R m×� based on a cut-off energy level or error in the Frobenius norm in the Y matrix reconstruction [26,31]. However, such a task usually requires the matrix Y to be loaded into the computer memory, which may prove a problem when processing a fully three-dimensional LES simulation data [32].
Require: matrix of snapshots Y , POD rank �, number of power iterations r 1: Generate Gaussian random P ∈ R n×� 2: ComputeỸ := Y P ∈ R m×� 3: for k = 1 to r do 4: Algorithm 1: Out-of-memory randomized POD modes computation In order to limit the memory demands of POD modes computation, we employed an outof-memory randomized SVD algorithm stemming from the works [33][34][35]. Utilizing this approach, the number of modes of interest � � n = rank (Y ) needs to be known a-priori. The POD modes computation is based on a generation of a Gaussian random matrix P ∈ R n×� and summarized in Algorithm 1. Note that all the computations involving the original matrix Y , i.e. the lines 2, 4 and 7 of Algorithm 1 may be decomposed into a series of outer vector products and the matrix Y does not have to be fully loaded into the computer memory. The PIV data are shown in Figure 4, while the simulation results are depicted in Figure 6. Corresponding data in the figures are always colored utilizing the same scales. There is a significant qualitative agreement between the experiment and simulation. Both the measured and computed velocity field comprise two recirculation zones in the near-wake directly behind the cylinder and ending approximately at x = 2. Furthermore, the simulated and measured wake dimensions are roughly equal, as indicated by the yellow color in the left hand side of Figures 4 and 6. However, compared to the PIV time-averaged streamwise velocity field, the simulation estimates a slightly mode pronounced recirculation zone and a slightly bigger wake, albeit with a less apparent velocity deficiency.
A similar information as from the comparison of U may be gathered from the TKE distribution. Namely, in both the experiment and simulation, there is a region of low turbulence kinetic energy right behind the cylinder and a local maximum close to the test-section plane of symmetry at x ≈ 2.5. Still, the simulated TKE is lower than the measured one and the maximum is moved by Δx ≈ 0.2 towards the cylinder. Moreover, the computed low-TKE region behind the cylinder is smaller and not bounded by high-TKE bands, cf. the right hand sides of Figures 4 and 6 at x ≈ 1.  Finally, to evaluate the quantitative agreement between the measured and computed timeaveraged fields, we plot the streamwise velocity component U and the turbulent kinetic energy k along the line σ : x ∈ [0. 5,7], y = 0, z = 0 shown in the left hand side of Figure 4. The local measured mean velocity minimum is placed at x 1.57 and has the value of minU σ exp −0.18. The computed local velocity minimum is minU σ sim −0.33 and it is placed at x 1.61. With respect to the ranges of measured data, the relative difference between the experiment and simulation is ΔU σ ,r min = 16 % and Δx r = 0.6 %. Still, the agreement between the measured and computed U is actually the worse in minU σ . The average relative error in U on σ is 1.2 %.
The situation with the turbulence kinetic energy is similar to the streamwise velocity component. Overall, the simulation under-estimates the amount of TKE in the system by approximately 12 %. At the same time, for x ∈ [0.5, 1.5], the computed TKE seems to be actually over-estimated by roughly 27 %, which is observable even via a comparison of the qualitative results in Figures 4 and 6,respectively. Also, similarly to the case of U, the biggest discrepancy in the data is encountered for the local TKE maximum observed at x 2.28 in experiment and at x 2.22 in simulation (relative error of 0.9 %) with the respective TKE values of 0.28 and 0.24, respectively (relative error of 14 %).
Despite the quantitative discrepancies in the data, all the qualitative trends in the timeaveraged variables are in agreement. Furthermore, the quantitative discrepancies between the experiment and LES are similar or lower than the ones reported for example in [2,8] or in [19]. The initial comparison of measured and computed flow dynamics is performed utilizing the spectra of streamwise and spanwise velocity components evaluated in the points p 1 and p 2 depicted in the left hand side of Figure 4 and located at (4, 0, 0) and (4.53, 1.15, 0), respectively. The measured and computed spectra are shown in Figures 8 and 9, respectively. Similarly to [11], all the results were processed using the standard method with the Hamming window.

Turbulence spectra
In both Figure 8 and 9, we added vertical lines representing the frequencies of the three most important harmonics. All the shown vertical lines, i.e. for both the experiment and simulation, are placed at the same positions. Namely, the fundamental harmonic f 1 = 71 Hz corresponds to the Strouhal number Sr = f 1 d/U in = 0.21 ≈ 0.2, which is a typical value of Sr in the broad range of Reynolds numbers starting with Re 200, see e.g. [36] and references therein.
The higher harmonics are located at f 2 = 142 Hz and f 3 = 213 Hz, respectively. All the three harmonics may be identified in both the measured and computed spectra. Still, the measured spectra are in general more pronounced. For example, after an identical postprocessing, the maximum value of the peak V -spectrum measured in p 1 connected to f 1 is about 100, while the same maximum in the computed data is approximately 20. Furthermore, all the spectra contain the turbulent part as follows from the decrease in the inertial part of the turbulence spectrum in agreement with the Kolmogorov hypothesis (−5/3 law). However, the decrease in the measured spectra for frequencies higher than f 3 seems a little bit steeper than the −5/3 law suggests. On the other hand, at the decrease in the computed spectra in the same region follows the Kolmogorov hypothesis almost exactly.
Finally, both the simulation and experiment provide considerably different spectra in p 1 and p 2 with the spectra corresponding to the same point being qualitatively similar. For example, the peak at f 1 is missing in both the U-spectra in point p 1 and it is present in all the remaining data, i.e. in the V -spectra in p 1 and in U and V -spectra in p 2 . Furthermore, the U-spectra in point p 1 , which is located on the σ line and inside the cylinder wake show identifiable peaks at f 2 and the V -spectra at the same location demonstrate peaks at f 1 and f 3 . The spectra in p 2 , placed close to the wake border, exhibit a distinct peak at the fundamental frequency only.

High-energy POD modes
The final part of the present contribution is devoted to the POD analysis of both experimental data and simulation results. The main comparison is performed on the plane of measurement (PoM) and utilizing the streamwise and spanwise velocity components. However, the section is concluded by a presentation of fully three-dimensional POD results in the whole region of interest (RoI) as obtained from the simulation. The main idea of utilizing the POD analysis for dynamic simulation validation is based on the following. First, the topological structure of the POD modes toposes should be similar between the modes representing similar fraction of the system total kinetic energy. Second, the power spectra of chronoses corresponding to similar toposes should exhibit the same peaks or lack of them.
In both the experiment and simulation, 4000 system velocity snapshots sampled at the frequency of 2 kHz were decomposed using the POD. The experimental data were processed via a standard POD procedure implemented in the DynamicStudio software [21]. The simulation data were processed utilizing the out-of-memory randomized POD, which is summarized in Algorithm 1. Although the out-of-memory approach is not necessary for the PoM data processing, it enables the full 3D POD based on the RoI data.  The analysis was limited to the two most energetic POD modes, which combined contain more than 60 % of the system total kinetic energy for both the experimental and numerical results on PoM. In particular, the experimental modes 1 and 2 account for 35 % and 33 % of the total kinetic energy, while their simulation counterparts represent 25 % and 24 % respectively. Such a situation indicates a quasi-periodic flow, the periodic part of which is covered by the POD modes 1 and 2. The higher modes than contain the periodic flow perturbations leading to the observed chaotic turbulent flow.
The most energy-rich modes for experiment and simulation are given in Figures 10  and 11, respectively. In both the figures, the mode topos is visualized in the right hand side and the power spectrum of the mode chronos is given on the left. The topos is represented by the vorticity distribution. Comparing the toposes of the POD modes 1 based on experiment and simulation, it is possible to identify strongly similar vortex structures symmetrical with respect to the line σ . Furthermore, the local vorticity maxima and minima are placed at similar positions. Still, the computed vortices in the POD mode 1 are slightly closer to the cylinder than the corresponding experimental results, which is to be expected based on the analysis of the time-averaged flow fields.
Focusing on the chronoses of the POD mode 1, in both the measured and computed power spectra, there is a distinct peak at f 1 = 71 Hz. Thus, both the real and simulated flow fields exhibit a strong underlying periodicity partially defined by the POD mode 2 oscillating at the system fundamental harmonics. In other words, the POD-based reconstruction of the original flow field as defined by the equation (7) is a superposition method where the time-stable toposes are superposed utilizing the chronoses. Thus, a strong harmonics in a POD mode chronos suggests an almost periodic component of the trajectory of the system η(t) = Ψ T n u(t), cf. equation (7). Consequently, approximating the original variable u as u � := Ψ T � u(t), where � is the number of modes exhibiting strong harmonics in the chronos leads to an (almost) periodic flow field approximation.  The situation is almost the same with the POD mode 2. Once more, there is a qualitative agreement between the experiment and simulation results, comprising the similar topology of the modes toposes and presence of a strong fundamental harmonics in the chronoses spectra. The biggest discrepancy in the structure of toposes may be identified in the region of negative vorticity close to x = 1.8. In the experiment, there is a single local minimum placed on the line σ , while in the simulation, there are two local minima located below and above σ . A similar situation may be observed even comparing the right hand sides of Figures 4 and 6, where the experiment suggests a single maximum of the turbulence kinetic energy located on σ at x = 2.28 and the simulation result shows two distinct maxima located at a similar streamwise position, but symmetrically below and above σ . The fully three-dimensional toposes of POD modes 1 and 2 in the region of interest as obtained from the simulation data are visualized in Figure 15. The dominant vortex structures are identified via contours of the Q-criterion and colored by the vorticity component oriented in accordace with the cylinder axis (z axis in the used coordinate system), i.e. the coloring is the same as in Figures 10-13. Combined, the modes 1 and 2 form an (almost) periodic base of the chaotic turbulent flow. Moreover, the identified structures are almost uniform in the z axis direction, suggesting that the most energy-rich part of the flow is close to the 2D structure. Note that the energetic prevalence of the modes 1 and 2 usually decreases with (i) the increasing flow Reynolds number and flow chaoticity, and (ii) changes in the system geometry promoting three-dimensional flow. For example, at Re = 100 and 2D flow approximation of a flow behind a circular cylinder, the modes 1 and 2 account for more than 90 % of the flow total kinetic energy [27].

Conclusion
A validation procedure for a large eddy simulation of the flow around a circular cylinder at Re = 4815 was presented. The validation steps include comparison of (i) time-averaged flow quantities, (ii) turbulence spectra, and (iii) POD analysis results, between the experiment and simulation. It might be assumed that a simulation returning results similar to the experiment for all the three steps reflects the real flow dynamics. However, the first step, i.e. a similarity of time-averaged fields, is only a necessary condition for the success of the simulation validation. At the same time, it is the only step that is easily quantified and for which exist standard validation procedures; albeit, developed for the RANS-based turbulence modeling. As for the other steps, arguably, a model that produces similar both turbulence spectra and POD modes as the validation experiment contains similar coherent structures that are evolving in time in a similar manner and can provide a representative insight into the real flow dynamics. Unfortunately, how to evaluate the "similarity" between experimental and simulation POD analysis is still an open question. As an example, let us mention that it is not clear how many simulated POD modes need to correspond to the experimental data in order for the model to accurately catch all the important flow characteristics. In the subsequent work, we plan to extend the POD analysis towards higher modes, to include the data measured on a cross-oriented plane of measurement (i.e. in the y − z plane, see [10]) into the validation dataset, and to devise some guidelines for using POD for transient simulations validation.