A dynamic mixed subgrid-scale model for large eddy simulation on unstructured grids: application to turbulent pipe flows

The paper presents a consistent large eddy simulation (LES) framework which is particularly suited for implicitly filtered LES with unstructured finite volume (FV) codes. From the analysis of the subgrid-scale (SGS) stress tensor arising in this new LES formulation, a novel form of scale-similar SGS model is proposed and combined with a classical eddy viscosity term. The constants in the resulting mixed model are then computed trough a new, cheaper, dynamic procedure based on a consistent redefinition of the Germano identity within the new LES framework. The dynamic mixed model is implemented in a commercial, unstructured, finite volume solver and numerical tests are performed on the turbulent pipe flow at Reτ = 320−1142, showing the flexibility and improvements of the approach over classical modeling strategies. Some limitations of the proposed implementation are also highlighted.


Introduction
In the last two decades, LES has emerged as the working horse of turbulence modeling approaches for several classes of flows and increasingly tends to be used for complex multi-physics applications [1]. However, while the numerical techniques and practices have strongly evolved in the years, recent LES approaches still tend to be strongly tied to the earliest formalization of the field [2]. This apparently innocuous trend can, nonetheless, become a limit when the so called implicit filtering approach is adopted. In implicitly filtered LES the separation between large, resolved scales and small, unresolved ones is implicitly determined by the grid and the numerical method in use; as a consequence, the features of the resolved flow field and the requirements for the necessary SGS model are directly affected.
A well known example of this comes from the implicit filtering approach with the FV method, where the flow equations are affected by the top-hat filtering implied by the integral formulation and a closure problem arises at the flux level. Nonetheless, this is not usually considered and most of the LES formulations for FV usually rely on the classical SGS closures originally developed in a continuous differential framework. Hence, as the SGS model is the only one responsible for representing the filter effects, a basic contradiction arises for implicitly filtered LES with FV: not only the discretization error is neglected but, more importantly, the form of the underlying implicit filter is assumed as completely different from the real one [3].
The aim of the present work is to introduce an LES framework which is a generalization of previous works [4,5] and allows a clear, consistent identification of filter effects which is independent from the specific numerical method adopted. The framework is then applied to the case of an implicitly filtered FV discretization and a Taylor series analysis of the resulting SGS stress tensor is used to derive a form of scale-similar model which is fully consistent with the numerical formulation. Finally, the scale-similar term is combined with a classical eddy viscosity closure in a new form of dynamic procedure which is substantially cheaper than the classical one. In extension of preliminary results obtained for the turbulent plane channel flow [6], the formulation is tested for the first time with unstructured hexaedral grids in the simulation of the turbulent pipe flow at Re τ = 320 − 1142 [7,8] using the commercial FV solver Fluent.

LES formulation
As in classical LES, the filtering operator G i is introduced: however, differently from classical LES, here the fluid domain Ω is considered finite and the filter kernel G i has a spatially varying width ∆ i (⃗ x). Notice that this is not a convolution filter, hence commutation with spatial derivative operators does not hold; however, the interaction with the domain boundaries can be avoided (as assumed here) by a proper reduction of the filter width. For future reference, we also introduce the conventions: To make the approach as general as possible, the classical Favre filtering is also introduced, ρ being the density of the fluid. With this notation, it is possible to express the Navier-Stokes equations at a generic filter level n ≥ 0; limiting the discussion to the momentum equations, we get: where τ n−0 ij is the SGS stress tensor arising in this new formulation. With respect to classical LES formulations, we notice that no commutation with spatial derivatives is ever required and, as a consequence, the resulting SGS stress tensor is completely different as it now also includes pressure and viscous terms. A classical criticism with the proposed formulation is the lack of Galilean invariance; however, this is a property of any quantity which is non-uniformly filtered in space, hence it is natural that equations based on such quantities inherit this property too [9].

Finite volume based LES
The proposed formulation has a clear advantage in the fact that the filter (1) can be readily interpreted, without any additional assumption, as a FV discretization of the equations:  3 being the volume of the generic grid cell. As the resulting equations now exactly fit those usually employed in general FV codes, two important consequences follow: a) FV codes are inherently based on fully consistent (i.e., without commutation errors) LES equations based on a top-hat filter and do not require specific modifications besides SGS modeling; b) the SGS stress tensor arising in the FV based LES approach is substantially different from the classical SGS stress tensor and, in theory, a dedicated modeling is required. With respect to point (a), notice that it holds not only for the present density-based formulation, used just for the sake of generality, but also for pressure-based ones; indeed, it is a matter of simple manipulations to show that the following exact pressure equation holds: which, once interpreted according to (4), exactly represents the usual FV discretization of the pressure equation (possibly assuming a constant density).
With respect to point (b), which is fully developed in the next sections, here we want to remark its importance in the implicitly filtered approach, as the SGS model becomes the only one representing the underlying filter operator and the resulting LES approach as a whole. However, even if this is not considered and a classical SGS model is used, an important difference still arises in the Germano Identity (GI) which holds in the present formulation and the resulting dynamic procedure [10]. Indeed, for any (n, m, k) combination of filter levels, the following generalized GI holds among the underlying SGS stresses: (6) which presents several differences with respect to the classical one usually adopted. First, it does not involve any test filtered tensor, which means that the resulting dynamic procedure does not require the constants to be arbitrarily extracted from the test filter operation; also, this implies a substantial saving in computational costs as only the single variables need to be test filtered (including pressure and, possibly the density), in contrast to the classical dynamic procedure which requires filtering products of variables and additional non-linear terms as present in the adopted SGS model. As an example, for incompressible flows, the computation of the constant in a dynamic Smagorinsky model would only require 3 applications of the test filter with the present method against the 15 of the classical dynamic procedure (the 3 velocity components plus the 6 independent components of two additional tensors). Second, due to the peculiarity of the SGS stresses, the present GI also involves the pressure and viscous terms. The saving in the computational costs of this new GI will be exploited by a dynamic, two parameters, mixed model whose scale-similar part is developed in the next section.

SGS tensor analysis
While the appearance of the SGS stress tensor in a different form would generally require a different form of SGS model, this does not seem to be a relevant issue when an eddy viscosity parameterization is used. Indeed, as long as the required role of the model is an added dissipation, its specific implementation can be considered of secondary significance with respect to the main modeling assumption. However, this reasoning loses validity when a structural modeling approach is adopted, like with scale-similarity models [11]. Indeed, structural SGS models make no assumptions on the nature of the interaction between resolved and non resolved scales but try to recover the full structure of the SGS stress tensor (i.e., the eigenvectors) from information in the resolved scales; this, in turn, allows the possibility for backscatter and anisotropy which, for example, are particularly important in wall-bounded turbulent flows. Among the structural SGS models, scale-similar ones are peculiar in the fact that they extract the required information from the smallest resolved scales trough repeated explicit filtering, as these scales are the most involved in the interaction with the sub-grid ones; also, according to a scale-similarity hypothesis, the SGS stress tensor at the basic filter level is approximated with an equivalent tensor based on the smallest resolved scales. Hence, in order for scale-similar models to be representative of the actual SGS stress tensor, a correct definition is required for the true SGS stress tensor and the underlying filter determining it. In order to leverage the several benefits offered by structural models with respect to simple eddy viscosity closures, in the present work we want to introduce a scale-similar SGS model which is suitable for the present FV based LES formulation, and take advantage of the reduced computational costs of the new GI by using it in a two parameters, dynamic mixed SGS model. In order to understand what is the proper form of a candidate scale-similar term, we analyze the Taylor series development for τ m−n ij assuming that, for m > n, the following hold [13]: where M k are filter dependent coefficients and ∆ m is the filter width resulting from the consecutive application of the filters G i for i = n, n + 1, . . . , m. Under this circumstance, which is especially relevant for a FV discretization, the substitution of the expansion (7) in the definition for τ m−n ij in (3) and additional manipulations (omitted for the sake of conciseness) lead to the following estimates: The resulting approximation for τ m−n ij can then be summarized as follows: where, following the product differentiation rule, f and h are functional forms due, respectively, to the constant and variable ∆ 2 m parts in the estimates above (h being due just to the viscous part in (10)). However, when the same analysis is performed for the true SGS stress tensor τ n−0 ij for sufficiently smooth variations of ∆ n , the fact thatφ m =φ n + O accuracy order) a substitution of the n-level variables with the corresponding m-level variables in the right hand side of equations (8-10); hence, the resulting estimate for τ n−0 ij is: The main conclusion to be drawn from this analysis is that a scale-similar model based on τ m−n ij , as originally proposed in [4], has a second order difference with respect to the true SGS term. In order to overcome such deficiency, we propose to scale the model according to the following estimate: From which it follows that a scale-similar tensor of the form accuracy; if, in addition, the ratio ∆ 2 n /∆ 2 m is constant in space, the accuracy is restored also for the diffusive part. This form of scale-similar term is thus a better candidate to represent the implicit filtering effects due to the FV formulation. However, it is worth noting that numerical error effects, due to the discretization of the fluxes on the faces of the FV cells, are not considered here; also, the relation (13) says nothing about its applicability and the relative importance of the neglected higher order terms. From the point of view of the practical implementation in general FV solvers, the n th filter level quantities in the model are just those already available at the basic implicit filter level determined by the FV discretization (i.e., n = 1 and ∆ n = ∆ 1 ). In contrast, the m th filter level quantities require explicit test filtering according to the assumption (7); while this form of differential filters is relatively easy to implement in unstructured solvers, both in explicit or implicit form, in the present work we adopt a simpler volume average filter over an extended neighbor stencil (see figure 1) in order to save computational time. While the present model is strictly dependent on the selected test filter, this choice can be justified by the fact that the model only considers the implicit filtering effects due to the FV formulation and not those due to the discretization, hence an additional approximation would still be present; also, as the model is implemented in a dynamic version, it is reasonable to expect a partial compensation effect for these errors due to the dynamic procedure. The limitations arising from this choice will be partially assessed in the following while more consistent options for the test filter will be the object of future works.

Dynamic mixed SGS modeling
Taking advantage of the reduced costs of the dynamic procedure within the FV framework, the new scale-similar term is combined with an eddy viscosity closure in a two parameter, dynamic mixed model. Introducing the test filter levels m and r (∆ r > ∆ m > ∆ n ), the proposed parameterization for the SGS stresses at the basic (n) and test (m) filter levels reads: where C ev and C ss are the dynamic constants to be computed and: If the relations (14) are expressed in compact form as: and inserted in the new GI (6), the constants in the model can be computed by a least-squares error minimization [10,11], leading to: . According to the considerations made in section 4, a standard mixed-scale model [14] is selected as eddy-viscosity model, the parameter θ being free to vary in order to obtain different closures. When compared to classical, two parameters, dynamic mixed SGS models [11], the present approach still provides a substantial saving in computational costs as, for the incompressible case, only 6 scalar fields need to be test filtered (3 velocity components per each test filter level involved, the isotropic pressure term being usually excluded from the modeling) against the 24 of a classical approach with a scale-similar model (3 velocity components per test filter level, plus the 6 components of the Leonard tensor per test filter level and the 6 components of the SGS model tensor at the basic filter level).

Flow description and computational setup
The test case considered is represented by the fully developed, turbulent, incompressible, isothermal flow in a straight circular pipe of radius R. The mean flow direction is along the z axis, while the r and θ directions define, respectively, the radial and tangential coordinates in planes normal to the z direction. The velocities in the three coordinate directions are, respectively, u z , u r and u θ . Two distinct cases at different Reynolds numbers are considered, respectively Re D = 10, 000 and Re D = 44, 000 based on the mean velocity U b and pipe diameter D. The corresponding friction Reynolds numbers, based on the friction velocity u τ and pipe radius R, are Re τ = 320 and Re τ = 1142. In both cases, a no-slip boundary condition is applied at the wall and a periodicity condition along the axial direction; this latter option, when applied to overly constrained domain (i.e., pipe) lengths deserves further attention. Indeed, in recent years, several DNS studies have put considerable attention to the description of large and very large-scale motions in turbulent pipe flows [15,16,17] in order to better clarify their role, as energy carrying structures, in the determination of higher-order flow statistics. Hence, when periodic boundary conditions are applied in a computational study, the issue of the minimum required pipe length to avoid numerical artifacts due to artificial boundary conditions clearly emerges. Previous experimental studies [18,19] evidenced how such large scale structures can have wavelengths as long as 16R in the outer region of the flow and this has led some DNS studies, e.g., [8], to adopt a pipe length L z = 15R. However, more recent computational works [20] have identified the more stringent requirement L z = 8πR in order for higher order statistics to converge within 1% of the respective values for longer domains (up to L z = 20πR); as a matter of fact, the current trend in DNS of turbulent pipe flows [15,16,17] is to adopt L z = 30R. Despite this evidence, in the present work we have adopted the axial domain lengths L z = 10R for the low Reynolds case [7] and L z = 15R for the high Reynolds one [8], both in accordance with the respective DNS references. This was motivated by two main reasons. The first, more practical one, is that there is abundant empirical evidence [7,8,16,17,20] of the fact that, for the low order statistics which are analyzed here, no substantial differences arise between the present reference DNS and computations on longer domains or experimental data; hence, the computational cost reduction resulting from a shorter domain length has been readily exploited. The second, more theoretical, reason stems from the fact that the selection of a reference DNS case allows to treat the comparison in purely mathematical terms, without any further reference to the relative target flow, as repeatedly stated in [21](pp. 111 and following). In the present case, this means that the proposed LES approach is tested toward the fully resolved periodic turbulent flow in a pipe with length, respectively, L z = 10R for the low Reynolds case and L z = 15R for the higher one; that is, the reference mathematical problem, the one which is actually discretized in the computations, exactly has the same boundary conditions for the present LES and the reference DNS. Also, these are exactly the unique boundary conditions which would allow the present implicitly filtered LES to converge toward the reference DNS with further grid refinements. A possibility, which would still not escape this reasoning, is that the large scale motions which cannot be captured by the present domain sizes (also known as super-grid scales) have such a strong interaction with the modeled sub-grid scales to require a specific SGS modeling approach. While this would invalidate the LES approach as a whole and, to the best of authors' knowledge, has never been investigated before 1 , here we simply assume this to be a possible deficiency of the SGS model which, in case, will be the object of future works.  For what concerns the computational grids, as the test case is selected to be also representative of more complex applications, unstructured hexaedral grids of different resolutions are used, the topology being sketched in figure 2. As summarized in table 1, a single grid resolution is used for 1 Such a test would require the comparison of LES computations made on domains of different lengths with the relative DNS data from computations on the same domains. the low Reynolds case, while three different resolutions are tested for the higher Reynolds case. Computational parameters for the reference DNS computations are also reported to highlight the difference, with respect to the present computations, in terms of total degrees of freedom represented. The computations are performed with the commercial, unstructured, FV solver Fluent [22] and the following numerical parameters are common among all the simulations: a non-iterative time advancement with the fractional step algorithm for the pressure-velocity coupling; an implicit second order algorithm for time advancement with a time step selected to maintain the maximum Courant number around 0.1; non-dissipative, second order, central schemes for the spatial discretization of convective and diffusive terms. According to equation (5), no specific modifications are needed for the pressure-velocity coupling, however the pressure is not included in the scale-similar part of the SGS model. The dynamic mixed model (DM) is tested with θ = 0 and compared with a dynamic Smagorinsky model (DS, C ss = 0 and θ = 1) and the reference DNS data [7,8]. In addition, computations with a classical dynamic Smagorinsky model [10] have been performed, but no appreciable differences with the present implementation were found and are not presented; nonetheless, this confirms that the new dynamic procedure introduced here can safely replace the classical one while attaining a significant reduction in the computational cost of the model. It is also worth mentioning that no spatial average of the dynamic constants is ever performed as it is not explicitly required by the formulation and in no case stability issues emerged during the computation. However, a clipping of the dynamically computed constants was adopted in order to limit the detrimental effects due to the approximate least square procedure; in particular, according to previous results [6], the constants were allowed to vary within the limits 0 ≤ C ev ≤ 0.18 2θ 0.2 1−θ and 0 ≤ C ss ≤ 2.
The flow is initialized with a power law velocity profile with superimposed random fluctuations. After the statistically steady state is reached, flow samples are collected at each time step for a total sampling time of 80R/u τ or, approximately, 1200R/U b and 1540R/U b for the low and high Reynolds cases; enough to allow a particle to travel more than 100 times trough the pipe axial length at velocity U b . Additional spatial averages are performed in the homogeneous direction θ over 4 perpendicular radial profile samples. All the results are made non-dimensional by the friction velocity u τ and the pipe radius R.

Numerical results
Flow statistics for the Re τ = 320 case are presented first, in figure 3. As both mean and fluctuation profiles show, no substantial differences are found between the proposed DM model and the simpler DS formulation. This could be explained by two main reasons: the first one is the low Reynolds number of the case, which implies a lack of inertial range structures in the flow and a relative low significance for the scale-similar part of the model; the second one is the relatively high resolution of the grid for this case (see table 1), which makes the contribution of the scale-similar term also unimportant with respect to the resolved part of the flow. However, both models provide a fairly good agreement with the reference DNS data.
Differences among the DM and DS modeling strategies are better highlighted in the Re τ = 1142 case, whose main statistics are presented in figure 4. This is especially true for the mean axial velocity and the axial velocity fluctuations. Indeed, while none of the models can correctly reproduce the mean DNS velocity profile, the proposed DM model provides a better agreement for both the medium and fine grids (the relative curves being completely overlapped), the error on the resulting Re D being 8% (it rises to 11% for both the DS Medium and DM coarse cases); the matching with the DNS data becomes almost perfect for the axial velocity fluctuations, as the DM model can substantially reduce the over-prediction of the inner-layer peak on all the grids. Despite some residual effects due to the convergence of the statistics, no particular differences are found on the remaining velocity fluctuations profiles except, maybe, for the tangential velocity fluctuations peak, which is better predicted by the DS model. What is worth noting is that the DM model can provide the same level of accuracy of the DS model at a fraction of its cost (the coarse grid employed here allowed a 30% reduction in computational costs with respect to the medium grid employed for the DS model).
However, these results also highlight a feature of the proposed model DM which, especially for the mean velocity profile, might appear clumsy, at least, for the present implicitly filtered code: a sort of grid convergence toward a non DNS result. A possible explanation is connected to the non optimality of the adopted grid refinement with respect to an expected grid convergence toward the DNS result; indeed, not only the three grids all have tangential and stream-wise spacings well above the Taylor microscale and within the inertial range (and are thus non-optimal [23] and expected to provide a sort of inertial range convergence [24]), but the near wall radial grid spacing has been kept constant among the three grids, possibly limiting the effectiveness of the adopted grid refinement. A second, less evident, possible explanation is connected to the clipping limits used for the constant C ss and the resulting effectiveness of the scale-similar part of the model. Indeed, the present clipping limits have been determined according to previous results for the turbulent channel flow [6] where, as shown in the left part of figure 5, allowed the full dynamic range for the constant C ss in the first layer of cells near the walls, where the scalesimilar model is more effective; in contrast, for the present turbulent pipe flow at Re τ = 1142, as shown in the right part of the same figure for the coarse grid (the same results hold also for the two additional grids not reported), the effect of the clipping has strongly limited the dynamic range for the constant near the wall. In particular, it is worth noting that the dynamic range of the unclipped constant for the turbulent pipe is much larger than the interval represented in figure 5 but here it has been restricted for a better graphical representation; for the same reason, the clipped constant for the same case has two strong peaks at 2 and 0 which have been removed from the picture. Hence, this behavior suggests that, most of the time, the scale-similar model is not active at the wall or not as much as would have been required. As the cut in the positive tail of the distribution has been observed also for the lower Reynolds case (not shown), without the same negative effects, it appears probable that the major problem of the present approach for the pipe resides in the clipping of negative constants, which are instead present only for the higher Reynolds case. However, none of these explanation has been completely confirmed yet and further investigations are required.
The nature of the two modeling strategies, DM and DS, can be better appreciated by  analyzing the stream-wise spectra of the underlying velocity fields, presented in figure (6) for the common medium grid. For both the velocity components, but especially for the radial one, the DM model provides a considerable recovery of energy at the smallest represented scales with respect to the DS model; however, the model does not seem capable of totally removing the energy pile-up which is evident for the radial velocity component of both the tested models. This fact can be explained by considering the nature of the test filter and its previously cited limitations. Indeed, the extended filter stencil adopted here is actually required in order to obtain a larger filter width and make the dynamic procedure sufficiently independent from the discretization error which, for the computational grids in use, tends to pollute a much wider range of scales with respect to more regular grids 2 . However, for the top-hat like filter used here, the use of a larger filter width also implies an incomplete removal of the energy content at the smallest resolved scales, as indeed appears from the velocity spectra. Nonetheless, the same figure also shows that the test filter applied on the resolved velocity field for the DM model can reproduce with a certain accuracy the velocity spectra of the DS model; this heuristically suggests that, in relation to the DM model used here, the present test filter can sufficiently represent the implicit filtering effect due to the FV formulation.  Figure 6. Stream-wise spectra at r = 0 for the Re τ = 1142 Medium case. Left: axial velocity, right: radial velocity. Legend: --DM Basic (n), · · · · · ·DM Test (n + 1), ----DS Basic (n).