Assessing OpenFOAM-based Large-eddy Simulation Using Decay Characteristics of An Isotropic Taylor-green Vortex

The aim of this study is to investigate OpenFOAM using isotropic TGV flows. The main focus is on the establishment of energy conservation laws, in particular using inviscid fields to assess energy constancy; Large-Eddy Simulation (LES) is applied and the results are compared with those obtained using spectral methods; unlike previous studies of OpenFOAM, this study uses an unsteady turbulent TGV flow. The computational methods used are a central difference method, a PISO algorithm and a Smagorinsky model. The study revealed errors in the conservation of kinetic energy in the OpenFOAM analysis. In particular, the energy decay was significant at high spatial resolution, confirming the different turbulence energy decay characteristics of DNS and OpenFOAM. However, when the model constant was 0.6, the decay characteristics were shown to be in good agreement with DNS.


Introduction
The ambient flow associated with industrial and transport equipment is incompressible turbulent flow [1].Understanding the phenomena of incompressible turbulence is an engineering necessity and numerical analysis is the method of choice.Large-eddy simulation (LES) is widely used because such a number of grid fractions of Re 9/4 or more in three-dimensional space is rarely possible in Direct numerical simulation (DNS) [2].Fluid dynamics solvers are also widely used in the numerical analysis of turbulent flows.Fluid analysis solvers are commonly used to solve unsteady flow fields in complex geometries.However, these solvers do not use accurate discretisation methods or spectral methods as shown in previous studies.Therefore, there is a need to investigate the accuracy of the numerical analyses performed by the fluid analysis solvers.
One of the generic fluid analysis solvers is OpenFOAM [2,3].OpenFOAM is widely used in both basic and applied research.The accuracy of unsteady turbulence analysis using OpenFOAM has been validated in previous studies.In these validations, the accuracy of OpenFOAM was tested using a channel flow as a high-resolution quasi-DNS.High resolution fluid analysis has also been performed previously.The most accurate method of analysis is a spectral method (e.g., [4]), but only a limited number of flows can be subjected to spectral methods.A finite volume method is often used for high fidelity analysis [5][6][7][8][9].This finite volume method allows unsteady turbulent flows to be analysed by accurately establishing the physical conservation laws using appropriate discretisation schemes [7].To validate the high-fidelity fluid analysis, the finite volume method is used in the analysis of unsteady turbulent flows.
2 Previous studies that have validated OpenFOAM have used channel turbulence (e.g., [2]) and not Taylor-Green vortex (TGV) flows [10][11][12] for validation; TGV flows are unsteady turbulence unlike channel turbulence.Therefore, the use of TGV flows allows OpenFOAM to be validated in terms of unsteady fields.The TGV flows in previous studies are globally anisotropic.Therefore, unlike isotropic turbulence fields, there is an energy exchange between components due to pressure in this flow field.The presence of global anisotropy can obscure the results of validation using TGV flows.This study believes that this OpenFOAM validation should have a physical conservation law perspective.In particular, it should be investigated whether the energy conservation laws in the inviscid state hold in the OpenFOAM analysis.
In light of these perspectives, the aim of this study is to investigate OpenFOAM using an isotropic TGV flow.The perspective of the energy conservation law, from which the governing equations are derived, is also considered in the validation.The turbulence characteristics of the isotropic TGV flow analysed by LES are compared with those obtained by the spectral method.This flow field has no energy dissipation due to viscosity, which allows the initial field energy to remain constant with respect to time.

Numerical Methods
In this study, OpenFOAM-based LES is investigated using the decay properties of an isotropised Taylor-Green vortex [10,11].In particular, the results of OpenFOAM-based LES are compared with those of DNS.This validation also takes into account the model constants perspective.The governing equations for this analysis are the continuity equation for incompressible fluids and the Navier-Stokes equations.These two equations are non-dimensionalised.The analytical object of this study is an unsteady incompressible flow.The computational domain is a 2π-sided cubic domain with periodic boundary conditions on each surface.Such a domain setup is often used in previous studies.
An isotropic Taylor-Green vortex was used for the initial field to validate OpenFOAM.The Taylor-Green vortex is an unsteady flow with symmetry in the initial conditions.Analytical solutions exist and are often used as benchmarks.It is also isotropic, as it is expected to simplify the domain of interest and not have to consider the contribution of static pressure fluctuations.Isotropisation is also often used to validate the analysis.The equation for velocity components, u, v, and w, of the isotropic Taylor-Green vortex is expressed by the following equation, where (u,v,w) are for coordiates (x,y,z): u(x,y,z) = (2/3 1/2 ) (-cos(x) sin(y) + cos(x) sin(z)), v(x,y,z) = (2/3 1/2 ) (-cos(y) sin(z) + cos(y) sin(x)), w(x,y,z) = (2/3 1/2 ) (-cos(z) sin(x) + cos(z) sin(y)). (1) The equation consists of a combination of anisotropic terms with coefficients set such that the magnitude of the initial kinetic energy is unity.The flow field also satisfies the continuity equation.
Previous studies have shown that the conservation property of the fluctuation intensities of discrete expressions used in numerical analysis is explicit, and that the conservation property of the squared quantities is reduced when the time integration method is implicit [13,14].The explicit methods implemented in OpenFOAM are not accurate sufficiently to validate the analysis [5][6][7][8][9].Therefore, the second-order accurate Crank-Nicolson method, which is an implicit method, was used.
In this analysis, the second-order accurate central difference method is used as the spatial difference method [7].Central difference is an appropriate spatial discretisation method for this study as it does not involve numerical viscosity.A collocated grid is used here for this discrete equation.The PISO algorithm is used to couple the governing equations [2].This is a pressure-velocity coupling method applied to transient problems.The analysis also uses large-eddy simulation, an analytical method that models small-scale vortices and directly calculates large-scale flows.The Smagorinsky model, one of the most commonly used models, was adopted as the SGS model for the LES analysis.This model is based on the assumptions of local equilibrium and eddy viscosity and does not require the solution of the transport equations to calculate the SGS stresses.The SGS model also includes model constants whose values must be set in advance.In the Smagorinsky model, a model constant value of about 0.1 is considered a typical value [1].
In this study, an inviscid field is used to investigate the kinetic energy conservation properties from the point of view of kinetic energy conservation properties [7,9].In an inviscid flow field, kinetic energy is analytically conserved.To set up the inviscid flow field, the viscous dissipation and the turbulent viscosity are set to zero from the governing equations.To set the viscous dissipation to zero, the kinematic viscosity coefficient is set to zero.Setting the kinematic viscosity coefficient to zero sets the inviscid limit.For the eddy viscosity, the value of the Smagorinsky model constant is set to zero to set the magnitude of the eddy viscous dissipation in the SGS terms to zero; in the Smagorinsky model, the model constant is originally set in Cs, but in the Smagorinsyk model implemented in OpenFOAM, the value of the model constant is set by two parameters, Ck and Ce.Here Cs and Ck, Ce are related as in this equation. ( In this study, the magnitude of the SGS eddy viscous dissipation is set to zero by setting Ck = 0 and Ce = 0.1.For the inviscid analysis, three spatial resolutions Nx = 10, 32 and 100 are set to investigate the spatial resolution dependence and time evolution up to time step t = 20.
In the viscous condition, the number of spatial grids was set to 32 in each direction.A time step width of 0.005 and a Reynolds number of 1600 were used.The values of the model constants were calculated from 0.1 to 0.6 in increments of 0.1 to calculate the turbulence energy values.In order to find the optimal model constant values for the SGS model, this study also uses the Fourier spectral method as an analysis to obtain reference results: the same flow field conditions as those analysed by the LES are calculated.The values of the model constants of the SGS model are calibrated using the turbulence energy values obtained from this reference analysis and the model constants to be set are calculated.The derived governing equations are developed in time using the six stage fourth order accurate Runge-Kutta method [9].The 2/3 rule is used for aliasing.The spatial resolution for spectral analysis is set to the cube of 256.This spatial resolution is fine enough for direct analysis without modelling.The initial field in the spectral method is the same as in this LES.

Examination Using the Inviscid Field
The time evolution of the kinetic energy of the inviscid flow field with an isotropised TGV for the initial field is shown in figure 1.For each spatial resolution condition, the distribution of the kinetic energy is shown at the top for the initial time t = 0 and at the bottom for the time step t = 20.For condition (a), the kinetic energy decays significantly at time step t = 20.In addition, under conditions (b) and (c), where the spatial resolution is increased, the kinetic energy also appears to decay significantly at time step t = 20.In an inviscid flow field, kinetic energy is inherently conserved.However, in the present OpenFOAM analysis, the kinetic energy appears to decay.As a result, it was confirmed that errors in the conservation of kinetic energy exist in the OpenFOAM analysis.The shape of the kinetic energy distribution at t = 20 retains the initial symmetry observed at all spatial resolutions.Thus, even if the kinetic energy is not conserved, the flow field appears to evolve with time while maintaining symmetry.
Figure 2(a) shows the time evolution of the spatially averaged value of the kinetic energy in the inviscid flow field.As in the previous visualisation, a comparison is made for three spatial resolutions Nx = 10, 32 and 100.In the inviscid field, the time evolution essentially retains the initial value.In the present analysis, the magnitude of the initial kinetic energy is set to unity, so that the analytical solution is shown as the blue dotted line in the diagram.For all spatial resolution conditions, the kinetic energy decays with time: for Nx = 100, the kinetic energy decays more significantly after t = 10 than for smaller spatial grid widths, and after t = 10 there is no spatial grid width dependence in the kinetic energy conservation property.The kinetic energy conservation property shows no spatial grid width dependence after t = 10.However, before t = 10, the kinetic energy increases under conditions with large spatial grid widths and a spatial grid width dependence is observed in the kinetic energy conservation property.Therefore, the focus is on the initial time before t = 10.As shown previously, there is a spatial grid width dependence of the decay of the conservation of kinetic energy at the initial time.The spatial grid width dependence of the initial deviation is shown in figure 2(b).The horizontal axis shows the grid width and the vertical axis shows the magnitude of the slope of the kinetic energy change at t = 0.The slope obtained here is expected to follow the order of accuracy of the spatial discretisation.As shown in figure 2 IOP Publishing doi:10.1088/1742-6596/2694/1/0120045 inviscid analysis is of second-order spatial accuracy.This result suggests that the OpenFOAM analysis is of second-order spatial accuracy.

Examination Using the Viscous TGV Flow
The time evolution of the spatially averaged turbulence energy in the DNS is shown by the solid black line in the figure 3(a).The time evolution was performed from t = 0 ~ 20. Figure 3(a) shows that the mean turbulence energy decreases with time: the slope of the figure is small and the decrease in turbulence energy is moderate up to t = 3.At t = 20 the turbulent kinetic energy finally reaches a value of about 0.2.The red points and dashed lines show the time evolution of the mean kinetic energy in OpenFOAM-LES with a model constant value of 0.1.On the other hand, the slope decreases with time, with OpenFOAM showing larger turbulence energy values than DNS from about time t = 15.Thus, comparing DNS and OpenFOAM-LES with a model constant value of 0.1.In this study, the magnitude of the eddy viscosity is varied by changing the values of the model constants and the expression for the eddy viscosity coefficient, which represents the magnitude of the eddy viscosity in the Smagorinsky model, is as follows.In OpenFOAM-LES the model constant value Cs is determined by setting two parameters, Ck and Ce.In this investigation, the values of the model constant are varied from 0.1 to 0.6 in increments of 0.1.The figure shows the average turbulence energy for each of the varied model constant values.The turbulent kinetic energy is plotted as a function of the model constant value, indicating that the turbulence energy continues to decrease as the model constant value increases.Also, the deviation of the turbulence energy values from the DNS results is more significant for all results.
The dissipation rate in the TGV has been found from previous studies to be the time derivative of the turbulence energy and is defined as follows [10,11]: Figure 3(b) shows the time dependence of the dissipation rate.The black line shows the results for the DNS, where the dissipation rate increases until about t = 4, oscillates from t = 4 to about 9, and then continues to decrease.The dashed lines show the OpenFOAM-LES results for each model constant value, where the dissipation rate in OpenFOAM increases rapidly until about t = 2, after which it decreases.Compared to the DNS results, the OpenFOAM dissipation rate is larger at the beginning, but the DNS dissipation rate is larger from about t = 4. Thus, the difference in decay between OpenFOAM 0 10 20 0.03 and DNS is quantitatively observed.Comparing the results for different values of the model constant, at the initial time the dissipation rate is also larger for larger values of the model constant, but at the developed time the dissipation rate is smaller for larger values of the model constant.Developed decaying turbulence is generally assumed to follow a power law [15,16].The power law equation is expressed as follows when the reference time is t0 and the energy at that time is K0.
Here n is the decay exponent.Differentiating this equation with respect to time gives the following equation.
dK/dt = -nK0 (t -t0) -(n + 1) . ( Furthermore, these equations can be used to derive an expression for turbulence time scale T, which is a linear function with respect to time, for the time region.
In this study, the turbulence time scale equation is used to derive the value of the decay exponent.Since the decay coefficient has disappeared, there is no need to calculate the decay coefficient.
A plot of the time scales is shown in figure 4(a).As shown in the figure, the time scale of the DNS is divergent in the initial time.It decreases overall until about t = 10, after which it increases macroscopically.In OpenFOAM, the time scale values after t = 10 oscillate when the model constant is between 0.1 and 0.4.On the other hand, for 0.5 and 0.6 there is no oscillation and in particular for 0.6 there is a linear increase after t = 12.Therefore, using the DNS results, a linear approximation was performed for the three conditions and at times after t = 12 for model constants 0.1 and 0.6 in OpenFOAM.This resulted in a virtual origin value of 0.495 for the DNS, 3.70 for Cs = 0.1 and 0.0135 for Cs = 0.6.The time scale results show that the decay exponent of the DNS is approximately 1.3, which is a generally accepted value for turbulent phenomena.In OpenFOAM, when the model constant is set to 0.1, the value is 0.67, which is about half the value of the decay exponent in the DNS.On the other hand, when the model constant is set to 0.6, the value is 1.35, which is in good agreement with the DNS results for the decay exponent.The results in figure 3(a) show that a model constant of 0.1 is closest to the DNS results, while a model constant of 0.6 shows the most significant deviation from the DNS results among the OpenFOAM results.However, the decay exponent results show that the turbulent phenomena with Therefore, although the turbulence energy values differ significantly from the DNS results, in terms of decaying characteristics, a model constant value of 0.6, which is a more eddy viscous value than 0.1, which is a common value for the Smagorinsky model of LES, is considered appropriate.

Conclusion
Previous studies of OpenFOAM have used turbulent channel flows rather than TGV flows, which are unsteady, turbulent and globally anisotropic.This study aims to investigate OpenFOAM from an unsteady perspective.The physical conservation laws using isotropic TGV flows are used to validate OpenFOAM with a focus on the energy conservation laws of the inviscid state; the turbulence characteristics of isotropic TGV flows obtained using Large-Eddy Simulation (LES) are compared with spectral methods.Investigate the decaying characteristics of the Taylor-Green vortex using an inviscid field and maintaining energy constancy; compare OpenFOAM LES results with DNS results, also considering model constant values.The governing equations use the continuity equation for incompressible fluids and the Navier-Stokes equations.The computational domain is a cube with 2π per side and periodic boundary conditions.The PISO algorithm and the Smagorinsky model are used with the central difference method as the spatial difference method.The time evolution of the kinetic energy of the inviscid flow field is studied.Using an initial setup of isotropic TGVs, the distribution of kinetic energy was observed under different spatial resolution conditions.Using an initial setup of an isotropic TGV, the distribution of kinetic energy was observed under different spatial resolution conditions.The results show that there are errors in the conservation of kinetic energy in the analysis using OpenFOAM.The results show that there are errors in the conservation of kinetic energy in the analysis using OpenFOAM.The time evolution of the spatially averaged kinetic energy in the inviscid flow field confirmed that the kinetic energy decreased with time for all spatial resolution conditions.The time evolution of the spatially averaged kinetic energy in the inviscid flow field confirms that the kinetic energy decreases with time for all spatial resolution conditions.Furthermore, a comparison of the time evolution of the turbulence energy in DNS and OpenFOAM-LES showed that the decay is different.It was shown that the decay of the turbulent energy increases with increasing values of the model constants.Using a time scale to evaluate the turbulence energy decay characteristics, it was observed that the results for DNS and OpenFOAM were different.In particular, a model constant of 0.6 was shown to be in good agreement with the decay characteristics of DNS.

Figure 1 .
Figure 1.Dependence of the spatial resolution on the time evolution of the kinetic energy in inviscid fields.Here this plane is obtained at z = 0 in the computational domain.

Figure 2 .
Figure 2. (a) Time evolution of the kinetic energy K in inviscid fields.This kinetic energy is spatially averaged at each time point.(b) Time derivative of the kinetic energy at the initial time, where this value should be zero in the inviscid field.
(b), the magnitude of the slope of the kinetic energy decrement is proportional to the square of the spatial grid width.This indicates that the present (b) N x = 32, t = 0 (a) N x = 10, t = 0 (c) N x = 100, N x = 32, t = 20 (d) N x = 10, t = 20 (f) N x = 100, Conference on Mechanical Engineering and Materials Journal of Physics: Conference Series 2694 (2024) 012004

Figure 3 .
Figure 3.Time evolution of the turbulent kinetic energy and its dissipation rate, K and e.These quantities are spatially averaged at each time.

Figure 4 .
Figure 4. (a) Time evolution of the turbulence time scale T. (b) Turbulent kinetic energy K divided by a power-law function based on the decay exponent, t -n .
of 0.1 have different decay characteristics from those of the DNS.When the model constant is 0.6, the decay characteristics are similar to those of the DNS as shown in figure 4(b).