Assessment of a Flow-dependent Subgrid Characteristic Length for Large-Eddy Simulation on Anisotropic Grids

This paper presents the latest results of a long track development activity in the context of low-dissipative finite volume method for compressible flows. Specifically, here we focus our attention on the Large-Eddy Simulation (LES) approach which can be considered a good candidate for turbulent flow simulations over the next decades. One of the key ingredients of LES models is the subgrid length scale which is typically evaluated based on the local mesh size. This standard approach suffers from loss of accuracy on anisotropic grids that are commonly employed to obtain sufficient wall-normal resolution, whilst keeping the total cell count to a minimum. In order to avoid this issue, we assess the effectiveness of a velocity-gradient-based length scale, referred to as least square length (LSQ) [1]. In this paper, we present for the first time results obtained with the LSQ length scale in the context of compressible LES. The superiority of the LSQ approach over the standard cubic-root length scale is demonstrated in terms of accuracy and overall time to solution.


Introduction
Direct numerical simulations (DNS) of transitional/turbulent flows are still not achievable for practical applications due to the required computational effort.For this reason, the Large-Eddy Simulation (LES) approach is commonly adopted as the best compromise between solution accuracy and required resources.LES models require the evaluation of a sub-grid length scale which is commonly related to the local mesh size.Standard, geometry-based length scales can strongly affect the accuracy of the solution when anisotropic grids are employed [1].This issue can be solved by increasing the overall number of cells, with the aim of reducing the cell anisotropy.However, this incurs a significant increase in the required computational resources.In order to propose a more affordable solution, Trias et al. [1] proposed a velocity-gradient dependent length scale, termed least square length (LSQ).The LSQ subgrid length scale, was demonstrated to minimize the effects of cell anisotropy on the accuracy of LES of incompressible flows.In this work, the effectiveness of the above described filter length was assessed in a compressible solver, for the well-established test case concerning the flow over a circular cylinder operating at a Reynolds number equal to 3900.This test case was selected thanks to the vast amount of available experimental and numerical results.The previously published results put in evidence that there is a flow bifurcation in the very near wake region.So, it can be considered a reliable test for the sake of benchmarking the proposed numerical simulation technique.This paper is organized as follows: in Sec. 2 the governing equations are presented, Sec. 3 provides a brief description of the discretization technique, Sec. 4 contains the numerical results and Sec. 5 summarizes the conclusions.

Governing equations
The flow model considered in this work is compressible LES in which continuity, momentum conservation and energy conservation equations are solved as follows: where the overbar denotes the Reynolds averaging the tilde is the Favre averaging.In eq. 1 ρ denotes the density, ũ is the velocity vector, p is the pressure and Ẽ denotes the total energy which is related to the other flow variables in the following fashion: In the implementation here presented, the last term in eq. 2 is not included in the solver since it produces a negligible impact on the obtained solution at low Mach numbers, [2,3].On the other hand, the closure of eq. 1 is achieved through the constitutive equation for Newtonian fluids, Fourier law and the equation of state for the perfect gas.Lastly, the ĖSGS contribution in eq. 1 is handled in the following way, [2]: since sub-grid scales (SGS) turbulent diffusion, and the SGS viscous diffusion can be considered negligible, [2].The anisotropic part of the SGS stresses, τ SGS , is parametrized using a standard approach: where S is the strain rate tensor and the SGS viscosity, µ t , is obtained using the following equation: In eq. 5 C s is the so-called Smagorinsky constant, ∆ is the sub-grid length scale and S is the magnitude of the strain rate tensor.By contrast, the dynamic procedure proposed by Germano-Lilly, [4], was adopted to compute C s .As regards to the sub-grid length scale, we have employed two different approaches.The first approach is based on the cube root of the cell volume: While the second approach is the least square length (LSQ), proposed in [1]: 40th UIT International Heat Transfer Conference (UIT 2023) Journal of Physics: Conference Series 2685 (2024) 012008 In the above equations V P is the cell volume.Differently, LSQ sub-grid scale is built starting from the velocity gradient tensor G and the length scales tensor ∆ D = (∆x, ∆y, ∆z) I.It is also very important to put in evidence that the length scales ∆x, ∆y, ∆z appearing in ∆ D should be the same as the discrete lengths used to obtain the discrete velocity gradient [5].Using Green's theorem, the generic cell P velocity gradient components are calculated as follows: with n x j P N = [n x , n y , n z ] T P N denoting the outward-pointing unit normal vector to the interface P N which is shared by the cell P and its neighbour cell N , The term S P N denotes the face area of the interface P N and N (P ) is the set of cell P 's direct neighbors.Lastly, the corresponding cell length scales, ∆x j = [∆x, ∆y, ∆z] are given by:

Numerical approximation
A standard, cell-centered finite volume method, available within the OpenFOAM framework, is adopted to discretize the governing equations described in Sec. 2. Convective terms are evaluated using Pirozzoli's second-order accurate, energy-conserving scheme [6], while a standard central-differencing scheme is adopted for the diffusive contributions.
For time integration, a compact-storage (2N) explicit Runge-Kutta (ERK) scheme was implemented thanks to its benefits for large scale computations.In particular we have considered a five-stage, fourth-order accurate, ERK scheme available in the open literature, [7].Our computational experience showed that this scheme is extremely appealing since it allows to use a maximum Courant number equal to 1. Thus, it allows to exploit the good parallel performance of the explicit approach, without imposing a severe time-step restriction.
All the aforementioned numerical methods were implemented in a freely-available OpenFOAMbased solver, termed caafoam, which is extensively described in D'Alessandro et al. [8]

Results
The test case considered in this work concerns the flow past a circular cylinder operating at a Reynolds number, Re, based on the cylinder diameter D, equal to 3900 [9, 10], and freestream Mach number equal to 0.2.This is a well-established benchmark for evaluating the simulation capabilities of separated and transitional flows.The flow is characterized by a separated shear layer that detaches from the cylinder surface in a laminar regime and becomes turbulent downstream of the bluff body, as depicted in Fig. 1.Moreover, a wide array of experimental and numerical results put in evidence that there is a flow bifurcation in the very near wake region: one flow state is characterized by the presence of a V-shaped profile of the stream-wise velocity [11], while a second state reveals the presence of a U-shaped profile, [12].Both states have been recognized in many LES and DNS computations.Recent literature studies, [13], confirmed that the U-shaped state is correct and that the V-shaped profile occurs both in experimental analyses and in numerical simulations because of an early laminar-to-turbulent transition of the separated shear layer.It is important to note that obtaining the correct solution of a U-shaped profile is not a trivial task.This specific flow problem has consequently become a popular case for testing LES solvers and associated numerical schemes.In this work the computational domain extends for 38 D in the radial direction starting from the  cylinder centre and π D in the spanwise direction.A set of fully structured O-type grids were generated using N θ = 300 cells in the circumferential direction, N r = 330 cells are employed in the radial one, while the span was discretized using a varying number of cells, N z , ranging from 32 to 64, as summarized in Table 1.Computational cells were clustered near the cylinder surface and the dimensionless height of the first cell next to the wall was imposed to be 10 −3 .An artificial sponge layer was also adopted as non-reflecting boundary treatment to damp the solution to a reference value at the far-field, avoiding unphysical wave reflections.In the considered configuration, the sponge layer starts from the external boundary of the computational domain; its thickness, L sp , is fixed to 13 D, leading to a dimensionless extension, L sp f /c ∞ ≃ 0.5, where f is the vortices detachment frequency and c ∞ is the free-stream speed of sound.Moreover, as suggested in Mani et al. [14], a sponge strength η target of 40 dB was chosen.Finally, all the computations described below were run using 3 nodes of The TaiShan 200 cluster hosted by Huawei.Each node in the cluster contains 2x 64-core Huawei Kunpeng 920-6426 CPUs operating at 2.6 GHz (128 cores per node).
A first baseline computation was conducted using the G1 grid and adopting the dynamic Smagorinsky SGS model, in which the sub-scale filter length, ∆, is computed as the cubic root of the cell volume.The obtained results are compared with three different sets of experimental data from the literature, [11,12,15].It is very clear that the above described approach is able to properly predict the U-shaped profile of the mean stream-wise velocity, Fig. 2, in the wake region at x/D = 1.06.Furthermore, the other velocity profiles are also in very good agreement with the reference literature data.Moreover, the effectiveness and the reliability of the proposed solution technique is also confirmed by observing the mean transverse velocity (see Fig. 3), and the mean resolved streamwise Reynolds stresses in the wake region (see Fig. 4).It is worth noting that flow statistics were collected after the flow field was fully developed, by averaging a period of time corresponding to 500 non-dimensional units.A second campaign of analysis was conducted adopting G2 and G3 grids, in order to investigate the possibility to reduce the computational load required to proper reconstruct the flow field developed around the 3D cylinder.The adopted SGS model was also the dynamic Smagorinsky in this case, while for what concerns the sub-scale filter length, it was evaluated both with the cubic root of the cell volume and the least square technique.Looking at the mean velocity profiles in Fig. 2, it is evident that the adoption of ∆ cbrt results in a V-shaped profile in the wake region at x/D = 1.06, and the produced solution does not match Parnaudeau's data for the mean transverse velocity.On the other hand, the least square filter length does not exhibit a similar sensitivity to spanwise resolution.The LSQ filter enables to increase the grid anisotropy without compromising the solution accuracy, thus resulting in a reduction of up to 50% in the total number of cells, with respect to those required with the cubic root filter length.The mean streamwise transverse resolved Reynolds stresses depicted in Fig. 4    Nevertheless, the overall agreement with the reference can be considered satisfactory.

Conclusions
In this paper we have presented several sets of LES computations of the flow field developing over a cylinder operating at Re = 3900 and a free-stream Mach number equal to 0.2.The main goal of the paper is to assess the impact of different sub-grid characteristic lengths on the compressible LES simulation accuracy.In this context, it was noted that the adoption of the flow-dependent, least square length introduced by Trias et.[1] offers a clear advantage over the standard, geometry-based (cube root) length scale, in terms of accuracy.Indeed, the evaluation of the LSQ length increases the number of floating-point operations (FLOPS) with respect to the simpler cube root length.However, it enables to significantly reduce the overall number of grid cells that are required to correctly capture the complex physics of separated flow.Overall, our measurements demonstrated that the total time-to-solution is in-fact reduced when using the LSQ length scale, considering the combination of the two aforementioned effects (i.e increase in FLOPS and reduction in cell count).This behavior may be extremely appealing especially in the context of large-scale, industrial LES.

Figure 2 .
Figure 2. Mean streamwise velocity at different locations in the wake region.

Figure 3 .
Figure 3. Mean transverse velocity at different locations in the wake region.

Figure 4 .
Figure 4. Mean transverse velocity at different locations in the wake region.Streamwise direction

Table 1 .
Computational grids features.Grid N θ N r N z

Table 2 .
Drag and lift coefficients.G1 − ∆ cbrt G2 − ∆ lsq G3 − ∆ lsq Lysenko et al. lsq when the spanwise resolutiong is decreased.In fact, the coarse G3 grid also produces acceptable results, but it requires more than double of the computed dimensionless time to converge to a stable solution, compared with G1 and G2.Finally, Tab. 2 summarizes the mean drag coefficient, ⟨C D ⟩, and root mean square of the lift coefficient, C Lrms , deriving from the above-described computations.Both ⟨C D ⟩ and C Lrms are slightly overestimated if compared with Lysenko's results, in all the adopted approaches.