LES mesh resolution requirements for particle-driven gravity currents

In this study we investigate the influence of grid resolution on a near-wall resolved LES model of a lock-exchange particle-driven gravity current. The simulations are performed using the finite volume Boussinesq code SnS with a Smagorinsky turbulence model for a buoyant Reynolds number of 60,000 on 4 grid sizes. According to previous studies, two-point correlations are most appropriate to estimate LES resolution. With the largest scales of the flow being resolved by more than 20 cells, well-resolved LES is obtained for grid resolutions of 1925×62×125 and finer. In addition, in order to apply the turbulence model correctly, we show that the velocity power spectrum densities provide useful information for the maximum cell size. The ratio of the subgrid scale viscosity to the molecular viscosity and the subgrid scale shear-stress to the resolved Reynolds stress show good convergence with grid refinement. The ratios 0.3 SGS    above the current and ( ' ') 0.05 SGS ave u v   inside the mixing layer, are chosen as threshold values, based on our evaluation study.


Introduction
Turbidity currents are an extended form of gravity currents where the individual inertial effects of each particle, which tend to settle down due to their larger size, become nonnegligible and impact the evolution, sustainability and sediment deposition of the particle wave, particularly during the later stage of its propagation [1,2].Despite an early identification of the existence of particle charged underflows inside the Rhone River by Forel [3], the interest in sediment-driven underwater flows has only started growing during the past forty years with the recognition of their impact on submarine infrastructure.Experimental studies have aimed to describe their large scales features [1,[4][5][6][7], as well as their behaviour in specific conditions such as when turbidity currents impinge on obstacles [8].To compensate for their intrinsic opacity, and the intrusive nature of the experimental measurements, computational fluid dynamics (CFD) has become increasingly used as a method to study the inner turbulent dynamics of turbidity currents.
Since the first Direct Numerical Simulation (DNS) by Härtel, et al. [9], the characteristics of gravity currents have been thoroughly investigated numerically [10][11][12].Necker, et al. [13] proposed an extended version of a model accounting for the particles settling within turbidity currents, which was implemented in the CFD code TURBINS [14] and applied to the study of turbidity currents interacting with complex sea-floors [15].The computational requirements limit the use of high resolution simulations to fairly low Reynolds numbers and Large Eddy Simulation (LES) has been shown to be an alternative method to reduce the computational times for high Reynolds number flows [16].
Grid resolution is a key parameter for a proper representation of the fluid features by an LES model, however the dependence of the LES subgrid scale (SGS) models on the cell size renders it difficult to specify general guidelines on how to assess the quality of the simulations.Several estimators have been developed to quantify the numerical errors such as the subgrid activity parameter s [17] or the relative Kolmogorov scale, SGS viscosity and resolved turbulence kinetic energy indices of resolution quality LES_IQ, LES_IQ and LES_IQ k [18].Their calculation requires either experimental or DNS data often unavailable, or the application of the Richardson extrapolation to estimate the potential DNS total kinetic energy [18,19].Therefore, the reliability of the latter method is questionable [20].It may lead to overestimation of the indices caused by the difficulty to represent the strong interdependency between the numerical error, due to the spatial discretization, and the modelling error, due to the turbulence modelling, which may sum up or negate the total error.
Davidson [21] proposed a case related approach based on the comparison of several turbulent statistics in areas where the turbulence is critical in the flow.It was shown that twopoint correlations and the ratio of the SGS shear-stress to the resolved Reynolds stress are the most suitable characteristics to assess the grid resolution, whereas the turbulence spectra do not seem to provide any conclusive results.As for the ratio of the SGS viscosity to the molecular one, its reliability is arguable.Despite showing good convergence for the LES of the dispersion around an isolated cubic building [22], it leads to incorrect evaluation of the resolved turbulence within the domain for a channel recirculating flow [23] and shows better resolved turbulence in areas where the other parameters indicate the opposite.
The purpose of this study is to qualify the applicability of these parameters to the assessment of grid resolution in the case of three dimensional lock-exchange turbidity currents at high Reynolds number.

Problem definition and numerical model
The 3D model is inspired by the standard lock-exchange setup of the experiment of Wilson, et al. [8], as shown in Fig. 1.The header box is filled with sediment-laden fluid and the flume channel is filled with the same fluid without sediment.At 0 t  , the sediment-laden fluid is released and the current starts flowing into the channel.The problem is described by the continuity and the momentum equations together with a transport equation for the dimensionless particle concentration.An in depth description of the coupling between the equations is given by Necker, et al. [13] and the approach is briefly summarised in this section.Contrarily to the standard gravity current formulation which rely on a Boussinesq approximation to qualify the buoyant forces caused by the fluids inner density variations, the driving forces are directly derived from the force applied by the particles on the fluid.In the case of spherical particles of diameter d p smaller than the Kolmogorov scale, the effort caused by the concentration distribution on the fluid fp is modelled by the product of the Stokes drag force on an individual particle Fp and the particle volume fraction m.The response of the particles to the flow motion is assumed to be faster than the time scale of the smallest flow structures, so that the particle velocity associated to the concentration transport can be corrected as + p S u u u  .u is the fluids velocity and S u the vertical settling velocity of a particle inside an immobile fluid derived from balancing the Stokes drag force with the Archimedes force.
( -3 ) The settling velocity is set to 0.0035

S u 
. The equations are presented in their dimensionless form.The characteristic quantities for scaling the variables are chosen to be half of the channel height H/2, the maximum particle volume fraction m r , the buoyant velocity and the time scale . We also define the buoyant Reynolds number the ratio between the two dimensionless characteristic numbers of the flow.The LES model is based on the filtered form of the continuity, momentum and particle concentration equations.Using as the filter operator, the equations are The simulations are performed using the structured non-staggered Cartesian finite volume SnS code developed by Norris [24].The filtered equations are solved for using a fractionalstep method, where the advection terms are discretized in time using an Adams-Bashforth scheme and the diffusion terms using the Crank-Nicolson scheme.A second order central differencing scheme is applied for the spatial discretization of the advection and diffusion terms.The momentum and pressure Poisson equations are solved using a bi-conjugate gradient stabilised (BICGSTAB) method with a strongly implicit procedure preconditioner.A Jacobi method is used for the concentration transport equation.The time step varies to maintain the value of the Courant-Lewy-Friedrich number ( i C L t F u x    ) within the range 0.15-0.25,and convergence is assured at each time step by ensuring a velocity divergence of less than 10 -6 .The code was originally designed for the efficient modelling of buoyancydriven convective flows, and has been applied to the modelling of natural [25] and mixed convection [26], and to the LES modelling of wind turbines [27].
Four grid resolutions were compared in this study, as listed in Table 1.A uniform grid size was chosen except in the near-wall region.To ensure a proper resolution at the bottom wall, the first vertical grid size was set to 0.00133 y  which ensured it was located in the near-wall viscous sublayer ( 5 y   ).The vertical mesh spacing was progressively increased over the layer 0 0.667 y   and the number of cells N y,w was chosen so the ratio between two consecutive vertical mesh size remain in the range 0.85 1.15 y    .

Results and discussion
The propagation of a gravity or turbidity current is typically described by the evolution of its front position and velocity, shown in Fig. 2. Limited sensitivity to the grid resolution is observed.The simulations predict the initial slumping phase with a similar quasi-constant velocity before decreasing following a t -0.56 power law.The transition occurs after 16 corresponding to eight lock height, which lies within the range of 5H to 10H, as generally observed in experiments [8,28].The deceleration rate lies halfway between the characteristic t -1/3 and t -4/5 power law observed in gravity currents.The former describes the self-similar phase dominated by buoyant and inertial effects and occurs just after the slumping phase [29].The latter represents the buoyant-viscous phase, dominated by viscous forces, distinctive of the current's later stage [28].This highlights the impact of the settling velocity on the efforts balance at early stage of the currents propagation [15].Fig. 3 shows the impact of the grid resolution on the internal structure of the current.Upon the gate release, the typical Kelvin-Helmholtz instability, characterised by large billows at the upper boundary of the current and large coherent vorticity structures, is observed as well as the lobe and cleft instability at the front of the current are observed.The structure of the billows is blurred by the presence of the smallest turbulent structures resolved at finer grid resolutions.Those structures have been shown to be the result of shear effects at the front propagating into the body induced by the lobes [16].However, refining the grid also reduces the size of the lobes while increasing their number, so that the turbulence production at the front is altered.The induced modifications inside the current can either enlarge or negate the simulations errors and shows a clear need for qualifying the lobes effect on the turbulence creation to prevent from misinterpreting the influence of grid refinement on the simulations.According to Davidson [21], [23] and Bazdidi-Tehrani, et al. [22], two-point correlations are most suited to assess the grid quality of a LES.From them, we can find the number of grid points needed to resolve the largest structures of the flow.It is given that a minimum of 8 cells is required to achieve a coarse LES and 16 cells for a well-resolved LES [21].The streamwise two-point correlations indicate that the largest structures are resolved with more than 20 cells for each case, therefore only the spanwise two-point correlations R uu (z) and R ww (z) averaged over the whole domain are displayed in Fig. 4. The correlations R uu and R ww decrease to 0.1 within only 11 and 8 cells respectively for grid G1, and 17 and 13 for G2 highlighting the coarseness of this mesh.A well-resolved LES is attained with G3 and finer, with a minimum of 30 and 22 cells used to resolve the largest structures.
LES requires the positioning of the filter's cut-off wavenumber of the turbulence model on the decaying k -5/3 inertial subrange of the turbulence spectra.While Davidson [21], [23] concludes that the energy spectra are of little use in the assessment of grid resolution for LES, they still provide important information to justify the latter condition.The spanwise power spectrum densities PSD z {U} and PSD z {W} are plotted on Fig. 5.They have a standard structure regardless of the resolution, even though the PSD z {W} exhibit a decaying range smoother than the typical t -5/3 law and surprising "pile-ups" at high wavenumbers.This behaviour is also observed for a channel flow by Davidson [21].He attributes it to inaccuracy of the turbulence model which fails to predict the SGS dissipation at very high wavenumber and transfers the dissipative functions at lower wavenumbers.The LES filter's width being defined by the cell size, the cut-off wavelength corresponds to the largest wavenumber representing the PSD z s.Both PSD z s do not reach the k -5/3 decay inertial range for mesh G1, indicating an inappropriate filter width and an overly coarse mesh size for the turbulence model.As for the two-point correlations, grid G2 barely reach the k -5/3 range and can be considered a coarse LES, while G3 and G4 show the correct decay range for both components.Particle-driven gravity currents are not quasi-steady.Time averaging, which is typically used for the statistical analysis of LES models, is thus not applicable.Furthermore, the prediction of the position of the current and its internal dynamics varies significantly between simulations, and so a comparison of the instantaneous fields is thus unreasonable.Consequently, instead of time-averaging the flow is spatially averaged in the spanwise direction and in a sufficiently long streamwise layer to have a sufficiently large number of points to give a meaningful average.The layer is chosen in the most energetic region of the current with developed turbulence and minimal vertical variation, and the extractions are taken at 29.1 t  corresponding to the time of current's peak of turbulence.Based on the vorticity fields (Fig. 3) and the turbulent kinetic energy fields (not shown), the layer is 4 2 The direct observation of the averaged streamwise velocity ave U and resolved Reynolds stress ( ' ') ave u v , see Fig. 6. do not show convergence.The streamwise velocity predicts similar results with differences at the area of rapid variation of its gradient.Likewise, ( ' ') ave u v evolves without converging with grid refinement.This can be explained by the flow's non-homogeneities and the limited number of averaging points.Small changes in the layer length or position can lead to significant differences and it is worth indicating that similar conclusions are obtained at different Reynolds numbers (not shown).The ratios of the SGS viscosity to the molecular one, and the SGS shear-stresses to the resolved Reynolds stresses (Fig. 7), are suited for grid resolution assessment [21][22][23].The high mixing region stress ratios (Fig. 6) correspond to the peak of ( ' ') ave u v .Both ratios show convergence to lower values when the grid is refined.High values indicate poor grid resolution.Care needs to be taken in the analysis of the peaks of ( ' ')

Conclusions
Several LES grid resolution assessment methods have been compared for the study of lockexchange sediment-driven gravity currents at a buoyant Reynolds number of 60,000.Unlike streamwise values, which do not give any particular information to assess the quality of the simulations, it is concluded that the most useful quantities are the spanwise two-point correlations.Those quantities inform on the number of cells to be used to resolve the largest scales of the current.It is the responsibility of the user to choose a suitable number of cells, therefore 16 cells is commonly considered a minimum for a well-resolved LES.It is found that the grid G2=1140×37×74 is too coarse for LES, whereas G3=1925×62×125 enables a well-resolved LES solution.The velocity power spectrum densities lead to similar conclusions.Only resolutions finer than G2 respect a proper positioning of the filter's cutoff wavenumber in the inertial subrange of the turbulence spectra.
For statistical assessment, difficult for a non-quasi-steady current, values are averaged on the spanwise direction for the strongly turbulent and energetic streamwise layer 4 2 . The ratios of viscosities SGS   and stresses ( ' ') SGS ave u v  decrease with grid refinement inside the current body.However, threshold values for a well-resolved LES are case related and no guidelines are generally given to quantify them.We show that recommendations can be made with the help of the two-point correlations and power spectrum densities.Thus, we conclude that a complete resolution of the flow in the lowly turbulent region above the current, should be reflected by

Fig. 3 .
Fig. 3. Concentration isosurfaces 0.1 c  and vorticity fields on the symmetry plane 0 y  after release of the lock ( 9.3 t  , left) and for fully developed turbulence ( 29.1 t  , right) for grid sets G2 (top) and G4 (bottom).

Fig. 4 .
Fig. 4. Averaged two-point correlation of the streamwise (left) and spanwise (right) components of velocity along the spanwise direction at 29.1 t  .

Fig. 5 .
Fig. 5. Averaged power spectrum densities of the streamwise (left) and spanwise (right) velocity components along the spanwise at 29.1 t  .

Fig. 6 .
Fig. 6.Vertical profiles of the streamwise resolved velocity component and the resolved Reynolds stress averaged over the spanwise direction at 29.1 t  .
and G4.The high values of G2 and G3 are not due to an increase in the SGS shearstress compared to the Reynolds stress, but to an artefact resulting from the shrinking of ( ' ') ave u v at the sides of the mixing layer.The minimum observed for G4 is the result of the increase of SGS  from near-null negative values inside the current to a positive peak inside the mixing layer.The ratios SGS   increase as expected in the mixing layer and decrease inside the ambient fluid.Although threshold values for the peak are difficult to choose as no reference value to remain in an acceptable range is available in the literature, one can assume that the flow should be entirely resolved in the low turbulence region above the current.This is indicated by 1 a limiting value, due to G3 being resolved enough.Similarly, one can consider good convergence when the SGS shear-stress represents roughly less than 5% of the Reynolds stress inside the mixing layer.

Fig. 7 .
Fig. 7. Ratio of SGS viscosities and molecular viscosities (left) and SGS shear-stresses and resolved Reynolds stresses (right), averaged over the spanwise direction at 29.1 t  .