Low-order isotropic spatial resolution characteristics of the viscous terms for LES analysis on iso/anisotropic steady turbulent flows

This study aims to clarify the effect of reducing the isotropic spatial resolution on the turbulence field obtained in the Large Eddy Simulation (LES). Here, isotropic or anisotropic steady turbulence is analysed in this LES. A linear forcing method is used to generate these steady turbulent flows. Although the spatial resolution can be enhanced by increasing the number of spatial grid points or by using higher-order discrete forms, the isotropic resolution is only sometimes maintained. This viewpoint is expected to be particularly evident in the present study, especially in the LES analysis of wall turbulence. The present study focuses mainly on the spatial resolution of the viscous terms, where the spatial resolution can be varied while maintaining isotropy. Here, the discretisation accuracy order of the viscous terms was set from second to eighth order. In this study, the discretisation accuracy order of the convection terms is also varied. As shown in this study, the effect of the viscous term on the turbulence fields is larger than that of the discretisation accuracy order of the convection terms. There is a difference in the characteristics of the turbulence field obtained when the discretisation accuracy of the viscosity terms is of second and fourth or higher order.


Introduction
Numerical analysis of the turbulence field requires the governing equations to be discretised in time and space.Increasing the spatial resolution when analysing turbulence fields improves the accuracy of the resulting turbulence fields [1].The spatial grid size, determined by the number of grid points, which affects the spatial resolution, should be sufficiently small to resolve the turbulence structures in the field of turbulence to be analysed.On the other hand, as seen in the definition of the Courant number, a small spatial grid width reduces the upper limit for setting the time increment of the time integration.Therefore, increasing the spatial grid points significantly increases the computation required.The governing equations used in the analysis of turbulence are the continuity equation and the Navier-Stokes equation.These equations involve first-and second-order derivatives.Increasing the accuracy of the derivatives that approximate these derivatives can improve spatial resolution.
The central difference method is often used to spatially discretise the governing equations for incompressible turbulence fields.The governing equations of kinetic energy conservation derived from the spatially discretised governing equations should also be discretised and satisfied.Secondorder central difference methods are often used to discretise the convection and pressure terms.
Higher-order difference schemes have also been proposed for discretising these terms [2].The spatial resolution of the central difference methods depends on the accuracy order, which in wavenumber space is expressed in effective wavenumbers.As seen in the wavenumber distribution of the effective wavenumbers, decreasing spatial resolution can be viewed as underestimating the output values in the high wavenumber range.The amount of this underestimation in the high wavenumber range is larger in the second-order derivative than in the first-order.Therefore, the spatial discretisation may be less accurate in the viscous terms than in the convection terms in the governing equations.As seen in previous studies [3,4], increasing the order of spatial accuracy of the viscous terms rather than the convection terms can improve the accuracy of the obtained turbulence field.
The accuracy of the resulting flow field can be improved by increasing the spatial resolution used to analyse unsteady turbulence.This increase in spatial resolution can be conducted through the number of grid points and the order of discretisation accuracy.In contrast to the periodic cubic computational domain, the spatial resolution used for the analysis of wall turbulence, which is a type of anisotropic turbulent flow, is often set to be anisotropic [3,5].Specifically, the number of grid points used and the order of discretisation accuracy are often different in the direction normal to the wall surface and the directions parallel to the wall surface.Therefore, to investigate the response of anisotropic turbulence to spatial resolution, aspects of the relationship between the isotropy of the computational grid and the anisotropy of the flow field need to be addressed.This study considers that the spatial resolution dependence of anisotropic turbulence under isotropic spatial resolution should be investigated.
The aim of this study is to reveal the effect of reduced isotropic spatial resolution in numerical simulations on isotropic and anisotropic steady turbulence fields using a large eddy simulation.The reduction of spatial resolution maintained isotropically is set up by reducing the discretisation accuracy of the viscous terms using an isotropic computational grid.Here, the present study uses second-to eighth-order forms to discretise the viscosity terms.The present study also focuses on the accuracy order of the discretisation of the convection terms.The present study approaches this problem using the Reynolds number dependence of the fundamental statistics.Therefore, this study uses the Vreman model [6] as the sub-grid scale (SGS) model for large-eddy simulation (LES).The Vreman model is a widely used SGS model and, in contrast to the Smagorinsky model, can be used for simulating disturbances at low Reynolds numbers.The isotropic and anisotropic steady-state turbulence fields in this study are established using a linear forcing method.

Governing equations and forcing terms
Effects of discretisation accuracy from low to high order on the steady turbulence fields are examined in this study by applying isotropic spatial resolution.The spatial resolution dependence on the Reynolds number is investigated using LES by analysing the flow fields in the low to high Reynolds number range.
The governing equations for the flow fields of interest are the non-dimensionalised continuity equation and the Navier-Stokes equations, details of these eqations are shown in the previous study [1].An external force term F due to the linear forcing method is included in the governing equations.Steady turbulence is generated by imposing this external force term.The discretised governing equations include a first-order difference formula for convection and a second-order difference for viscous terms.This study focuses on increasing the spatial resolution of these two terms to high-order accuracy.In this study, a periodic cube computational domain with a side of 2π is used.The reason for using this computational domain is that previous studies have often used it.The number of grid points is 32 in each direction of the computational domain.Since errors in spatial discretization decrease with decreasing grid width, effects obtained in this study are expected to possibly decrease with decreasing grid width.
The present study adopts the linear forcing method [7][8][9][10], in which isotropic and anisotropic external force terms are given to generate steady turbulence.The linear forcing method generates steady turbulence by giving the external force term components in proportion to velocity components.The external force F is given by F = Cu.Here the velocity components of u, (u x F , u y F , u z F ), are set as those of the analytical solution, and C is a constant.The following Equation (1) shows the components given to produce the isotropic steady turbulence, and Equation (2) shows the components to produce the anisotropic steady turbulence.u x F =  cos(x) sin(y), u y F =  cos(x) sin(y), and u z F = 0. (1) and

Simulation methods
For the convection terms, Morinishi's discretisation method [2] is used.This study uses the skewsymmetric type, where the squared velocity components are essentially conserved.The viscous terms are the diffusion terms in the equations of motion.Most differential approximations of the secondorder differential coefficients of the viscous terms control the velocity fluctuations and dissipate the kinetic energy.The grid size can be reduced for more accurate analysis, but higher-order accuracy of the difference method is required.The Vreman model is used to model the LES.The Vreman model is capable of analysing not only well-developed turbulent flow fields but also flow fields in transition regions.A value of the model constant C v in the Vreman model is obtained from that of the Smagorinsky model constant C s by C v  2.5 (C s 2 ).Central difference methods are used for the discretisation of the convection and viscosity terms.The convection term, which is a first-order difference term, is analysed using second-and fourth-order difference methods [2][3][4][11][12].The viscous terms are discretised using the second-order difference.Second-, fourth-, sixth-and eighth-order accuracy methods are used to discretise the viscous terms.In this study, the same accuracy order was used in the x, y and z directions and isotropic spatial resolution was used for the validation.For the fourth-order accuracy of the convection terms, viscous terms with second-, fourth-, sixth-and eighth-order accuracy are analysed.A total of eight sets of conditions are used to examine the effects of reducing resolution.

Numerical conditions
The computational conditions used in this study were the staggered grid and the computational domain (2) 3 .The time increment t was set to t = 0.005, and time evolved up to t = 3000.The constant in the linear forcing method C was set to C = 1.Reynolds numbers were set to Re = 20, 30, 50, 100, 200, 300, 500, 700, and 1000.The number of spatial grid points is set to N 3 = 32 3 .The six-stage fourthorder accurate Runge-Kutta method [13] was used for the time integration method, and the fractional step method was applied for coupling the continuity equation and the pressure field.In this fractional step, the Navier-Stokes equation without the pressure term is first integrated over a small time period.The pressure is then calculated by solving the Poisson equation for pressure with the term based on the divergence of this obtained velocity field as the right-hand side.The velocity field obtained by taking the gradient of this pressure is then imposed on the previously obtained velocity field to obtain the velocity field at the next instant in time.This procedure is carried out in the present study at each stage of the Runge-Kutta method.The residuals of the equation of continuity are of machine-zero order in double-precision floating-point computing.The overall computation time for one case is approximately tens of hours.In this study, fast computation is achieved by solving the Poisson equation using the Fast Fourier Transform.
Since the results of this study are related to spatial resolution in numerical analysis, this study sees the possibility that experimental studies such as wind tunnel experiments may not provide reference results.This study considers that further calculations should be performed as future research to validate the results of this study.The numerical analysis code used in this analysis has already been validated in previous studies [3][4][11][12], and relevant previous studies are cited in the revised manuscript.This study only limits itself to presenting numerical results.The reason for the limitation is that large-eddy simulation is used in this numerical analysis in view of engineering applications.Different from a direct numerical simulation, this study considers that in the large-eddy simulation, the results obtained are influenced by the sub-grid scale model.Future direct numerical calculations should be applied to fully elucidate the results of this study.

Turbulence fields and kinetic energy
Figure 1 shows the visualisation results for Re = 50 and Re = 500 using isotropic steady turbulence.All results are obtained by applying the second-order accuracy form to the convection terms, while the viscous terms are discretised using second-or fourth-order accuracy.Here, the green isosurface is characterised by a negative value of pressure P [14][15][16][17], indicating a large-scale turbulence structure.As seen in the results for Re = 50, the results obtained by setting a second-order accuracy form for the viscous terms can differ slightly from the results obtained by applying a fourth-order accuracy form for those terms.These results suggest that turbulence characteristics may be quantitatively different between the second-and fourth-order accuracy forms of discretising the viscous terms for the Re = 50 condition.
The results for Re = 500 also show the difference in the order of accuracy of the discretisation of the viscous terms.The present study finds that the difference in the spatial accuracy order for discretising the viscous terms is visualised in the large-scale eddy structure of the isotropic field.Figure 2(a) shows results for the dependence of the turbulent kinetic energy on the Reynolds number for isotropic steady turbulence.Figure 2(a) shows that in the low Reynolds number range seen as Re = 20 to 50, the red plots where the viscous terms discretisation is second order agree with the blue and black plots obtained by applying higher-order accurate discrete expressions for the viscous terms.On the other hand, around Re = 100 ~ 700, values shown in the plots where the viscous terms are second-order accurate are smaller than those shown in the plots shown in blue and black.In this Reynolds number condition, the blue plots where the viscous terms are at least fourth-order accuracy are in good agreement with the plots in black.The plots agree with the second and fourth orders of accuracy, discretising the convection terms.This result indicates that the effects of the accuracy orders of different convection terms are smaller than those of the viscous terms.
Figure 2(b) shows results for the dependence of turbulent kinetic energy on the Reynolds number using anisotropic steady turbulence.As shown in the figure, the plots agree well when higher-order accuracy discrete forms are applied to the viscous terms.The values of the red plots obtained using the second-order accuracy forms for the viscous terms are smaller than those obtained by applying the high-order accuracy discrete formula to the viscous terms.For different orders of the accuracy of the convection terms, the values of the plots for the second-and fourth-order accuracy forms are in close agreement in the low Reynolds number range, while a disagreement is observed in the high Reynolds number range.Figure 3 shows the visualisation results for Re = 500 with anisotropic steady turbulence.For both results, the discretisation accuracy order of the convection terms is second order, and the orders of the viscous terms are second and fourth order.The white isosurface is characterised by a negative pressure value, while the light blue isosurface is characterised by zero pressure.As shown in the figure, twodimensional eddy structures can be seen as a characteristic feature of anisotropic turbulence.Focusing on the diameter of the central eddy tube, the present study shows that the diameter obtained with the viscous terms set to second-order accuracy is larger than that obtained as fourth-order accuracy.On the other hand, the overall shape of the eddy structure agrees between the viscous terms discretisation accuracy orders.As shown in the figure, in the high Reynolds number range, the distribution tendency of the turbulent kinetic energy is similar between the different orders of viscous terms discretisation accuracy, and the values differ.This study considers that the visualisation results obtained can be related to the result found in the turbulent kinetic energy distribution.

Turbulent Reynolds number
Figure 4(a) shows the results for the dependence of the turbulent Reynolds number on the computational Reynolds number using isotropic steady turbulence.As shown in the figure, the values of the results are in close agreement between the conditions where the viscous terms are set to fourthorder accuracy or higher.However, the plots obtained with the viscous terms set to second-order accuracy do not agree in value with the other plots in the low Reynolds number range.This deviation decreases with increasing the Reynolds number.For the difference in the convection terms, the plots obtained using the second-order accuracy form for the convection terms were nearly agreed with those obtained using the fourth-order accuracy form.The effect of the reduction in spatial resolution was more evident when the viscous terms were set to second-order accuracy.
Figure 4(b) shows the results of the dependence of the turbulent Reynolds number on the computational Reynolds number in anisotropic steady-state turbulence.The blue and black plots obtained with the viscous terms discretised to fourth-order accuracy or higher are in close agreement with each other.On the other hand, the red plots obtained where the viscous terms were discretised in second-order accuracy form were not in agreement with the other two types of plots.Between the second and fourth-order discretisation schemes of the convection terms, the results illustrated in the plots agreed with each other.Therefore, the decreasing spatial resolution effects are more evident for the viscous terms, even for the different accuracy orders of the convection terms.Similar to the results for anisotropic fields, the plots for the viscous terms discretised in the second-order accuracy form approach the results of the other two types in the high Reynolds number range.In the turbulent Reynolds number results, for both isotropic and anisotropic fields, the effects of the different spatial accuracy orders of the convection terms are smaller than those of the viscous terms.

Small-scale turbulence characteristics
Figure 5 shows results for the Reynolds number dependence of small-scale isotropy with isotropic and anisotropic steady turbulence.Here the parameter shown in the figure is defined as follows.
This measure should be unity in the isotropic and anisotropic turbulence fields analysed in this study.
Therefore, this condition can be used to validate the results of the present analysis.As shown in the figure, the values of the parameter obtained in the isotropic and anisotropic turbulence fields are unity for all conditions.The present study considers that this agreement validates the present analysis.Figure 6 shows the values of the anisotropy parameters for small-scale fluctuations in the anisotropic field.Here, two parameters describing the anisotropy of small-scale fluctuations are defined as follows.
The measure shown in Equation ( 4) is characterised as the ratio of the intensity of the fluctuations of the derivatives between the x-direction and the z-direction, for which the field is statistically homogeneous.The measure shown in Equation ( 5) is the ratio of the intensity of fluctuations of the derivatives of the velocity components in the x-direction and in the z-direction, where the field is statistically homogeneous with respect to the y-direction.As shown in the figure, the values of the plots are in agreement between the accuracy orders of the viscous terms under high Reynolds number conditions.
As the Reynolds number is decreased, the different discretisation accuracy orders of the viscous terms can affect the values of the parameters.Specifically, when the second-order accuracy discretisation formula is applied to the viscous terms, the values of both indicators are small at low Reynolds numbers.This underestimation is commonly found between the two indicators characterising the anisotropy of small-scale fluctuations.
Figure 7 shows visualisation results for the intensity of derivative fluctuations constituting the small-scale anisotropic parameter defined by Equation ( 5) in anisotropic steady turbulence at Re = 50.Here, for these visualisation results, the convection terms discretisation accuracy is second order, and the accuracy order of the viscous terms is second or fourth order.As shown in the figure, the visualised structure is fragmented and characteristic structures are not seen when the viscous terms discretisation accuracy order is second-order.On the other hand, when the discretisation accuracy order of the viscous terms is fourth order, distinctive structures can be seen in the isosurfaces with both fluctuation intensities.These results may be related to the difference in small-scale anisotropy parameter values between accuracy orders of the viscous terms.The present study considers that the lack of distinctive structures reproduced in the derivative fluctuation intensities may be related to the underestimation of parameter values found when the discretisation accuracy order of the viscous terms is second order.A spatial grid point size of 32 3 has been set as the number of spatial grid points in this study.This study considers that the effects of the isotropic reduction in spatial resolution, which this study focuses on, are more significant in the small-scale domain.Improving the spatial resolution by increasing the number of spatial grid points improves the resolution in the small-scale region related to spatial discretisation and thus reduces the difference in the magnitude of underestimation of the difference values in the small-scale region due to differences in spatial accuracy order.Therefore, this study considers that increasing the number of spatial grid points can reduce the differences between the set spatial resolutions that represent the impact approached by this study.Increasing the spatial grid points may reduce the magnitude of SGS eddy viscosity.
As this study uses the Vreman model as the SGS model, the results of this study are expected to apply not only to the analysis of fully developed turbulence but also to the analysis of phenomena related to laminar-turbulent transition.Therefore, the results of the present study may be applied to the analysis of flow around an aerofoil.The laminar-turbulent transition can also be observed in wake turbulence.The results of the present analysis can also be used to analyse non-equilibrium decaying turbulence involving wake turbulence in LES, which has been the focus of previous studies [21][22][23][24][25][26][27][28][29][30].

Conclusions
The aim of this study is to clarify the effects of isotropic reduction of spatial resolution in numerical analysis on isotropic and anisotropic steady turbulence fields using large-eddy simulation.In this study, the spatial resolution of the viscous terms was set to be the same in the x-, y-and z-directions, and second-order to eighth-order differences were used in the analysis, as the previous studies often assume that the accuracy order of the viscous terms is the same for the convection and the viscous terms.Second-and fourth-order first-order differences were also used to discretise the convection terms.Validation was carried out under eight different conditions to discretise the convection and viscous terms.The computational Reynolds number was set between 20 and 1000, and sub-grid scale fluctuations were modelled using the Vreman model.Steady turbulence was generated using isotropic and anisotropic external force terms.The effects of the isotropic reduction in the spatial resolution were examined using turbulence statistics and visualisation results.
A difference existed in the turbulent Reynolds number results between the second and higher orders of the viscous terms discretisation accuracy.The turbulence energy also observed the effects of reducing the spatial resolution by setting the viscous terms to second-order accuracy.Visualizing large-scale eddy structures in the isotropic fields examined the effects of the difference in viscous terms discretisation accuracy.Large-scale structures in anisotropic fields were two-dimensional, and the diameter of the eddies depended on the discretisation accuracy order of the viscous terms.Smallscale anisotropy significantly depended on the discretisation accuracy order of the viscous terms.This result was discussed using the intensity of the spatial derivatives of the velocity fluctuations.

Figure 1 .
Figure 1.Large-scale turbulence structure in the isotropic fields compared between viscous term discretization orders.Here, isosurface characterised by the negative pressure value is used for visualising the turbulence structure.

Figure 2 .
Figure 2. Dependence of mean turbulent kinetic energy on computational Reynolds number.

Figure 3 .
Figure 3. Large-scale turbulence structure in the anisotropic fields compared between viscous term discretization orders.Here, isosurface characterised by zero and the negative pressure values is used for visualising the turbulence structure.

Figure 4 .
Figure 4. Dependence of turbulent Reynolds number on computational Reynolds number.

Figure 5 .
Figure 5. Validation of turbulence fields using anisotropic parameter values expected to be isotropic.

Figure 6 .
Figure 6.Dependence of the anisotropic parameters on the computational Reynolds number.

Figure 7 .
Figure 7. Visualization of small-scale fields using velocity gradient components between the discretization orders of the viscosity terms.Here, isosurface characterised by the gradient intensity values is used for the visualization.4.Discussions and Conclusions4.1.DiscussionsIn this study, the Vreman model is used as the SGS model, which can be used to analyse turbulence fields in high Reynolds number conditions and low Reynolds number conditions.The fact that the Reynolds number dependence of the turbulence fields is used in this study to solve this problem corresponds to the application of the Vreman model.On the other hand, there are other SGS models, such as the conventional Smagorinsky model and the implicit LES model.In particular, the Smagorinsky model is often used in commodity fluid solvers such as OpenFOAM[18][19][20].As shown in this study, the results obtained using these two SGS models are in qualitative agreement with the results of the present study.This agreement indicates that the effect of the different SGS models on the results obtained in this study may be small.The Smagorinsky model used by this study as an SGS model includes one model constant.The value of this model constant needs to be set in order to use this SGS model.The value of this model constant directly determines the magnitude of SGS eddy viscosity and SGS dissipation.In the present study, a widely used value of 0.1 is set for this model constant.The study considers that the uncertainty in setting the value of this model constant is not large enough to qualitatively influence the conclusions of the study.Results for the flow fields obtained by optimising the value of the model constant could be quantitatively more accurate.Optimising the value of the model constant may contribute to a quantitatively more rigorous investigation of the effects of isotropic spatial resolution