On the Evaluation of Mesh Resolution for Large-Eddy Simulation of Internal Flows Using Openfoam †

: The central aim of this paper is to use OpenFOAM for the assessment of mesh resolution requirements for large-eddy simulation (LES) of ﬂows similar to the ones which occur inside the draft-tube of hydraulic turbines at off-design operating conditions. The importance of this study is related to the fact that hydraulic turbines often need to be operated over an extended range of operating conditions, which makes the investigation of ﬂuctuating stresses crucial. Scale-resolving simulation (SRS) approaches, such as LES and detached-eddy simulation (DES), have received more interests in the recent decade for understanding and mitigating unsteady operational behavior of hydro turbines. This interest is due to their ability to resolve a larger part of turbulent ﬂows. However, veriﬁcation studies in LES are very challenging, since errors in numerical discretization, but also subgrid-scale (SGS) models, are both inﬂuenced by grid resolution. A comprehensive examination of the literature shows that SRS for different operating conditions of hydraulic turbines is still quite limited and that there is no consensus on mesh resolution requirement for SRS studies. Therefore, the goal of this research is to develop a reliable framework for the validation and veriﬁcation of SRS, especially LES, so that it can be applied for the investigation of ﬂow phenomena inside hydraulic turbine draft-tube and runner at their off-design operating conditions. Two academic test cases are considered in this research, a turbulent channel ﬂow and a case of sudden expansion. The sudden expansion test case resembles the ﬂow inside the draft-tube of hydraulic turbines at part load. In this study, we concentrate on these academic test cases, but it is expected that hydraulic turbine ﬂow simulations will eventually beneﬁt from the results of the current research. The results show that two-point autocorrelation is more sensitive to mesh resolution than energy spectra. In addition, for the case of sudden expansion, the mesh resolution has a tremendous effect on the results, and, so far, we have not capture an asymptotic converging behavior in the results of Root Mean Square (RMS) of velocity ﬂuctuations and two-point autocorrelation. This case, which represents complex ﬂow behavior, needs further mesh resolution studies. a consistent trend in the results. These behaviors in the results show that further mesh reﬁnement is still required to get a converged value for the integral length scales, especially in the streamwise direction, since the difference in the results for Mesh 3 and 2 is still big.


Introduction
Hydraulic-turbines are considered a highly reliable power source that can cover an extensive range of operating conditions in response to electricity demand. The importance of the large-eddy simulation (LES) for off-design operating conditions is due to the fact that RANS (Reynolds-averaged Navier-Stokes) studies are mainly capable of accurate flow simulations near the Best Efficiency Point, since optimal swirling flow occurs and the level of turbulence is low. However, at off-design operating conditions, swirling vortices are the main phenomena inside hydraulic turbine draft-tubes and the flow is highly turbulent. Therefore, it is necessary to use high fidelity methods, such as large-eddy simulation, to capture the turbulent fluctuations in the flow. Analysis of flow phenomena for different components of hydraulic turbines has been performed by a large number of researchers in the last few decades. The URANS (Unsteady Reynolds-averaged Navier-Stokes) is the most common approach used by researchers for numerical analysis of the flow field in hydraulic turbines [1][2][3]. The limitation of URANS in predicting self-induced vortex rope was reported by Foroutan and Yavuzkurt (2012), and the importance of time-dependent boundary conditions with turbulent fluctuations at the inlet of the draft tube was shown [4]. SAS (Scale Adaptive Simulation), due to its capability to dynamically adjust the length scale in the turbulent flow and hence provide more detailed predictions of turbulent flow structures, is often employed for draft tube flow field behavior investigations [5,6]. For instance, in a study by Neto et al. (2012), SAS-SST (Shear Stress Transport) results of velocity components inside a draft tube were validated against experimental data for a Francis turbine at part load [7]. DES (Detached Eddy Simulation) as a Hybrid RANS/LES approach, applies the RANS approach for near wall modeling of turbulence and resolves the core flow region using LES. This approach was employed by Sentyabov et al. (2014) for the analysis of vortex core precession of a Kaplan turbine [8]. A comparison of URANS and DES potentials to predict turbulence statistics for the draft tube was presented by Paik et al. (2005). Although both methods agree in mean velocity field, the results of turbulence statistics show significant discrepancies for URANS [9]. DES was also used for the prediction of pressure pulsation in high-head Francis turbines by Minakov et al. (2015), and an accuracy of 10% was reported for the simulations [10].
A comparison of LES and SAS results for a sudden-expansion test case which resembles the vortex rope of a draft tube at part load was performed by Javadi and Nilsson (2014) [11]. Rotating stall mechanism of a pump-turbine was investigated by Pacot et al. (2014) using LES [12]. The results of efficiency prediction of a Kaplan turbine using Zonal LES (ZLES) is compared to SAS in a study by Morgut et al. (2015) [13]. LES with a cavitation model for Francis turbine is also performed at part load and high load in a study by Yang et al. (2016). In this study, however, the LES mesh requirements were not well-satisfied [14]. Comparison of LES results with DES and URANS results for the Francis-99 draft tube at part load, Best Efficiency Point (BEP), and high load was performed by Minakov et al. (2017) for mean velocity profiles and pressure fluctuations [15]. Compressible LES is performed by Trivedi and Dahlhaug (2018) for a whole configuration of the Francis-99 turbine to compare the high-amplitude stochastic fluctuations at speed no-load [16] and runaway [17].
In the papers featuring LES studies for hydraulic turbines, the effect of the mesh has barely been studied. Moreover, there is no consensus on the space and time resolution requirements for the LES or DES studies. The validation studies in most of the references are less than adequate, since most of the studies avoid to validate the turbulent characteristics of the flow. This points to a choice of academic test cases with proper Direct Numerical Simulation (DNS) or experimental validation data. In order to circumvent the complexities related to the geometry and other requirements for the boundary conditions of the flow simulation inside hydraulic turbine draft-tube, two academic test cases including internal non-swirling and swirling flows with available DNS and experimental data were chosen. First, a case of channel flow which has a simple geometry and includes complex near wall turbulent behavior was selected. Second, a sudden-expansion test case which resembles the swirling flow at part load operating condition of a hydraulic turbine draft-tube was investigated. In this paper, LES results are presented and analyzed for both test cases.

Simulation and Modeling of Turbulent Flow
In LES, the velocity field (U i ) is decomposed into a filtered (or resolved) component (U i ) and a residual (or sub-grid scale, SGS) component (u ). By formally applying the filtering operation to the continuity equation and Navier-Stokes equations, it is possible to derive conservation laws for the filtered flow variables. Due to the linearity of the continuity equation, applying the filtering is straightforward. The form of the equation remains unchanged.
The resulting Navier-Stokes equations are of the same form as the original Navier-Stokes equations and also a residual stress tensor or SGS stress tensor (τ r ij ) appears: The equation of motions must be closed by modeling the SGS stress tensor. This is usually performed using an eddy viscosity model. A one-equation model and the WALE (Wall-Adapting Local Eddy-viscosity) model are used in the present research to model eddy-viscosity. In the WALE model, the eddy-viscosity reads as: where and the constant is in the range 0.55 ≤ C w ≤ 0.60. A one-equation model can be used to model SGS turbulent kinetic energy: ν sgs = c k ∆k 1/2 sgs , Very close to the wall, the van Driest damping function is used to correct the behavior for the ν sgs . This function has the following form: where A + = 26. The final length scale is given by: where ∆ g is a geometric-based filter length, such as the element cube-root volume delta. Several approaches are presented in the literature to evaluate the resolution of LES. Following   [18] and   [19], some of these methods can be used by a prior RANS simulation or one LES calculation. For example, a good practice is to obtain the integral length scale L int from a prior RANS simulation as L int = k 3/2 and, therefore, to estimate the initial grid resolution for LES [20]. Therefore, they can also be called single-grid estimators, which means that, by running one LES calculation, they can be estimated. The multi-grid estimator methods require a number of LES calculations and some sort of Richardson Extrapolation. In this research, two-point autocorrelation is used for the evaluation of mesh resolutions in LES.
Correlation means the tendency of two values to change together, either in the same or opposite way. If we think of u(x) as a component of the velocity along a line in statistically homogeneous turbulence, the autocovariance or two-point correlation in space is: where r is the separation distance between two points. As the two points move closer to each other, R increases. If the points move further away, R will go to zero. R(r) is often normalized by the root-mean-square of velocity. Integral length scale L int can be computed from the autocovariance which is the integral of R(r) over the separation distance r, i.e., The two point correlations are calculated by Davidson ( , 2011 [20,21]. In this study, the separation distance is varied with the grid resolution. The number of cells needed for the autocovariance to approach zero is verified, thereby checking the grid resolution. It is recommended that in a good LES, the integral length scale should cover 8-16 cells. In this study, it is reported that the energy spectra and the ratio of the resolved turbulent kinetic energy to the total one is not a good measure of LES resolution. Two-point correlations are suggested as the best measures for estimating LES resolution.
The ways to create turbulent inlet boundary conditions vary and include the two main categories of synthesis inlets and precursor simulation methods which are briefly described in the following sections.
Divergence-free synthetic eddy method (DFSEM) was proposed by Poletto et al. (2013) [22] to reduce near-inlet pressure fluctuations and the development length in the original SEM method by Reference [23]. The steps in DFSEM are the following:

1.
User selection of inlet surface Ω.

4.
Definition of the number of eddies.

5.
Assigning random positions x k and intensities α k to all the eddies. 6.
Eddies being convected through the eddy box, u (x) calculated and superimposed to u to generate the inlet condition. 8.
Repeat steps 6-7 for all the subsequent time steps.
The fluctuating velocity field is: where: • N: the number of eddies introduced into the DFSEM; • r k = x−x k σ k with σ k being the eddy length scale for the k th eddy; • q σ (|r| k ) is a suitable shape function; and • α k i are random numbers with zero averages which represent eddy intensity. The fluctuating velocity field, taking into account the turbulence anisotropy, is: where: • d k = (r k j ) 2 ; and • βjl is the Levi-Civita symbol.
The expression for the Reynolds stresses is also expressed as follows: Therefore, σ k β and σ k γ are the eddy length scales for the k th eddy in β and γ directions. This method in comparison with original SEM methods is able to provide a divergencefree velocity field and also to reproduce all possible state of Reynolds stress anisotropy as a function of the characteristic ellipsoid eddy shapes described by the aspect ratio

Channel Flow Test Case
The turbulent flow in the channel is simulated at a moderate Reynolds number. The geometry of the channel is shown in Figure 1. The flow in the channel is considered statistically homogeneous in the spanwise direction (z) and statistically developing in the streamwise (x) and channel-height (y) directions. Domain size for the channel is (20π, 2, π) [m]. The physical modeling is done based on friction Reynolds number which is Re u τ = u τ δ ν = 395, where δ = 1 [m] is a characteristic length (Channel half-height) and friction velocity, u τ = τ w ρ is assumed to be equal to 1 [ m s ]. This research took advantage of OpenFOAM built-in mesh generators. Capabilities and meshing strategies of OpenFOAM are presented by Concli et al. (2020) [24]. The initial number of nodes for the mesh was chosen as (500, 46, 82) with a total of 1,866,000 cells. Therefore, the wall units which are the normalized distance of the first mesh vertex next to the wall are (x + , y + , z + ) = (49.36, 1.1, 15.1). This orthogonal mesh results in a maximum aspect ratio of 23. Two sets of boundary conditions are considered for the generation of the inflow turbulence, periodic and Divergence Free Synthetic Eddy Method (DFSEM). The boundaries in the spanwise direction are considered periodic. k-equation is used as the SGS model. In the initial verification study, we want to investigate the influence of the inlet eddy sizes on the development of the flow inside the channel, and we want to verify that we can obtain the mean velocity and Reynolds stress field independent of the size of eddies generated at the inlet and transported into the domain. Therefore, the effect of the size of the inlet eddies will be evaluated based on the development length of the mean streamwise velocity profile inside the channel. Our hypothesis is that further increment in the size of the inlet eddies will not further exhibit changes in the flow development length. Starting the simulation by considering 1 cell per eddy simulation, and then increasing to 3 and 5, we expect that simulating each eddy with an adequate number of 3 cells would allow capturing non-decaying turbulence fluctuations in the domain.
The results of this verification study are provided in Figures 2 and 3. In these figures, the downstream development of the profiles of the mean streamwise velocity and u u component of the Reynolds stress is presented, respectively, for several cross sections and for different values of the parameter nCellsPerEddy representing the eddy sizes at the inlet. Moreover, the results for the periodic boundary condition case is also provided in the same figure, as a reference for the comparison.
The results of these simulations show enormous changes in the profiles of the mean streamwise velocity and u u component of the Reynolds stress, which are shown in Figures 2 and 3. In these figures, we can see how the flow develops across the channel using different conditions for the inlet eddies. Changing the value of the number of mesh cells per eddy to a value of 3, which means that each eddy is represented by at least 3 mesh cells in the domain, although requiring smaller time steps, improves the results of the development of the flow in the channel. The influence of further increasing this value to 5 does not exhibit a further improvement in the flow characteristics and only imposes more expensive computations in terms of time steps. In the aforementioned figures, the results at each streamwise section is averaged in the spanwise direction and time; therefore, the notation and {YZ} index in the axis titles imply averaging in time and in the {YZ} plane, respectively.

Sudden-expansion test case
The Dellenback Abrupt Expansion test case resembles the swirling flow inside the draft-tube of hydraulic turbines at part load condition. Measurements of mean and fluctuating velocities were performed in a water flow with a laser Doppler anemometer [25]. The measurements were taken at the cross-sections shown in Figure 7. The simulation are performed at Reynolds number 30000 and swirl number 0.6. The swirl number may be physically interpreted as the ratio of axial fluxes of swirl to linear momentum divided by a characteristic radius.
The largest uncertainties in the experiments were reported to be about 2% in Reynolds number, 8% in swirl number, and 1% in probe volume positioning. Uncertainties in mean and RMS velocities stemming from the many possible biases and broadening errors were estimated to be about ±3% and ±10%, respectively. The mesh is shown in Figure 8. This mesh has 1567944 cells. This mesh corresponds to a maximum aspect ratio of 21, maximum non-orthogonality of 28, average non-orthogonality of 5 and and maximum skewness of 0.7, which are within the acceptable range. For inlet boundary conditions, velocity profiles from the measurements are specified as axisymmetric profiles using the profile1DfixedValue boundary condition in OpenFOAM. The k − ω SST model is used as the turbulence model for the RANS and WALE model is used as the SGS model for the LES. For the case of simulations at swirl number 0.6, no inlet turbulence generation method is used, as it is stated in [26], the inlet turbulence is only important when the swirl in the flow is low. The validation study is performed to compare the results of the mean velocity field and RMS of velocity fluctuations against experimental data. These results for the axial and tangential components of the mean velocity at the aforementioned cross sections are presented and compared with experimental data in In Figure 3, the cross sections with the values of x/δ = 0.06 and x/δ = 62.7 represent the center of the first cell immediately after and before the inlet and outlet boundaries. It seems that a minimum number of cells in the streamwise direction and away from the boundaries is required for the attenuation of the effect of boundary conditions. Moreover, it seems that, although both boundary conditions with 3 and 5 mesh cells are representing shorter development length in comparison with the one with one mesh cell per eddy, which is evident by the results at downstream locations of x/δ of 12.5, 50.2, and 56.5 being almost on top of each other, the results of the 5 cells per eddy is worse than with 3 cells per eddy. A possible explanation could be that there should be an optimum number of cells per eddies for a given mesh, and increasing this value does not necessarily provide more accurate results. In other words, the optimum value for this parameter could be mesh-dependent. This hypothesis needs further investigation and analysis which is beyond the scope of this research.
For mesh resolution analysis, two other meshes are generated. The details of these meshes are presented in Table 1. The results of the verification analysis for the three mesh resolutions including components of mean velocity and RMS of velocity fluctuations are presented in Figure 4. In this figure, the results are averaged in both {YZ} and {XY} planes, meaning both span-wise and stream-wise directions. The results show different flow behavior when moving farther from the wall in the outer layer (y + > 50). The finest mesh gives results that are in very good agreement with DNS across the whole region. We can conclude that grid refinement has a positive effect which is reflected in the gradual increase in the accuracy of the obtained results. The results of the RMS of velocity fluctuations show that all three components converge towards the DNS data as the mesh gets refined. Closer to the core region of the channel (y/δ > 0.5), mesh 2 shows generally better prediction for the three components of RMS of velocities. The verification analysis is performed to investigate the impact of the mesh resolution on the LES results. To this aim, two other meshes one coarser and one finer than the one for which results were presented in the previous paragraph were generated. The details about these meshes are given in Table 3. The results of these simulations are provided in Figure 13 and Figure 14 for the mean velocities for axial and tangential components and in Figure 15 and Figure 16 for RMS of velocity fluctuations, again for axial and tangential components. These results show higher sensitivity of the velocity profiles to mesh resolution for both mean and RMS values. A higher mesh resolution does not necessarily give less deviation from experimental data, due to the complex nature of the flow, especially for RMS of velocity fluctuations. However, in general, results for Mesh 3 appear to be slightly more consistent with experimental data. The results of the two-point correlations of the components of the velocity fluctuations at the centerline of the cylinder downstream of the expansion are calculated for the three meshes. These results are presented in Figure 17. This figure shows that the complex vortical structures also represent a complex interactions which does not linearly go to zero as it changes randomly further downstream of the expansion.  The second verification study for the mesh resolution is done for the two-point autocorrelation of velocity fluctuations. The results of this analysis are presented in Figure 5 for y + = 5 . It should be noted that the results of autocorrelation in the x direction are averaged in a cross section plane (z-plane) and the results of autocorrelation in the z direction are averaged in the stream-wise plane (x-plane). In general, the profiles indicate that the chosen size of the domain is large enough to accommodate all the relevant turbulent structures, since the profiles reach near zero value or an asymptotic value in each computed direction. These profiles show higher sensitivity of this indicator to the mesh resolution, especially in the x-direction. The results computed with DFSEM boundary conditions, overall, show better performance than the periodic case, almost for all the cases, even for the coarsest mesh. In the integral length scale results in the Table 2, we do not observe a consistent trend in the results. These behaviors in the results show that further mesh refinement is still required to get a converged value for the integral length scales, especially in the streamwise direction, since the difference in the results for Mesh 3 and 2 is still big.   The third verification study for the mesh resolution is done based on the energy spectra. The objective of this verification analysis is to evaluate the influence of the mesh resolution on the energy spectra and to compare the spectra in the inertial subrange with the Kolmogorov spectrum. The results of this verification analysis are presented in Figure 6. These results show a much less sensitive behavior of the energy spectra to the mesh resolution; therefore, this parameter does not prove to be useful as an indicator of the LES mesh resolution. two-point autocorrelation and hence the integral length scales showed that even for a simple case of channel flow, obtaining mesh-independent results is still challenging and further mesh-dependency analysis is required. The verified framework for LES on the channel flow was then applied to study the flow in a sudden expansion pipe flow at a Reynolds number of 30000 and swirl number of 0.6. As with this case, the validation analysis of the LES WALE results versus RANS k-ω SST results show a large improvement of the LES results over the RANS ones for the mean and RMS of axial and tangential velocities. These improvements were also more evident at the cross sections which were immediately after the expansion. The influence of mesh resolution was also studied and showed a higher sensitivity of the results of this case to the mesh resolution and characteristics. The energy spectrum of the pressure fluctuations also did not exhibit much sensitivity to the mesh resolution, as it was also evident in the channel flow test case.

Acknowledgments
This research was supported by the Natural Sciences and Engineering research Council of Canada through Discovery Grants

Sudden-Expansion Test Case
The Dellenback Abrupt Expansion test case resembles the swirling flow inside the draft-tube of hydraulic turbines at part load condition. Measurements of mean and fluctuating velocities were performed in a water flow with a laser Doppler anemometer [25]. The measurements were taken at the cross-sections shown in Figure 7. The simulation are performed at Reynolds number 30,000 and swirl number 0.6. The swirl number may be physically interpreted as the ratio of axial fluxes of swirl to linear momentum divided by a characteristic radius.
The largest uncertainties in the experiments were reported to be about 2% in Reynolds number, 8% in swirl number, and 1% in probe volume positioning. Uncertainties in mean and RMS velocities stemming from the many possible biases and broadening errors were estimated to be about ±3% and ±10%, respectively.
The mesh is shown in Figure 8. This mesh has 1,567,944 cells. This mesh corresponds to a maximum aspect ratio of 21, maximum non-orthogonality of 28, average non-orthogonality of 5, and and maximum skewness of 0.7, which are within the acceptable range. For inlet boundary conditions, velocity profiles from the measurements are specified as axisymmetric profiles using the profile1DfixedValue boundary condition in OpenFOAM. The k − ω SST model is used as the turbulence model for the RANS and WALE model is used as the SGS model for the LES. For the case of simulations at swirl number 0.6, no inlet turbulence generation method is used, as it is stated in Reference [26], and the inlet turbulence is only important when the swirl in the flow is low.  The verification analysis is performed to investigate the impact of the mesh resolution on the LES results. To this aim, two other meshes one coarser and one finer than the one for which results were presented in the previous paragraph were generated. The details about these meshes are given in Table 3. The results of these simulations are provided in Figures 13 and 14 for the mean velocities for axial and tangential components and in Figures 15 and 16 for RMS of velocity fluctuations, again for axial and tangential components. These results show higher sensitivity of the velocity profiles to mesh resolution for both mean and RMS values. A higher mesh resolution does not necessarily give less deviation from experimental data, due to the complex nature of the flow, especially for RMS of velocity fluctuations. However, in general, results for Mesh 3 appear to be slightly more consistent with experimental data.          The third verification analysis of the mesh resolution is performed for the results of the energy spectra of pressure fluctuations. The objective of this analysis, as for the case for the channel flow, is to investigate the effect of the mesh resolution on the convergence behavior of the energy spectra and also to compare the spectra in the inertial subrange with the Kolmogorov spectrum. Time history of the pressure fluctuations at several points at the wall is recorded and monitored during the simulations. Power spectral density of the pressure fluctuations at one point at z/D = 2 is calculated for the three meshes, and results are presented in Figure 18. This case also shows less sensitivity of the results to the mesh resolution, and this criteria can again not be used as an indicator of the mesh resolution adequacy.

Conclusions
This work is dedicated to developing an expertise on how to use, validate, and verify the OpenFOAM CFD code for the large-eddy simulation of the flow types similar to the ones which occur at off-design operating condition, such as part load inside hydraulic turbines. Two main test cases, a case of turbulent channel flow and a case of sudden-expansion, were considered in this study. The first one brings the classical problem of near-wall turbulence complexities, and the second one resembles the swirling flow inside a hydraulic turbine draft-tube at part load off-design operating condition. The calculated profiles of two-point autocorrelation and hence the integral length scales showed that, even for a simple case of channel flow, obtaining mesh-independent results is still challenging, and further meshdependency analysis is required. The verified framework for LES on the channel flow was then applied to study the flow in a sudden expansion pipe flow at a Reynolds number of 30,000 and swirl number of 0.6. As with this case, the validation analysis of the LES WALE results versus RANS k-ω SST results show a large improvement of the LES results over the RANS ones for the mean and RMS of axial and tangential velocities. These improvements were also more evident at the cross sections which were immediately after the expansion. The influence of mesh resolution was also studied and showed a higher sensitivity of the results of this case to the mesh resolution and characteristics. The energy spectrum of the pressure fluctuations also did not exhibit much sensitivity to the mesh resolution, as it was also evident in the channel flow test case.