Large eddy simulation of anisotropic turbulence investigating thermal field effects

Abstract

Keywords: Large eddy simulation; return to isotropy of turbulence; pressure scalar gradient correlation; turbulent heat flux modeling

1. Introduction

The return to isotropy of anisotropic turbulence is a fundamental phenomenon studied extensively to analyze the structure of turbulent flows [1,2,3]. It also represents a key step in understanding the complexity of homogeneous turbulent flows and serves as a cornerstone for calibrating models of the finite part of the pressure-strain correlation and the pressure-scalar gradient correlation [4,5,6]. The concept of return to isotropy in turbulence refers to the process by which a turbulent flow, initially exhibiting anisotropic behavior, tends to become isotropic over time, particularly at larger scales or in the absence of significant directional influences.

Within the realm of thermal fluid dynamics, the isotropization process governs the efficiency of numerous engineering applications. In turbulent convection—whether driven by buoyancy or external forcing—anisotropy typically emerges from sharp temperature gradients. As these flows evolve, the progression toward isotropy offers a pathway for simplifying heat transfer parameterizations. This return to isotropy is a bottleneck for better heat transport forecasts in boundary layers and for scaling industrial exchangers or atmospheric models [7]. Effective thermal management in high-gradient systems such as furnaces and combustion reactors, depends on isotropization dynamics. Similarly, in geophysical fluids, buoyancy and stratification forces dictate how eddies recover isotropy—a fundamental requirement for reliable climate mapping and

Corresponding Author E-mail Adress: mounir.bouzaiane@fsb.ucar.tn *

Submitted: 28 December 2024 ; Accepted: 12 May 2025 This paper was recommended for publication in revised form by Editor-in-Chief Ahmet Selim Dalkılıç Published by Yıldız Technical University, İstanbul, Türkiye This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).

weather forecasting [8]. This scalability is also vital for environmental engineering, where plume isotropy determines the accuracy of pollutant dispersion and heat mitigation strategies. The work of Rotta in 1951 was the first one [9], which introduced a linear relaxation framework to describe Reynolds stress evolution in homogeneous turbulence. This established theoretical scaffolding for decades of experimental and numerical validation. More than a decade later in (1966) Comte-Bellot and Corrsin utilized grid-generated turbulence to confirm that anisotropic states naturally decay toward isotropy when mean velocity gradients are absent [10]. By (1980), Gence and Mathieu investigated the impact of successive plane strains, revealing that the isotropization rate depends heavily on the directional distribution of turbulent energy, noting that specific anisotropic configurations may inhibit the process [11]. Choi and Lumley (2001) further advanced these findings by examining turbulence subjected to plane distortion and axisymmetric expansion [12].

are the pressure strain correlation term and dissipation term respectively

The resuls of direct numerical simulation from Iida and Kasagi [13] in the case of the return to isotropy of homogenoues turbulence imposed on a themal field with a temperature gradient can be considered the principal results during the last two decades. Our study builds on this by using Large Eddy Simulation (LES) [1416] to investigate the same fundamental flow. We specifically aim to fill the gap in LES-based exploration of the effects of thermal fields on isotropization. While most models prioritize the velocity field, our work treats thermal gradients as a primary driver of the transition. By resolving large-scale structures and modeling sub-grid scales, we provide a more granular view of how thermal dynamics drive the return to isotropy. LES which serves as an intermediate approach between DNS and Reynolds-Averaged Navier-Stokes (RANS) simulations [17-21], allows for the deterministic simulation of large-scale motions [22], while statistically modeling the influence of small-scale turbulence. This methodological choice motivates the structure of this work [23], which is organized as follows: The next section introduces the fundamental Reynolds-Averaged Navier-Stokes equations and the governing equations of the LES method, along with the computational conditions. The third section presents a comparison with the DNS results obtained by Iida and Kasagi. Finally, in Section 4, after a brief review of our previous work on a new formulation of a pressure–scalar gradient correlation model based on realizability, an optimization of the model constants is carried out numerically.

Figure 1. Coordinate system and imposed temperature gradient With the same flow field, a constant mean temperature gradient as shown in figure 1, which would create a turbulent heat flux was imposed in the y-direction (or x2 direction) Then the budget equation of the turbulent heat flux is derived as follows: 2ui i = Pii + { ii - f ii (4) 2t

are production, scalar pressure gradient correlation and dissipation terms respectively.

3. Results and discussions

Figure 3 illustrates the time evolution of the Reynolds stress tensor components, with solid lines representing our simulation results and circles denoting the DNS results of Iida and Kasagi. A qualitative agreement is observed between the two sets of results. The Reynolds stress components decrease with increasing time from t=0 to t=5. However, our simulation tends to overestimate values relative to the DNS results reported by Iida and Kasagi.

The spectral equation associated with the filtered temperature field, as discussed in [14, 26], is written as follows: ti 2it 2 t ij ) + vr e k 2 it = -ut j D ij + Ik j (-u j it + v' e 2 + v' e D xj 2t

r o e - o) ol o (K 2T L v e = Pr + Prt , vl e = Prte and K ij = correspond to the 2x j mean thermal diffusivity, fluctuating thermal diffusivity, and the mean temperature gradient, respectively. (Pr and Prt represent Prandtl number and turbulent Prandtl number respectively) where ν =0.00294 and Pr =0.7 In previous works, Schiestel [24, 25] developed a Large Eddy Simulation (LES) method based on a pseudo-spectral approach to study kinematic homogeneous turbulence with initial anisotropy. Building upon this framework, the present study [24] extends the code to account for a passive scalar introduced by a constant mean temperature gradient in the x2 direction. The flow remains kinematic and is initially anisotropic. This extension enables a more detailed investigation of the coupling between thermal and velocity fields, offering insights into scalar transport mechanisms driven by turbu3 lence. Numerical simulations are carried out on a 64 grid (Figure 2 – Annex 1), within a cubic domain of side length L=2π, with pe-

Figure 4 presents the time evolution of the principal component of the pressure-strain { 22 correlation (the redistribution term). The same legend as in Figure 3 is used to ensure consistency. Good qualitative agreement is observed between our results and the DNS of Iida and Kasagi. However, the quantitative agreement is only evident for times t greater than 2. For times t less than 2, our results underestimate the evolution predicted by the DNS. This discrepancy suggests that the Smagorinsky model employed in this study to model the subgrid-scale terms is unable to fully capture the influence of low-Reynolds-number effects, which are more clearly observed in the DNS results. The classical Smagorinsky model does not account for fluid viscosity, which limits its accuracy at low Reynolds numbers.

Figure 6 presents the time evolution of the redistribution term in the pressure-temperature - gradient correlation. Excellent agreement is observed between our results and the DNS results of Iida and Kasogi for times greater than 2.5. However, for t less than 2.0, only qualitative agreement is observed. Both sets of results exhibit an initial decay phase, with the simulations showing similar trends up to t=1.

Figure 4. Time Evolution of U 22 We now turn our attention to the thermal field. Figures 5, 6, 7, 8, and 9 illustrate the time evolution of the turbulent fluxes and their associated budgets. Figure 7. Time evolution of production term In Figure 7, we show the time evolution of the production term P2i . Again, we observe excellent quantitative and qualitative agreement between our results and those reported by Iida and Kasagi for t greater than 2. However, for times less than 2, our results slightly underestimate the DNS results. Figure 5. Time Evolution of u 2 i

Figure 5 presents the time evolution of the turbulent flux. Our LES results show strong agreement with the DNS data from Iida and Kasagi, particularly at later times. The asymptotically constant values reached in our simulations confirm the trends predicted in their study. However, for early times (t<2), noticeable discrepancies appear, reflecting the limitations of the classical Smagorinsky subgrid-scale model in capturing low Reynolds number effects. During this phase, the flow remains highly anisotropic, and the model fails to reproduce the DNS values with sufficient accuracy—this constitutes one of its main shortcomings. As time advances (t>2), the flow becomes increasingly isotropic, and the agreement between LES and DNS improves significantly.

Figure 8 presents the time evolution of the dissipation term f 2i . Only qualitative agreement is observed between our results and those of Iida and Kasagi. Both simulations show an initial decay phase. However, in the DNS results, there is a slow increase from t=1 to t=2.5 whereas our results quickly reach asymptotic equilibrium values starting from t=1.5

Figure 8. Time evolution of the dissipation term The results indicate that for t<2, the flow remains strongly anisotropic, and the Smagorinsky model fails to accurately reproduce the DNS results, highlighting its main limitation. As time progresses (t>2), the turbulence becomes increasingly isotropic, leading to improved agreement between LES and DNS. Figure 6. Time evolution of the temperature-gradient correlation term

This temporal evolution illustrates the gradual return to isotropy. Anisotropic features inherited from the initial conditions dominate the early flow. With time, nonlinear interactions and energy

redistribution mechanisms promote the decay of these structures, ultimately resulting in a statistically isotropic state. This behavior is consistent with classical theories of return to isotropy, such as Rotta’s model and rapid-distortion turbulence frameworks. Quantitative analysis confirms this trend, with decay rates and spectral energy redistribution aligning with theoretical predictions. A key limitation of the subgrid-scale (SGS)model, particularly at low Reynolds numbers, lies in its inadequate representation of small-scale turbulence dynamics. Due to limited grid resolution, the model often assumes isotropy in regions where anisotropic effects prevail, resulting in inaccurate predictions of turbulent stresses. Moreover, the closure assumptions in SGS models may fail to capture the complex interactions between large and small eddies, especially in structured, weakly dissipative turbulence typical of low Reynolds number flows.

4. Calibrations of models constants

The results of the large eddy simulations (LES) are examined in the context of a new calibration of the model coefficients associated with the return-to-isotropy term in the pressure-scalar gradient correlation. In a previous study, Bouzaiane and Lili [27], as well as Gaaliche et al. [28], proposed and tested a nonlinear model for the return-to-isotropy term. At this stage of our work, we will briefly review key aspects of these prior contributions. We have developed a model for the return-to-isotropy of the pressure-scalar gradient correlation in the following form [26]: { ii = -

This model can be considered an extension of the return-to-isotropy model for the pressure-strain correlation, originally proposed by Sarkar and Speziale [6]: II { ij = -f (a 1 b ij - a 2 (b ik b kj - 3b d ij))

The realizability condition for the thermal field is written as: dD 11 dt $ 0

After a straightforward development, we arrive at the final condition that our coefficient model must satisfy: C 2 4 3 C 1 - r - 1 + 3 C 2 + 9 C 3 - (C 2 + 33 + 2 (a 1 - 1)) I D + II 1

1. II 27II

Here, II b, II D are scalar invariants derived from tensors b and D, and r represents the ratio of the mechanical to the scalar time scale.

The optimization process for the model coefficients C1, C2, C3 subject to the realizability condition, involves the integration of the heat flux transport equation using a finite difference scheme: 2ui i = Pii + { ii - f ii 2t

In this equation, our model is introduced, and the defined criterion is to find the coefficients C1, C2, C3 that minimize the difference err = R ( u i i N - u i i LES) between the values obtained from the numerical integration of the equation and those from the LES method

After an extensive optimization process, the final values for the model coefficients are obtained as: C1=7.75, C2=−3,C3=−3.

These values also satisfy the strong form of the realizability condition.

where the coefficients a 1 = 3.4 and a 2 = 3 (a 1 - 2) = 4.2 are derived from a stability analysis of the fixed points and the strong form of the realizability condition for the kinematic field, which is expressed as follows: du 21 dt $ 0

The concept of realizability was initially introduced by J. Lumley [1] and was used by Sarkar and Speziale [6] to construct their model. In a similar manner, we applied the strong form of the realizability condition for the thermal field in previous work [26, 27], as introduced by J. Lumley, through the positive tensor D [1], where its components are given by: D ij = i 2 u i u j - iu i iu j

Figure 9: Time evolution of turbulent heat flux components, comparing LES with the optimized model Figures 9 and 10 illustrate the time evolution of the three components of the turbulent heat flux (i=1,2, 3) and the pressure-scalar gradient correlation (i=1,2, 3). The symbols represent LES predictions, while the lines correspond to results from the models with optimized coefficients. For the turbulent heat flux, excellent agreement is observed for the dominant component, which aligns with the temperature gradient (dT/dx2​). The remaining two components of the turbulent heat flux are weak, showing only qualitative agreement.

For the components of the pressure-scalar gradient correlation, the model with the new optimization provides a good approximation of the LES values.

Figure 10. Time evolution of pressure-scalar gradient correlation, comparing LES with the optimized

In this work, we employed large eddy simulation techniques to investigate the return to isotropy of turbulence in the presence of a thermal field subjected to a constant mean temperature gradient. The numerical code developed by R. Schiestel was extended to capture thermal quantities such as heat flux, the pressure-scalar gradient correlation, and the dissipation term. The results were compared with those from direct numerical simulations by Iida and Kasagi, and good overall agreement was found for most kinematic and thermal turbulence parameters. Divergences were observed for time evolutions below 2.5, primarily in the components of heat flux and in those of the pressure-scalar-gradient correlation. We believe that a finer mesh of the order of 128³ or more could improve our results. Additionally, incorporating corrections to the Sub-CS) Scale (SGS) models such as the Dynamic Smagorinsky Model, which automatically adjusts the Smagorinsky coefficient (CS)—could significantly improve the results and represents a promising direction for future work. Large Eddy Simulation was also used to optimize the model constants of the return-to-isotropy term in the pressure–scalar gradient correlation model, subject to strong realizability constraints. The optimized constants (C1, C2, and C3) minimized deviations between LES results and model predictions. Moreover, the evolution of the second and third invariants remained within Lumley’s realizability triangle, confirming the physical consistency of the approach.

IIID Figure 12. Time Evolution of the second invariant IID on the third invariant IIID The choice of the coefficients C1 , C2 and C3 is now justified based on the realizability condition. Lumley introduced the realizability triangle for any positive-definite tensor (such as b and D — see figure 11 Annex 2), in which the second invariant is plotted against the third invariant. According to Lumley [1] , any physically realizable turbulence solution must produce values that lie within this triangle. The boundaries of this triangle correspond to limiting cases, namely axisymmetric turbulence and two-dimensional turbulence. Figure 12 presents the second invariant of tensor D (tr (D.D)) plotted against the third invariant (i.e., tr (D.D.D)). It is evident from this figure that the selected values of the coefficients satisfy the strong form of the realizability condition, as they lead to predictions that remain within the bounds of the Lumley triangle.

5. Conclusion

The return to isotropy in homogeneous turbulence is widely regarded as a canonical flow that provides a framework for studying various turbulence mechanisms. This phenomenon is relevant to numerous applications, including thermal convection, meteorology, oceanography, and various industrial and engineering processes.

This work enhances the accuracy of turbulence models, particularly in thermal flows, and supports the optimization of turbulent flow simulations in practical engineering applications, such as heat transfer in combustion chambers and climate modeling. Future research will explore homogeneous sheared turbulence, rotating flows, and stratified turbulence to extend the proposed model to more complex turbulent systems.

Acknowledgement

We are grateful to Mr. R. Schiestel for providing the LES code developed at IRPHE (Institut des Phénomènes Hors Équilibre, Marseille), without which this study would not have been possible.

Nomenclature

b ij = ( u i u j /2K) - (d ij /3) Reynolds stress anisotropy tensor 2 2 K Turbulent kinetic energy [m /s ] P Production due to mean kinetic and scalar gradients 2 p Pressure[N/m ] t Time [s] T Temperature [°C] ui i-th component of the fluctuating velocity [m/s] i-th component of the mean velocity [m/s] U i 2 2 Reynolds stress tensor [m /s ] u i u j Turbulent scalar flux [m°C /s] u i i

Thermal diffusivity of fluid [m /s] Kronecker symbol 2 3 Dissipation due to molecular effects [m /s ] Fluctuation of the scalar temperature [°C] 2 Temperature variance [°C ] -1 -1 Dynamic viscosity [kg.m .s ] 2 Kinematic viscosity [m /s] -3 Density of the fluid [kg.m ] 2 3 Pressure-strain correlation [m /s ] 2 Pressure-scalar gradient correlation [°C m/s ]

Lumley, J. L. Computational modeling of turbulent flows. Advances in Applied Mechanics, vol. 18, pp. 123–176, 1978. https://doi.org/10.1016/S0065-2156(08)70266-7 Choi, K. S., & Lumley, J. L. The return to isotropy of homogeneous turbulence. Journal of Fluid Mechanics, vol. 436, pp. 59–84, 2001 https://doi.org/10.1017/S002211200100399X Crespo, M., Jamme, S., & Chassaing, P. On the return to isotropy of compressible anisotropic turbulence. École Nationale Supérieure d’Ingénieurs de Constructions Aéronautiques(ENSICA),Toulouse, France.2005 https://doi. org/10.2514/6.2005-5026, Attou, Y., Bouhafs, M., & Feddal, A. Numerical analysis of turbulent flow and heat transfer enhancement using V-shaped grooves mounted on the rotary kiln’s outer walls. Journal of Thermal Engineering, vol. 10 no. (2), pp. 350–359., 2024 https://doi.org/10.18186/thermal.1448621 Yang, P. F., Pumir, A., & Xu, H. Return to isotropy of homogeneous shear-released turbulence. Physical Review Fluids, vol. 6, 044601, 2021 https://doi.org/10.1103/PhysRevFluids.6.044601 Sarkar, S., & Speziale, C. G. A simple nonlinear model for the return to isotropy in turbulence. Physics of Fluids A: Fluid Dynamics, vol. 2 no.(1), pp. 84–93. 1990 https://doi. org/10.1063/1.857686 Akaje, T. W., Taiwo, M. A., Olajuwon, B. I., Raj, M. T., & Akinleye, S. Double-diffusive nonlinear buoyancy force significance on free convective chemically reacting fluids flow past vertical Riga surface Journal of advanced research in fluid Mechanics and thermal sciences vol. 102 no.(1), pp.59–72 2023. https://doi.org/10.37934/arfmts.105.1.5975 Akaje, T. W., Olajuwon, B. I., & Raji, M. T. Computational analysis of the heat and mass transfer in a Casson nanofluid with a variable inclined magnetic field. Sigma J Eng Nat Sci, vol.41 no.(3), pp.512−523. 2023 https://doi.org/10.14744/sigma.2023.00057 Rotta, J. C.. Statistische Theorie nichthomogener Turbulenz. Zeitschrift für Physik vol. 129, pp.547–572. 1951 https://doi. org/10.1007/BF01330059

Comte-Bellot, G., & Corrsin, S. The use of a contraction to improve the isotropy of grid-generated turbulence. Journal of Fluid Mechanics, vol. 25, pp. 657–682. 1966 https://doi. org/10.1017/S0022112066000338 Gence, J. N., & Mathieu, J. The return to isotropy of a homogeneous turbulence having been submitted to two successive plane strains. Journal of Fluid Mechanics, vol 101 no. 3, pp. 555–566, 1980 https://doi.org/10.1017/S0022112080001796 Iida, O., & Kasagi, N. Redistribution of the Reynolds stresses and destruction of the turbulent heat flux in homogeneous decaying turbulence. 9th Symposium on Turbulent Shear Flows, Kyoto, Japan, Paper 24-4-1 to 24-4-6.1993 Lesieur, M., Métais, O., & Comte, P. Large-Eddy Simulations of Turbulence. Cambridge University Press. 2010 https://doi. org/10.1017/CBO9780511781254 Métais, O. Large-Eddy Simulations of Turbulence. In: New Trends in Turbulence, Springer-Verlag, pp. 112–186.2001, https://doi.org/10.1007/3-540-45691-8_3 Dianat, M., Yang, Z., McGuirk, J. J. Large Eddy Simulation of Scalar Mixing. Engineering Turbulence Modelling and Experiments, vol. 6, pp. 823–832. 2005 https://doi.org/10.1016/ S1479-4408(05)80087-0 Madi Arous, F. Numerical simulation with a Reynolds stress turbulence model of flow and heat transfer in rectangular cavities with different aspect ratios. Journal of Thermal Science and Technology, vol. 11 no.(1) ,JTST0012. 2016 https://doi. org/10.1299/jtst.2016jtst0012 Naifer, S., & Bouzaiane, M. Effect of rotation and stratification on homogeneous sheared turbulence through second-order model. Turbulence, Heat and Mass Transfer, 9, 1069–1072. 2018 https://doi.org/10.1615/THMT-18.1220 Chebbi, B., & Bouzaiane, M. On the effects of rotation on the passive scalar and kinematic fields of homogeneous sheared turbulence. Journal of Applied Fluid Mechanics, vol. no. 2, pp 55–67. 2012 https://doi.org/10.36884/JAFM.5.02.12168 Naifer, S., & Bouzaiane, M.. On the turbulent Prandtl number in stably stratified turbulence by second-order models. Journal of Thermal Engineering, vol. 6 no 3, pp.369–380. 2020 https://doi.org/10.18186/thermal.733250 Thamri, L., Naffouti, L., & Bouzaiane, M. A study of the asymptotic equilibrium behavior in stratified turbulence submitted to horizontal shear. Journal of Applied Fluid Mechanics, vol. 7 no.3, pp. 401–413, 2014 https://doi.org/10.36884/ jafm.7.03.20588 Kumar, R., & Singh, N. K. Large Eddy Simulation of Flow Over Elliptic Cylinder Array in Square Configuration at Subcritical Reynolds Numbers. Journal of Thermal Engineering, vol. 7 no.1 pp.,204–219. 2021 https://doi.org/10.18186/thermal.850286 Ben Hamida, L., & Bouzaiane, M. Simulations aux grandes échelles du retour à l’isotropie d’une turbulence homogène en présence d’un champ thermique. JITH 2024, Paris, 29–31 octobre 2024.2024

Schiestel, R. . Les écoulements turbulents : Modélisation et simulation, 2e éd. revue et augmentée, Paris, Hermès Science Publications, ISBN : 2-86601-681-5.1998 Schiestel, R. . Sur la modélisation des écoulements turbulents en non-équilibre spectral. Comptes Rendus de l’Académie des Sciences. Série II, vol.302 no.11, pp. 327–330.1986 Lilly, D. K. The representation of small-scale turbulence in numerical simulation experiments. In Proceedings of the IBM Scientific Computing Symposium on Environmental Sciences, 195–206.1967 Bouzaiane M, Lili T. Réalisabilité d’un schéma non linéaire de retour à l’isotropie de la corrélation pression-gradient d’un scalaire. Rev. Gen. Therm., vol. 37, pp 371–383. 1998 Gaaliche S, Bouzaiane M, Lili T. Réalisabilité d›un modèle non-linéaire de la corrélation pression-température. Mécanique & Industries, vol.7, pp. 565–572. 2006https://doi. org/10.1051/meca:2006027

Figure 11. Triangle of realizability of Lumley Here the limits of the triangle correspond to the limit case of axisymmetric turbulence and bi-dimensional turbulence

Appendice 1 Discretization of the cubic domain with periodic boundary conditions

Figure 11. Triangle of realizability of Lumley Here the limits of the triangle correspond to the limit case of axisymmetric turbulence and bi-dimensional turbulence Figure 2. Discretization of the cubic domain with periodic boundary conditions Periodic boundaries conditions mean that at the edges of the domain (i.e., the faces of the cube), the values of the variable (the velocity) are replicated. This emplies that : U(x+L,y;z)=U(x,y,z) U(x,y+L,z)=U(x,y,z) U(x,y,z+L)=U(x,y,z) To match the conditions of Iida and Kasagi, the dimensionless initial parameters were set as follows: • Kinematic viscosity: ν=0.00294 • Kinetic energy : E=0.29 • Energy peak: kp​=4.12

Figure Captions Figure 1 Coordinate system and imposed temperature gradient Figure 2. Discretization of the cubic domain with periodic boundary conditions Figure 3. Time Evolution of u 22

Figure 6. Time evolution of the temperature-gradient correlation term Figure 7. Time evolution of the production term Figure 8. Time evolution of the dissipation term

Figure 9. Time evolution of turbulent heat flux components, comparing LES with the optimized model Figure 10. Time evolution of pressure-scalar gradient correlation, comparing LES with the optimized Figure 11. Triangle of realizability of Lumley Figure 12. Time Evolution of the second invariant IID on the third invariant IIID

Share and Cite