Performance of the Finite Element and Finite Volume Methods for Large Eddy Simulation in Homogeneous Isotropic Turbulence

Department of Mathematics, Shahjalal University of Science & Technology, Sylhet-3114, Bangladesh Institute of Industrial Science, the University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 1538505, Japan Graduate School of Engineering, Hokkaido University, Kita 13 jou Nishi 8 chome, Kita-ku, Sapporo-shi, Hokkaido 060-8628, Japan Department of Mechanical and Aerospace Engineering, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8552, Japan


Introduction
During the past decades, LES has been demonstrated to be an accurate and sophisticated predictive method for flows of engineering interest.The position of LES approach is intermediate between DNS and RANS (Reynolds-Averaged Navier-Stokes) techniques.Although recent development of supercomputers enable to carry out the DNS [1-5], which is considered as the exact approach of turbulence simulations but the grid dependence is very high (proportional to Re 9/4 ) and calculation is fairly time consuming, so that the DNS is not appropriate to the practical use.On the other hand, LES is less expensive and can simulate very complex flow fields in turbulence.Unlike the full-scale turbulence modeling of RANS technique, in LES method, large-scale motion is exactly calculated and the effects of subgrid-scale (SGS) motion is modeled.Although LES being superior to RANS, it still has some theoretical and practical drawbacks; and still today, people are putting their efforts to eliminate this drawback for LES.
However, in recent days the important issues for LES are numerical method and subgrid-scale (SGS) modeling.It is known that the numerical methods which are widely used for LES are either spectral or the conventional finite difference method with structure grids [6][7]; but for the case of complicated flow the use of structure grids method is often unsuitable.Since Finite Element Method (FEM) is based on unstructured grids, this method is very useful for engineering applications to the complicated flow fields.
The objective of our present study is to develop "Front Flow" next generation fluid simulator based on LES using FEM and FVM in order to apply to the engineering and practical problems.However, before applying there, it is necessary to examine the effectiveness and performance of our numerical solver through the benchmark problem which other researchers have examined.
In this study, LES is performed in homogeneous isotropic turbulence using FEM and FVM formulation and the results are compared with DNS result which is calculated by a spectral method.The results with SSM and DSM in both FEM and FVM are also compared.Our interest is in whether any effects of the subgrid-scale model should appropriately be damped out.In both cases, we also discussed about vortical structures in the computed flow fields by LES comparing with those in the instantaneous DNS data by visualization of flows.

Numerical scheme and SGS model
In this study Front Flow numerical solvers based on FEM and FVM are used for the computations.The detailed mathematical formulation and developing procedure of these numerical solvers are too long.Since the purpose of this study is to show the performance and effectiveness of these two codes in a benchmark problem so the mathematical formulations are not given here.However, a brief description of these Front Flow numerical codes is given as follows:

Front Flow/FEM code
The Front Flow/FEM numerical code used in this study is general-purpose fluid simulation code which calculates incompressible unsteady flows in arbitrary shaped geometry that involves moving (but not deforming) boundaries.It is particularly designed for computing unsteady flows on turbomachinery and simulating sound pressure spectra that result from unsteady fluid motion.In order to obtain an accurate sound spectra in general, it is of up most importance to simulate the source, i.e., the fluid fluctuations accurately in terms of their spatial and frequency spectrum.
This code is based on an explicit streamline-upwind finite element method with second order accuracy both in time and space.In this numerical scheme, the spatial discretization is performed by an explicit hexahedral FEM and the coordinate system is the three dimensional Cartesian.
The pressure algorithm is based on the ABMAC (Arbitrary Boundary Marker and Cell) method proposed by Viecelli [8], where both the velocity components and static pressure are simultaneously corrected until the maximum divergence of the flow field decreases to less than a prescribed critical value.The subgrid-scale model for LES supported in this code are the SSM [9] incorporating the wall-damping function and the DSM proposed by Germano et al. [10].The details of the numerical methods and developing procedures of this code have been given in the previous studies by Uddin et al. [11] and Kato et al. [12].

Front Flow/FVM code
In this code, the spatial discretization is performed by a finite volume method and the coordinate system is the three-dimensional Cartesian.The usable elements in this code include hexahedrons, triangular prisms, square pyramids, and tetrahedrons, however, in this study we use hexahedrons elements only.In this numerical solver the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) scheme is adopted for computations of compressible flows and a formula based on the low Mach number approximation can be used, however, using this code the incompressible flow computation is also possible.Like FEM code, in this solver the adopted SGS models include the SSM and DSM for SGS modeling.The detailed numerical methods and developing procedures of this code have been given in the previous study by Unemura et al. [13].

The computation
The computational domain of the fine mesh was selected to a periodic box (2π×2π×2π) and the computation was performed using 64 3 grids.The present FEM and FVM codes assume the solution at t = nΔt is known, where Δt is the time step increment, t is the nondimensional time and n is the time step, and the solution at the next time step, t = (n+1)Δt, can be calculated using the residual terms described in the finite element and finite volume formulation [11][12][13].Since we are interested to compare our present results with the result of spectral DNS, the initial flow field for LES is calculated with same condition and procedure that is done for spectral DNS by Tanahashi et al. [2][3].The calculation presented here was done with nondimensional Δt = 0.00316 and a ν = 1/Re = 8.26×10 -3 .In the case of SSM, the Smagorinsky coefficient, C s = 0.2 is used in both FEM and FVM simulations [14].

DNS database
The reference DNS is performed at 64 3 resolutions by using a spectral code developed by Tanahashi et al. [2][3].At the end of calculation, the Taylor microscale Reynolds number, Re λ , based on root mean square of velocity fluctuations (u rms ) and Taylor microscale (λ) of the DNS data is Re λ = 30.5 (t =3.792), and the maximum possible Reynolds number (Re) of the flow is 121.1.

Numerical Results
In this section we compare the numerical results in LES with the results of DNS to understand the decay and statistical behavior of homogeneous isotropic turbulence.We will also discuss about vortical structures in LES data comparing with DNS data.

Decay of turbulence
The comparisons of the three-dimensional energy spectra in DNS and LES data at the end of calculation (t = 3.792) for both FEM and FVM with SSM and DSM models are presented in Fig. 1.The energy spectrum is calculated by the definition given as: In Fig. 1, we can observe that the DNS spectrum shows the power decay close to k -5/3 .The energy spectrum in LES calculation using FEM with DSM model shows almost similar decay as in DNS calculation at the whole wave number range (in Fig. 1 LES but the DSM result suggests that the decay of turbulence in LES follows the k -5/3 power law, and the numerical accuracy is quite good.On the other hand, although the decay of turbulence using FVM collapse with DNS data in the low wave number range but it underestimates DNS data in the high wave number range for both SSM and DSM models (in Fig. 1b).Moreover, the SSM results show little faster decay than DSM results.The difference between SSM and DSM results may happen due to Smagorinsky constant or model itself.However, in both FEM and FVM codes, it reveals that DSM model works better than SSM model comparing with DNS data (in Fig. 1c, d).It is also revealed that among these two codes, the performance of FEM code is better than FVM code in these simulations.

Turbulence statistics
In this sub-section we discuss about the turbulent statistics in LES comparing with DNS results.Fig. 2 shows the decay of total resolved energy, E k , where 2 2 It reveals that FEM results with DSM model coincide with the DNS through out the entire analysis (Fig. 2(b)) and SSM result is dissipative in the early stage but collapse with DNS data at the end of the calculation (Fig. 2(a)).On the other hand, FVM results do not collapse with DNS results through out the entire analysis for both SSM and DSM.
The decay of resolved enstrophy is presented in Fig. 3, where enstrophy is defined as: In this case, the FVM results with SSM and DSM models differ significantly from the DNS through out the entire analysis and show too much dissipation in the initial stage of the calculation.The FEM results with SSM model are also dissipative in the initial stage and it shows relatively higher dissipation from DNS throughout the entire analysis (Fig. 3a).However, FEM results with DSM are in good agreement with DNS (Fig. 3b).However, the performances of these two codes with DSM is better than that of SSM.The decay of root mean square of velocity fluctuations (u rms ) is presented in Fig. 4, where, The decay of u rms for FEM with DSM fully collapses with DNS data in the whole analysis.Although, initially FEM results with SSM is little faster but it is in good agreement with DNS after time t=2.5.On the other hand, FVM results with SSM and DSM do not collapse with DNS through out the entire analysis.Hence it is revealed that the FEM-LES result gives good agreement with DNS result in terms of both spatial spectra and decay of turbulence statistics [11].Skewness and flatness factors of velocity are important statistical properties to rep r resent characteristics of turbulence.The production of the rate of dissipation of turbulent kinetic energy or, equivalently, the production of enstrophy is directly related to skewness in isotropic turbulence [15].Skewness and Flatness of a random variable are presented in Fig. 5, and Fig. 6, respectively.Skewness and Flatness of random variables, respectively, are defined as: The agreement of skewness (in Fig. 5) of LES results with DNS data in both FEM and FVM simulations for both SSM and DSM cases is good until t =2.0 and it overestimate DNS data hereafter.In this case, particularly it is revealed that at the end of calculation the profile of skewness for FVM code with DSM goes very close to the DNS results.However, skewness for all cases is almost zero at t =0 and the LES results for both FEM and FVM solvers either with SSM or DSM agree well with DNS data.On the other hand, flatness (in Fig. 6) of LES results for both FEM and FVM simulations collapse with DNS data well throughout the entire analysis.At time t =0, flatness of DNS and LES for all cases is almost 3.0 and shows nearly a constant value at the end of the calculation.The behavior of flatness suggests that turbulence velocity at the end of the calculation does not include the effects of initial condition and reaches to fully developed state.

Visualization of flows
In this section we shall discuss about the vortical structures in LES comparing with DNS result by visualization of flows.There are several methods for identification of vortical structures and their visualization in turbulence with significant differences [16] and most of them show threshold dependence.In this study, by direct use of 'local flow pattern method' we discuss about vortical structures in LES flow fields.From the distributions of second and third invariant of the velocity gradient tensor one can define 'local flow pattern method' which does not depend on the thresholds of the variables.Details of this method can be seen in the previous studies by Tanahashi et al. [2][3].
The second invariant of the velocity gradient tensor is defined as: ( ) where, 1 2 are the symmetric and asymmetric part of the velocity gradient tensor, Fig. 7 shows the contour surfaces of the positive second invariant given in Eqn.(6) of the velocity gradient tensor in DNS and LES data for FEM and FVM simulations as well as for SSM and DSM cases at the end of analysis (t = 3.792).In this figure the visualized region is the whole calculation domain and the view point is same for all cases.The level of the isosurface is selected to be Q =10 for all cases.These figures show that lots of tubelike vortical structures are randomly oriented in DNS data as well as in LES data for all cases.In this visualization we considered Q =10 only to show the vortical structures in DNS and LES.However, if we increase or decrease the value of Q, we can also show distinct tube-like structures exists in DNS and LES, of course, it will be little bit different from the present visualized structures.Our present study reveals that in actual LES we can have the coherent tube-like structures somewhat similar to the structures in DNS and these structures are quite unique and distinct as we can see in Fig. 7.For FEM code, appearance of the vortical structures using DSM case seems higher than that of SSM case and close to DNS data.On the other hand, for FVM code, appearance of these vortical structures for SSM and DSM cases are similar but differ from the DNS, even it decrease from the SSM case in FEM code.This observation again suggests that the accuracy of LES calculation using FEM with DSM is better than that of FVM. this study it is revealed that the coherent structures in FEM-LES are similar to the structures in the DNS data which suggest that FEM-LES results may have suitable role to tackle turbulence intermittency like DNS.

Conclusions
Larg eddy simulations in homogeneous isotropic turbulence have been done using FEM e and FVM formulation in order to show the performance of these numerical methods as well to assess its spectral accuracy.The results in both cases are compared with those from DNS based on a spectral method.It is shown that the results of FEM simulation are in better agreement with DNS than the results of FVM simulation.Regarding SGS modeling it revealed that the performance of DSM is better than SSM in both FEM and FVM formulation.In FEM formulation DSM gives good agreement with DNS results in terms of both spatial spectra and decay of turbulence statistics.
Visualization of second invariant of velocity gradient tenso d by LES for both FEM and FVM simulations reveal the existence of distinct, coherent and tube-like vortical structures similar to those found in instantaneous flow field computed by the DNS.In this visualization it is also revealed that the appearance of the vortical structures in FVM simulation is lower than that of FEM simulation for both DSM and SSM modeling.Moreover, the appearance of these vortical structures using DSM in the resolved turbulence by FEM seems higher than that of SSM, and very close to DNS data.

A
T and Technology (MEXT), Japan under an IT research program "Frontier Simulation Software for Industrial Science".

Fig. 1 .Fig. 2 .
Fig. 1.Comparisons of three-dimensional energy spectra of velocity fluctuations in DNS and LES data.(a) SSM and DSM results using FEM; (b) SSM and DSM results using FVM; (c) SSM results using FEM and FVM, (d) DSM results using FEM and FVM.The legend FEM and FVM stand for finite element method and finite volume method, respectively.
Development of skewness of a velocity component.(a) SSM and (b) DSM results using FEM F and FVM, respectively. of flatness of a velocity component.(a) SSM and (b) DSM results using FEM and FVM, respectively.